What is a simple code for the chapman cycle involving only hydrogen and oxygen.This code should produce a curve but it only produces straight lines?

function Chapman_cycle
y0= [2e17 1e-6 5e12];
tspan=[0 100];
%tspan=[0,1.0e-5,1.0e-3,0.1,0.5,1,2,4,6,8,10,20,30,40,50,60];
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,y]=ode15s(@oz,tspan,y0,options);
disp('   iteration       t           y1           y2          y3  ');
for i=1:length(t)
fprintf('%12.1e %12.1e %12.2e %12.2e \n',i,t(i),y(i,1),y(i,2),y(i,3));
end
%Plot(t,y(:,1),t,y(:,2));
semilogy(t,y(:,1),t,y(:,2),t,y(:,3));
title(['Solution of the Chapman Cycle']);
xlabel('t');
ylabel('solution y');
%axis([tspan(1) tspan(end) 1.0e-12 10]);
%-------------------------------------------------------------
function f = oz(t,y);
f=zeros(3,1);
k1= 1.03e-94;     %y(1) = O2, y(2) = O = sum , y(3) = O3,
k2= 5.92e-34;     
k3= 4.38e-26;

f(1)= -k1*y(1)-k2*y(1)*y(2)+k3*y(3);
f(2)= 2*k1*y(1)-k2*y(1)*y(2)+k3*y(3);
f(3)= k2*y(2)*y(1)-k3*y(3);

[Walter Roberson]

If it produces straight lines in logy then that implies that it produces curves in y.

[Matthew Worker]

so how do I alter my code to show this? thank you.

[Walter Roberson]

Change

semilogy

to

plot

[Matthew Worker]

I did this and no lines appear on my figure now??

[David Goodmanson]

I think you need to take a look at the units of concentration and rate constants this expression. [note added: By 'take a look at' I mean 'find out what is wrong with' because the rate constants are inconsistent with the time scale of 0-100]. At t=0, dy/dt for all three components of y is due to the k(3)*y(3) term which is the largest. That term is equal to 2.19e-13. Your time interval is from 0 to 100 so of course you see no change in any of the y components.

[John D'Errico]

Yes. It could easily be that there are some numerically incorrect constants here. If not, then much higher precision would be needed to solve this, as well as a much longer simulation time.

Answer 1

[John D'Errico]

You have massive problems with the dynamic range of your numbers here, to the extent that you literally cannot solve this in double precision arithmetic.

Look here:

k1= 1.03e-94;     %y(1) = O2, y(2) = O = sum , y(3) = O3,
k2= 5.92e-34;     
k3= 4.38e-26;

y0= [2e17 1e-6 5e12];

See that the k coefficients vary by 64 powers of 10.

A simple rule of double precision arithmetic is that if we try to add (or subtract) numbers together that span more than 16 powers of 10, then we will have a numerical problem.

Next, look at how we will use those numbers. For example:

f(1)= -k1*y(1)-k2*y(1)*y(2)+k3*y(3);

k1 is SO immeasurably tiny that the contribution from the first term will never be significant, until y(1) grows by dozens of powers of 10.

And as David points out, at t == 0, these terms are dominated numerically by the k3*y(3) term, which is itself tiny.

k3*y(3)
ans =
     2.19e-13

f
f =
     2.19e-13     2.19e-13    -2.19e-13

f is a vector of differential components. They are all so small that over a time span from 0 to 100, you will only see almost no change. At best, it will be linear, and the only y component that will show even the tiniest amount of difference will be y(2).

y(1) and y(3) are both numbers that are more than 16 powers of 10 larger than this. So over a time span of only 0-100, y(1) and y(3) will be EXACTLY constant values in double precision. Even if you used higher precision arithmetic, that would be irrelevant here. 100 time units will simply not be a long enough time to develop any significant chance in y(1) and y(3).

And while y(2) will see SOME change in that period of time, the difference in y(2) will be TINY. Thus on a period of time that goes from 0-100, we can expect y(2) to change by only approximately 2.19e-13*100 = 2.19e-11. And since y(2) is initially 1e-6, the change in y(2) will be seen down in only the 5 significant digit in that period of time.

As I said, you CANNOT solve this problem in MATLAB using double precision arithmetic. The variation in dynamic range of the numbers involved will cause your arithmetic to fail. Whoever posed this problem for you to solve would surely have known and expected that.