# How to frequently run ODE45?

2 views (last 30 days)
Hamid Ahmadi Eshtehardi on 16 Nov 2019
Edited: Walter Roberson on 16 Nov 2019
Hi everyone. I have written the following code in order to solve a set of ODE equations. As you can see in my function dy=equation(z,y) I have a constant named T (T=300 %temperatue).
function bimolecular_reaction
Kb = 1.3806488e-23; %Boltzmann Constant
(J/K)
end
function Kdes =
desparameter(A,T,m,tetarot,sigma,Edes)
Kb = 1.3806488e-23; %Boltzmann Constant
(J/K)
h = 6.2606957e-34; %Planck Constant (J/s)
R = 8.3144621; %universal gas constant
(J/mol.K)
Kdes =
((Kb*(T^3))/(h^3))*((A*2*pi*m*Kb)/(sigma*tetarot))*
exp((-Edes)/(R*T));
end
function Karr = surfaceparameter (T,Eac,nu)
%Kb = 1.3806488e-23; %Boltzmann Constant
(J/K)
%h = 6.2606957e-34; %Planck Constant (J/s)
R = 8.3144621; %universal gas constant
(J/mol.K)
Karr = nu*exp((-Eac)/(R*T));
end
function dy = equations(z,y)
T=300; %temperature
P = 20e5; %total presure of system (Pa)
pa = (2/3)*P; %partial pressure of CO (Pa)
pb = (1/3)*P; %partial pressure of O2 (Pa)
pc = 0; %partial pressure CO2 (pa)
ma = 28 * 1.66054e-27; % mass of CO (Kg)
mb = 32 * 1.66054e-27; %mass of O2 (Kg)
mc = 80 * 1.66054e-27; %mass of CO2 (Kg)
k1des = desparameter(1e-20,T,ma,2.8,1,80e3)*1e-
13;
k2des = desparameter(1e-20,T,mb,2.08,2,40e3)*1e-
13;
kf = surfaceparameter(T,120e3,1e13)*1e-13;
kb = surfaceparameter(T,180e3,1e13)*1e-13;
k4des = desparameter(1e-
20,T,mc,0.561,1,10e3)*1e-13;
r1b = k1des*y(3);
r2b = k2des*(y(2)^2);
r3f = kf*y(1)*y(2);
r3b = kb*y(3)*y(4);
r4b = k4des*y(3);
dy = zeros (4,1);
dy(1) = r1f-r1b-r3f+r3b;
dy(2) = (2*r2f)-(2*r2b)-r3f+r3b;
dy(3) = r4f-r4b+r3f-r3b;
dy(4) = r1b-r1f-(2*r2f)+(2*r2b)+r3f-r3b-r4f+r4b;
end
[z,y] = ode45(@equations,[0 1],[0 0 0 1]);
rco2 = desparameter(1e-20,300,80 * 1.66054e-
27,0.561,1,10e3)*y(end,3);
plot(z,y(:,2),'r')
end
Firstly, I am going to construct a loop by which I could solve my ode system based on different values of T (between 300 to 1000) but I do not know how can do so. Secondly, since, the constants of my ode system is of order of 10^13 it will takes me too much time to reach to an answer for my ode system how can I also solve this problem?
##### 2 CommentsShowHide 1 older comment
Hamid Ahmadi Eshtehardi on 16 Nov 2019
Thank's for your comment. I have already read the parameterzing functions but I did not understand that how can I use them in my code. Could you please give me a little bit more advise in this? About the solver you were right my system was stiff so I had to use one of stiff solvers such as ode15s or ode23s.

Walter Roberson on 16 Nov 2019
Edited: Walter Roberson on 16 Nov 2019
for T = 300:100:1000
[z,y] = ode23s(@(t,y) equations(t,y,T), [0 1], [0 0 0 1]);
rco2 = desparameter(1e-20,T,80 * 1.66054e-27,0.561,1,10e3)*y(end,3);
plot(z,y(:,2),'r', 'DisplayName', sprintf('T=%g', T));
legend show
hold on
drawnow();
end
function dy = equations(z,y,T)
%T=300; %temperature
etc

R2019b

### Community Treasure Hunt

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

Start Hunting!

Translated by