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

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

[Walter Roberson]

Hypothesize that one of the volunteers looking at your question is not familiar with chemical reaction problems. What would you advise the volunteer to pay attention to in order to determine whether the system is or is not working properly?

[Phwey Gil]

I added more details. I hope that it is more clear now.

[Torsten]

I'm surprised that MATLAB does not complain because Fv has to be a column vector.

Best wishes

Torsten.

[Star Strider]

It may pick that up from ‘initialvalues’ being a column vector.

Answer 1

[Phwey Gil]

I solved the problem. I needed to call the Fv function with a function handle.