How to frequently run ODE45?

1 view (last 30 days)
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
Walter Roberson on 16 Nov 2019
I wonder if you should be using one of the "stiff" solvers such as ode23s ?
Hamid Ahmadi Eshtehardi
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.

Sign in to comment.

Accepted Answer

Walter Roberson
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

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!