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

Question

Phwey Gil 2016-2-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/267041-using-ode-to-solve-a-system-of-19-equations-there-seems-to-be-something-wrong-with-the-system

回答： Phwey Gil 2016-2-9

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