ODE 45 Error in line 256 and line 343

I have written a code as shown in below to describe a reactor, however, it gives an error on th ODE.

function dFdW=sabatier(catalyst,F)
 
%defining elements of output
 
Fh2o=F(1);
Fh2=F(2);
Fch4=F(3);
Fco=F(4);
Fco2=F(5);
T=F(6);
 
%Parameters
R= 8.314;%J/molK
 
P=55;
 
k1=8.90*10^8*exp(122.4/(R*T));
 
k2=3.42*10^6*exp(93.1/(R*T));
 
k3=9.22*10^-5*exp(104.8/(R*T));
 
K1eq=1.198*10^17*exp(-206/(R*T));
 
K2eq=1.767*10^-2*exp(41/(R*T));
 
K3eq=2.117*10^15*exp(-165/(R*T));
 
%Rate equations
Ft=Fh2o+Fco+Fh2+Fco2+Fch4;
Vt=Ft*R*T/P;
 
r1=k1/(Fh2*R*T/Vt)^2.5*(((Fch4*R*T/Vt)*(Fh2o*R*T/Vt))-((Fh2*R*T/Vt)^3*(Fco*R*T/Vt)/K1eq));
 
r2=k2/(Fh2*R*T/Vt)*(((Fco*R*T/Vt)*(Fh2o*R*T/Vt))-((Fh2*R*T/Vt)*(Fco2*R*T/Vt)/K2eq));
 
r3=k3/(Fh2*R*T/Vt)^3.5*(((Fch4*R*T/Vt)*(Fh2o*R*T/Vt)^2)-((Fh2*R*T/Vt)^4*(Fco2*R*T/Vt)/K3eq));
 
%Kinetic parameters (Cp)
%H2O
funone=@(Temp)32.243+(19.238*10^(-4)*Temp)+(10.555*10^(-6)*Temp.^2)+(-3.596*10^(-9)*Temp.^3);
%H2
funtwo=@(Temp)27.143+(92.738*10^(-4)*Temp)+(-1.381*10^(-5)*Temp.^2)+(76.451*10^(-10)*Temp.^3);
%CH4
funthree=@(Temp) 19.251+(52.126*10^(-3)*Temp)+(11.974*10^(-6)*Temp.^2)+(-1.132*10^(-8)*Temp.^3);
%CO
funfour=@(Temp)30.869+(-1.285*10^(-2)*Temp)+(27.892*10^(-6)*Temp.^2)+(-1.272*10^(-8)*Temp.^3);
%CO2
funfive=@(Temp)19.795+(73.436*10^(-3)*Temp)+(-5.602*10^(-5)*Temp.^2)+(17.153*10^(-9)*Temp.^3);

deltaH1o=-206.3; %kJ/mol
deltaH2o=41;%kJ/mol
deltaH3o=-165.1;%kJ/mol
 
deltaH1= -(integral(funfour,298,T)+integral(funtwo,298,T))+deltaH1o+ integral(funone,298,T)+integral(funthree,298,T);
 
deltaH2= -(integral(funfive,298,T)+integral(funtwo,298,T))+deltaH2o+ integral(funone,298,T)+integral(funfour,298,T);
 
deltaH3= -(integral(funfive,298,T)+integral(funtwo,298,T))+deltaH3o+ integral(funone,298,T)+integral(funthree,298,T);
 
%H2O
Cph2o=32.243+(19.238*10^(-4)*T)+(10.555*10^(-6)*T.^2)+(-3.596*10^(-9)*T.^3);
%H2
Cph2= 27.143+(92.738*10^(-4)*T)+(-1.381*10^(-5)*T.^2)+(76.451*10^(-10)*T.^3);
%CH4
Cpch4=19.251+(52.126*10^(-3)*T)+(11.974*10^(-6)*T.^2)+(-1.132*10^(-8)*T.^3);
%CO
Cpco=30.869+(-1.285*10^(-2)*T)+(27.892*10^(-6)*T.^2)+(-1.272*10^(-8)*T.^3);
%CO2
Cpco2=19.795+(73.436*10^(-3)*T)+(-5.602*10^(-5)*T.^2)+(17.153*10^(-9)*T.^3);

FCptot=Fh2o.*Cph2o+Fco.*Cpco+Fh2.*Cph2+Fco2.*Cpco2+Fch4.*Cpch4;
 
dFco2dW=r2+r3; 
dFh2dW=3*r1+r2+4*r3;
dFh2odW=-r1-r2-2*r3;
dFcodW=r1-r2;
dFch4dW=-r1-r3;
dTdW=(r1.*deltaH1+r2.*deltaH2+r3.*deltaH3)/(FCptot);
dFdW=[dFh2odW;dFh2dW;dFch4dW;dFcodW;dFco2dW;dTdW];

Also ı use below program to run the function

catalystspan=[0 1];
F0=[0.1 250 0.1 0.1 1000 500];
 
[catalyst,F]= ode45(@sabatier,catalystspan,F0);
Fh2o=F(:,1);
Fh2=F(:,2);
Fch4=F(:,3);
Fco=F(:,4);
Fco2=F(:,5);
T=F(:,6);
plot(catalyst,Fh2o,'g',catalyst,Fh2,'b',catalyst,Fch4,'r',catalyst,Fco,'y',catalyst,Fco2,'o',catalyst,T);

But it keeps giving error.

[Jan]

Please share the error message with us. You see this valuable information on the screen and we do need it for understanding the problem. "Error in line 256 and line 343" is impossible, because Matlab stops at the first error - there cannot be two errors.

[Bora Kader]

Thank you for your interest,

I do not get any error message in the command window. However, the program gets into debugging mode showing that it stuck on line256 for ODE45 and when ı press continue it stops again in line 343. It does not proceed further than that.

[Jan]

Why does your code stop, when you enter the debugging mode? What do you call "entering the debugging mode"?

Answer 1

[Jan]

For me your code works fine, but it is very slow. I've inserted some simplifications:

• 0.123*10^(-8) is an expensive power operation, while 0.123e-8 is a cheap constant.
• integral is a very expensive operation. In your code the same integrals are computed repeatedly:

deltaH1 = -(integral(funfour,298,T) + integral(funtwo,298,T)) + deltaH1o + ...
            integral(funone,298,T)+integral(funthree,298,T);
deltaH2 = -(integral(funfive,298,T) + integral(funtwo,298,T)) + deltaH2o + ...
            integral(funone,298,T)+integral(funfour,298,T);
deltaH3 = -(integral(funfive,298,T) + integral(funtwo,298,T)) + deltaH3o + ...
            integral(funone,298,T) + integral(funthree,298,T);

Avoid the repeated work:

% H2O
fcn1 = @(t) 32.243   + (19.238e-4*t)      + (10.555e-6*t.^2)   + (-3.596e-9*t.^3);
% H2
fcn2 = @(t) 27.143   + (92.738e-4*t)      + (-1.381e-5*t.^2)   + (76.451e-10*t.^3);
% CH4
fcn3 = @(t) 19.251   + (52.126e-3*t)      + (11.974e-6*t.^2)   + (-1.132e-8*t.^3);
% CO
fcn4 = @(t) 30.869   + (-1.285e-2*t)      + (27.892e-6*t.^2)   + (-1.272e-8*t.^3);
% CO2
fcn5 = @(t) 19.795   + (73.436e-3*t)      + (-5.602e-5*t.^2)   + (17.153e-9*t.^3);

Int1 = integral(fcn1, 298, T);
Int2 = integral(fcn2, 298, T);
Int3 = integral(fcn3, 298, T);
Int4 = integral(fcn4, 298, T);
Int5 = integral(fcn5, 298, T);

deltaH1 = -(Int4 + Int2) + deltaH1o + Int1 + Int3;
deltaH2 = -(Int5 + Int2) + deltaH2o + Int1 + Int4;
deltaH3 = -(Int5 + Int2) + deltaH3o + Int1 + Int3;

As next step the numerical integration of a polynomial can be accelerated massively, because it is trivial to perform the integration symbolically by hand:

Integral(a + bx + cx^2 + dx^3) = ax + bx^2/2 + cx^3/3 + dx^4/4:

% H2O
% fcn1 = @(t) 32.243   + (19.238e-4*t)      + (10.555e-6*t.^2)   + (-3.596e-9*t.^3);
fcn1_I = @(t) 32.243*t + (19.238e-4*t.^2/2) + (10.555e-6*t.^3/3) + (-3.596e-9*t.^4/4);
% H2
% fcn2 = @(t) 27.143   + (92.738e-4*t)      + (-1.381e-5*t.^2)   + (76.451e-10*t.^3);
fcn2_I = @(t) 27.143*t + (92.738e-4*t.^2/2) + (-1.381e-5*t.^3/3) + (76.451e-10*t.^4/4);
% CH4
% fcn3 = @(t) 19.251   + (52.126e-3*t)      + (11.974e-6*t.^2)   + (-1.132e-8*t.^3);
fcn3_I = @(t) 19.251*t + (52.126e-3*t.^2/2) + (11.974e-6*t.^3/3) + (-1.132e-8*t.^4/4);
% CO
% fcn4 = @(t) 30.869   + (-1.285e-2*t)      + (27.892e-6*t.^2)   + (-1.272e-8*t.^3);
fcn4_I = @(t) 30.869*t + (-1.285e-2*t.^2/2) + (27.892e-6*t.^3/3) + (-1.272e-8*t.^4/4);
% CO2
% fcn5 = @(t) 19.795   + (73.436e-3*t)      + (-5.602e-5*t.^2)   + (17.153e-9*t.^3);
fcn5_I = @(t) 19.795*t + (73.436e-3*t.^2/2) + (-5.602e-5*t.^3/3) + (17.153e-9*t.^4/4);

Int1 = fcn1_I(T) - fcn1_I(298);
Int2 = fcn2_I(T) - fcn2_I(298);
Int3 = fcn3_I(T) - fcn3_I(298);
Int4 = fcn4_I(T) - fcn4_I(298);
Int5 = fcn5_I(T) - fcn5_I(298);

deltaH1 = -(Int4 + Int2) + deltaH1o + Int1 + Int3;
deltaH2 = -(Int5 + Int2) + deltaH2o + Int1 + Int4;
deltaH3 = -(Int5 + Int2) + deltaH3o + Int1 + Int3;

Now a call of the function to be integrated takes 0.064 ms instead of 4.65 ms: 70 times faster.

Even then the intergration take a huge time. If I reduce the timestep to [0, 1e-7], ODE45 calculates 32'929 time steps. For the time span [0, 1] we can expect 1e7 * 33000 time steps, which take at least 0.064ms for calling the function to be integrated: 21120000 seconds or 244 days.

Well. Did you check if your problem is stiff? Then ODE45 is not a suiting integrator.

Using ODE15S instead needs 0.08 seconds for the complete time span [0 1]. Wow. A speedup of the factor 260e6 is the best I got in this forum and it will be hard to beat in the future!

The resulting diagram looks boring, because the values reach a stead state in the first microsecond. Maybe there is a bug in the formula?

By the way, do you see how the readability of the code is improved by using spaces? I'm too old to find typos in codes like this:

r1=k1/(Fh2*R*T/Vt)^2.5*(((Fch4*R*T/Vt)*(Fh2o*R*T/Vt))-((Fh2*R*T/Vt)^3*(Fco*R*T/Vt)/K1eq));
r2=k2/(Fh2*R*T/Vt)*(((Fco*R*T/Vt)*(Fh2o*R*T/Vt))-((Fh2*R*T/Vt)*(Fco2*R*T/Vt)/K2eq));
r3=k3/(Fh2*R*T/Vt)^3.5*(((Fch4*R*T/Vt)*(Fh2o*R*T/Vt)^2)-((Fh2*R*T/Vt)^4*(Fco2*R*T/Vt)/K3eq));

I've attached my complete code.

[Bora Kader]

Thank you soooo much. Your explanation had became useful in many ways. Additionally, thank you for the fast responce. I will analyze the code and work on the ways to correct that small error in the formula.