Solving a system of ordinary differential equations

14 visualizzazioni (ultimi 30 giorni)
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 Commenti
James Tursa
James Tursa il 6 Nov 2023
Modificato: James Tursa il 6 Nov 2023
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 il 7 Nov 2023
Hi @Emily,
Is it expected for the O3 concentration to increase without limit when O2 concentration is consistently supplied at 210,000?

Accedi per commentare.

Risposte (3)

Sulaymon Eshkabilov
Sulaymon Eshkabilov il 6 Nov 2023
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 il 6 Nov 2023
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 il 7 Nov 2023
Modificato: William Rose il 7 Nov 2023
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 Commento
William Rose
William Rose il 7 Nov 2023
Modificato: William Rose il 7 Nov 2023
[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.

Accedi per commentare.

Categorie

Scopri di più su Chemistry in Help Center e File Exchange

Prodotti


Release

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by