Solving a system of ordinary differential equations

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));

[James Tursa]

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]

Hi @Emily,

Is it expected for the O3 concentration to increase without limit when O2 concentration is consistently supplied at 210,000?

Answer 1

[Sulaymon Eshkabilov]

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

Answer 2

[Torsten]

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

Answer 3

[William Rose]

@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

[William Rose]

@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.