ODE 45 Error in line 256 and line 343
2 views (last 30 days)
Show older comments
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.
3 Comments
Jan
on 1 May 2021
Why does your code stop, when you enter the debugging mode? What do you call "entering the debugging mode"?
Accepted Answer
Jan
on 1 May 2021
Edited: Jan
on 1 May 2021
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.
More Answers (0)
See Also
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!