MATLAB Answers

0

How to frequently run ODE45?

Asked by Hamid Ahmadi Eshtehardi on 16 Nov 2019 at 17:30
Latest activity Edited by Walter Roberson
on 16 Nov 2019 at 21:08
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
function Kads = adsparameter(P,A,T,m)
Kb = 1.3806488e-23; %Boltzmann Constant
(J/K)
Kads = (P*A)/(sqrt(2*pi*m*Kb*T));
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)
k1ads = adsparameter(pa,1e-20,T,ma)*1e-13;
k1des = desparameter(1e-20,T,ma,2.8,1,80e3)*1e-
13;
k2ads = adsparameter(pb,1e-20,T,mb)*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;
k4ads = adsparameter(pc,1e-20,T,mc)*1e-13;
k4des = desparameter(1e-
20,T,mc,0.561,1,10e3)*1e-13;
r1f = k1ads*y(4);
r1b = k1des*y(3);
r2f = k2ads*(y(4)^2);
r2b = k2des*(y(2)^2);
r3f = kf*y(1)*y(2);
r3b = kb*y(3)*y(4);
r4f = k4ads*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 Comments

Walter Roberson
on 16 Nov 2019 at 20:07
I wonder if you should be using one of the "stiff" solvers such as ode23s ?
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.

Sign in to comment.

Products


Release

R2019b

1 Answer

Answer by Walter Roberson
on 16 Nov 2019 at 21:08
Edited by Walter Roberson
on 16 Nov 2019 at 21:08
 Accepted Answer

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

  0 Comments

Sign in to comment.