Solving a system of ordinary differential equations

18 次查看（过去 30 天）
Emily
Emily 2023-11-6
编辑： William Rose 2023-11-7
I am trying to plot the O3 concentration over time, according to a system of differential equations that describe chemical rates and concentrations. I don't think the graph should be linear, but I keep getting a linear result. I am not experienced in MATLab and am looking for a solution. I tried using ode45 at first, but the function seemed to be too large to solve in a timely manner, so i switched over to ode15s. My code is below. Any help would be greatly appreciated. Please be descriptive in your answer, as I am new to MATLab code.
function f = ODEfun(t,Y)
NO=Y(1);
NO2=Y(2);
O=Y(3);
O3=Y(4);
C8H18=Y(5);
H2O=Y(6);
C8H17=Y(7);
O2C8H17=Y(8);
OC8H17=Y(9);
C7H17CHO=Y(10);
%k values
k1=0.4;
k2=21.83;
k3=26.6;
k4=12000;
k5=5272.43;
k6=1476.87;
k7=8551.09;
k8=3662.64;
%Constants
O2=210000;
OH=4.59E-08;
HO2=4.07E-06;
%Differential Equations
dNOdt = k1*NO2 - k3*O3*NO - k4*HO2*NO - k7*O2C8H17*NO;
dNO2dt = k3*O3*NO + k4*HO2*NO + k7*O2C8H17*NO - k1*NO2;
dOdt = k1*NO2 - k2*O2*O;
dO3dt = k2*O2*O - k3*O3*NO;
dC8H18dt = -k5*C8H18*OH;
dH2Odt = k5*C8H18*OH;
dC8H17dt = k5*C8H18*OH - k6*C8H17*O2;
dO2C8H17dt = k6*C8H17*O2 - k7*O2C8H17*NO;
dOC8H17dt = k7*O2C8H17*NO - k8*OC8H17*O2;
dC7H17CHOdt = k8*OC8H17*O2;
f = [dNOdt; dNO2dt; dOdt; dO3dt; dC8H18dt; dH2Odt; dC8H17dt; dO2C8H17dt; dOC8H17dt; dC7H17CHOdt];
end
tspan = [0 60]; %timespan
y0 = [0.0675 0.0075 0 0 0.15 0 0 0 0 0]; %initial conditions
[t y]=ode15s(@ODEfun,tspan,y0);
plot(t,y(:,4));
  2 个评论
James Tursa
James Tursa 2023-11-6
编辑：James Tursa 2023-11-6
There is no way we can help unless you show us the differential equations you are trying to solve so we can compare that to your actual code. E.g., post an image of the equations.
Sam Chak
Sam Chak 2023-11-7
Hi @Emily,
Is it expected for the O3 concentration to increase without limit when O2 concentration is consistently supplied at 210,000?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Sulaymon Eshkabilov
Ran in:
In fact, your solution code is ok but if you look at a bit closer by taking a shorter time span [0, 7] you will see the nonlinear behavior of the system at the start of the simulation. Besides, you may also try to set up ode settings with tighter error tolerances.
tspan = [0 7]; % Shorter timespan
y0 = [0.0675 0.0075 0 0 0.15 0 0 0 0 0]; %initial conditions
OPTs = odeset('RelTol',1e-9, 'AbsTol',1e-13);
[t, y]=ode15s(@ODEfun,tspan,y0, OPTs);
figure(2)
plot(t,y(:,4), 'r-', 'linewidth', 2);
grid on
function f = ODEfun(t,Y)
NO=Y(1);
NO2=Y(2);
O=Y(3);
O3=Y(4);
C8H18=Y(5);
H2O=Y(6);
C8H17=Y(7);
O2C8H17=Y(8);
OC8H17=Y(9);
C7H17CHO=Y(10);
%k values
k1=0.4;
k2=21.83;
k3=26.6;
k4=12000;
k5=5272.43;
k6=1476.87;
k7=8551.09;
k8=3662.64;
%Constants
O2=210000;
OH=4.59E-08;
HO2=4.07E-06;
%Differential Equations
dNOdt = k1*NO2 - k3*O3*NO - k4*HO2*NO - k7*O2C8H17*NO;
dNO2dt = k3*O3*NO + k4*HO2*NO + k7*O2C8H17*NO - k1*NO2;
dOdt = k1*NO2 - k2*O2*O;
dO3dt = k2*O2*O - k3*O3*NO;
dC8H18dt = -k5*C8H18*OH;
dH2Odt = k5*C8H18*OH;
dC8H17dt = k5*C8H18*OH - k6*C8H17*O2;
dO2C8H17dt = k6*C8H17*O2 - k7*O2C8H17*NO;
dOC8H17dt = k7*O2C8H17*NO - k8*OC8H17*O2;
dC7H17CHOdt = k8*OC8H17*O2;
f = [dNOdt; dNO2dt; dOdt; dO3dt; dC8H18dt; dH2Odt; dC8H17dt; dO2C8H17dt; dOC8H17dt; dC7H17CHOdt];
end
Torsten
Torsten 2023-11-6
Ran in:
It's just the starting phase where you observe a linear profile (see below).
tspan = [0 60000]; %timespan
y0 = [0.0675 0.0075 0 0 0.15 0 0 0 0 0]; %initial conditions
[t y]=ode15s(@ODEfun,tspan,y0);
plot(t,y(:,4));
function f = ODEfun(t,Y)
NO=Y(1);
NO2=Y(2);
O=Y(3);
O3=Y(4);
C8H18=Y(5);
H2O=Y(6);
C8H17=Y(7);
O2C8H17=Y(8);
OC8H17=Y(9);
C7H17CHO=Y(10);
%k values
k1=0.4;
k2=21.83;
k3=26.6;
k4=12000;
k5=5272.43;
k6=1476.87;
k7=8551.09;
k8=3662.64;
%Constants
O2=210000;
OH=4.59E-08;
HO2=4.07E-06;
%Differential Equations
dNOdt = k1*NO2 - k3*O3*NO - k4*HO2*NO - k7*O2C8H17*NO;
dNO2dt = k3*O3*NO + k4*HO2*NO + k7*O2C8H17*NO - k1*NO2;
dOdt = k1*NO2 - k2*O2*O;
dO3dt = k2*O2*O - k3*O3*NO;
dC8H18dt = -k5*C8H18*OH;
dH2Odt = k5*C8H18*OH;
dC8H17dt = k5*C8H18*OH - k6*C8H17*O2;
dO2C8H17dt = k6*C8H17*O2 - k7*O2C8H17*NO;
dOC8H17dt = k7*O2C8H17*NO - k8*OC8H17*O2;
dC7H17CHOdt = k8*OC8H17*O2;
f = [dNOdt; dNO2dt; dOdt; dO3dt; dC8H18dt; dH2Odt; dC8H17dt; dO2C8H17dt; dOC8H17dt; dC7H17CHOdt];
end
William Rose
William Rose 2023-11-7
编辑：William Rose 2023-11-7
Ran in:
@Emily,
It looks reasonable to me, especially when I simulate for a longer time.
tspan = [0 500]; %timespan
y0 = [0.0675 0.0075 0 0 0.15 0 0 0 0 0]; %initial conditions
[t y]=ode15s(@ODEfun,tspan,y0);
plot(t,y(:,1:4))
legend('NO','NO2','O','O3'); grid on
function f = ODEfun(t,Y)
NO=Y(1);
NO2=Y(2);
O=Y(3);
O3=Y(4);
C8H18=Y(5);
H2O=Y(6);
C8H17=Y(7);
O2C8H17=Y(8);
OC8H17=Y(9);
C7H17CHO=Y(10);
%k values
k1=0.4;
k2=21.83;
k3=26.6;
k4=12000;
k5=5272.43;
k6=1476.87;
k7=8551.09;
k8=3662.64;
%Constants
O2=210000;
OH=4.59E-08;
HO2=4.07E-06;
%Differential Equations
dNOdt = k1*NO2 - k3*O3*NO - k4*HO2*NO - k7*O2C8H17*NO;
dNO2dt = k3*O3*NO + k4*HO2*NO + k7*O2C8H17*NO - k1*NO2;
dOdt = k1*NO2 - k2*O2*O;
dO3dt = k2*O2*O - k3*O3*NO;
dC8H18dt = -k5*C8H18*OH;
dH2Odt = k5*C8H18*OH;
dC8H17dt = k5*C8H18*OH - k6*C8H17*O2;
dO2C8H17dt = k6*C8H17*O2 - k7*O2C8H17*NO;
dOC8H17dt = k7*O2C8H17*NO - k8*OC8H17*O2;
dC7H17CHOdt = k8*OC8H17*O2;
f = [dNOdt; dNO2dt; dOdt; dO3dt; dC8H18dt; dH2Odt; dC8H17dt; dO2C8H17dt; dOC8H17dt; dC7H17CHOdt];
end
  1 个评论
William Rose
William Rose 2023-11-7
编辑：William Rose 2023-11-7
@Emily,
[edit: adding a plot at the end, which I meant to include before]
Your system of equations leads to negative concentrations, which are of course impossible, at long time scales. Negative concentrations in the simulation suggests that (a) there is an error in the differential equations, or (b) an important limiting feature has not been accounted for. You have not provided the underlying chemical reactions, and I have not attempted to back-calculate them from your differential equations. Therefore I have not checked possibility (a). For possibility (b), I recommend that you consider adding a differential equation for [O2], instead of holding it constant. @Sam Chak asked about this above.
Now for the details. I have added code to plot the concentrations of all ten species.
The first plot below shows the results for tmax=10^5 (about 28 hours, assuming t is in seconds). It shows that [O3] (top panel) is increasing, and the other nine species have reached steady concentrations. (One possible concern: Around t=7e4-8e4, the concentrations of C8H17 and OC7H17 dip below zero by about 1e-20, then come back to tiny positive values. These negative concentrations are so tiny and relatively brief that I will not worry about them for now.)
In the next plot (below), we increase the simulation duration to 10^8 (about 3 years). Little has changed, except [O3] has continued to rise. The other species have not changed.
When we increase tmax from 10^8 (above) to 10^9 (below, about 32 years), we see a transient increase in [NO2] and [O], followed by a gradual decline in both. This is surprising, since everything seemed stable (if you view a slow, steady rise in [O3] as "stable") through t=10^8. See below.
When we increase tmax from 10^9 (above) to 10^10 (below, about 320 years), we see that [O3] is finally starting to decline, for the first time. Note that [NO2] becomes negative, and shortly thereafter, [O] becomes negative. When [O] becomes zero or negative, [O3] starts to fall. The differential equation for [O] suggests that [O] becomes negative because [NO2] went negative. And the differential equation for [O3] shows that it is the negative value of [O] that is causing the decline in [O3]. If you extend to tmax=10^11 (not shown), [O3] also enters negative territory.
Negative concentrations in the simulation suggest that (a) there is an error in the differential equations, or (b) an important limiting feature has not been accounted for. You have not provided the underlying chemical reactions, and I have not attempted to back-calculate them from your differential equations. Therefore I have not checked possibility (a). In possibility (b), I recommend that you consider adding a differential equation for [O2], instead of holding it constant.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Chemistry 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by