Using ODE to solve a system of 19 equations. There seems to be something wrong with the system.

3 views (last 30 days)
I am trying to solve a system of 19 differential equations. It is a chemical reaction problem, and each of the 19 equations represent the change in concentration of that species as a function of time, dC/dT.
When I use ODE45 to solve the system, I get unexpected results. The concentration of two of the components become abnormally high. Usually, high concentration of a component favors the reactions that consume the component. The current results are showing the two components that obtain abnormally high concentration. I don't think this should happen.
There is one thing I noticed, and I hope that you can explain. When I change the value of k1 to something absurdly large (i.e. 100000), the results do not change.
function Yv = runODE()
tic;
%Initialize concentrations.
initialvalues = [0.12;3.75;0;3.29;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
timespan = [0 400];
[tv,Yv]=ode45('funsys',timespan,initialvalues);
toc;
end
function Fv=funsys(t,Y)
% Initialize rate constants
k1 = 4;
k2 = 1;
k3 = 0.6;
k4 = 0.4;
k5 = 0.01;
k6 = 0.0008;
k7 = 0.001;
k8 = 0.002;
k9 = 0.0001;
k10 = 5;
k11 = 0.1;
k12 = 0.9;
k13 = 0.1;
k14 = 100;
k15 = 0.05;
k16 = 100;
k17 = 0.1;
k18 = 0.01;
k19 = 0.001;
k20 = 500;
k21 = 100;
k22 = 10;
k23 = 100;
Fv(1) = -k1*Y(1)*Y(2)+k2*Y(3)+k6*Y(6)*Y(4)+k8*Y(5)-k10*Y(1)^2+k11*Y(11)-k14*Y(12)*Y(1)-k16*Y(11)*Y(1)+k17*Y(16)+k24*Y(19)*Y(13);%%%
Fv(2) = -k1*Y(1)*Y(2)+k2*Y(3)-k4*Y(3)*Y(2)+k5*Y(6)-k15*Y(13)*Y(2)^2-k18*Y(15)*Y(2)^2-k19*Y(14)*Y(2);
Fv(3) = k1*Y(1)*Y(2)-k2*Y(3)-k3*Y(3)*Y(4)-k4*Y(3)*Y(2)+k5*Y(6)+k7*Y(5);
Fv(4) = -k3*Y(3)*Y(4)-k6*Y(6)*Y(4)+k7*Y(5)-k20*Y(18)*Y(4)+k21*Y(8)*Y(19);
Fv(5) = k3*Y(3)*Y(4)-k7*Y(5)-k8*Y(5)-k9*Y(5);
Fv(6) = k4*Y(3)*Y(2)-k5*Y(6)-k6*Y(6)*Y(4);
Fv(7) = k6*Y(6)*Y(4);
Fv(8) = k8*Y(5)+k20*Y(18)*Y(4)-k21*Y(8)*Y(19);
Fv(9) = k9*Y(5);
Fv(10) = k9*Y(5);
Fv(11) = k10*Y(1)^2-k11*Y(11)-k12*Y(11)+k13*Y(12)*Y(13)-k16*Y(11)*Y(1)+k17*Y(16)+k19*Y(14)*Y(2)-k22*Y(19)*Y(11)+k23*Y(14);
Fv(12) = k12*Y(11)-k13*Y(12)*Y(13)-k14*Y(12)*Y(1);
Fv(13) = k12*Y(11)-k13*Y(12)*Y(13)-k15*Y(13)*Y(2)^2-k24*Y(19)*Y(13);%%%
Fv(14) = k14*Y(12)*Y(1)-k19*Y(14)*Y(2)+k22*Y(19)*Y(11)-k23*Y(14);
Fv(15) = k15*Y(13)*Y(2)^2-k18*Y(15)*Y(2)^2;
Fv(16) = k16*Y(11)*Y(1)-k17*Y(16);
Fv(17) = k18*Y(15)*Y(2)^2;
Fv(18) = k19*Y(14)*Y(2)-k20*Y(18)*Y(4)+k21*Y(8)*Y(19);
Fv(19) = k20*Y(18)*Y(4)-k21*Y(8)*Y(19)-k22*Y(19)*Y(11)+k23*Y(14)-k24*Y(19)*Y(13);%%%
end
  4 Comments
Torsten
Torsten on 8 Feb 2016
I'm surprised that MATLAB does not complain because Fv has to be a column vector.
Best wishes
Torsten.

Sign in to comment.

Accepted Answer

Phwey Gil
Phwey Gil on 9 Feb 2016
I solved the problem. I needed to call the Fv function with a function handle.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!