How to function 𝑎𝐴 + 𝑏𝐵 → 𝑝P in ODE89
Mostrar comentarios más antiguos
𝑎𝐴 + 𝑏𝐵 → 𝑝P
𝑑𝐴/𝑑𝑡 = −𝐾 ∗ 𝐴 ∗ 𝐵 𝑑𝐵/𝑑𝑡 = (𝑏/𝑎) ∗ (𝑑𝐴/𝑑𝑡) = −𝑌𝐵 ∗ (𝐾 ∗ 𝐴 ∗ 𝐵) 𝑑𝑃/𝑑𝑡 = −(𝑝/𝑎) ∗ (𝑑𝐴/𝑑𝑡) = 𝑌𝑃 ∗ (𝐾 ∗ 𝐴 ∗ 𝐵)
Respuestas (1)
Davide Masiello
el 7 de Mzo. de 2022
Editada: Davide Masiello
el 7 de Mzo. de 2022
This should work:
clear,clc
tspan = [0,10];
y0 = [1,1,0];
[t,y] = ode89(@yourODEsystem,tspan,y0);
plot(t,y)
legend('A','B','P','Location','best')
function out = yourODEsystem(t,y)
% Coefficients
K = 1;
a = 2;
b = 1;
p = 0.5;
% Variables
A = y(1);
B = y(2);
P = y(3);
% Time derivatives
dAdt = -K*A*B;
dBdt = -(b/a)*K*A*B;
dPdt = (p/a)*K*A*B;
% Output
out = [dAdt;dBdt;dPdt];
end
Just replace you actual values of stoichiometric coefficients and kinetic constants.
6 comentarios
Torsten
el 7 de Mzo. de 2022
Have in mind that chemical kinetics ODEs are usually stiff.
Use ode15s instead of ode89.
Navaneetha Krishnan Murugadoss
el 7 de Mzo. de 2022
Davide Masiello
el 7 de Mzo. de 2022
Editada: Davide Masiello
el 7 de Mzo. de 2022
Please show all your code. Moreover, please indicate what version of Matlab you are using.
Navaneetha Krishnan Murugadoss
el 7 de Mzo. de 2022
Navaneetha Krishnan Murugadoss
el 7 de Mzo. de 2022
Davide Masiello
el 7 de Mzo. de 2022
The function call in ode89 must be equal to the function name. Write this
clear,clc
tspan = [0,12];
y0=[0 1 3];
[t,y] = ode89(@DEdef,tspan,y0);
plot(t,y)
legend('CL','NOM','DBP','Location','best')
function Ddv_div = DEdef(t,y)
% Coefficients
K = 5E-5;
YB=1;
YP=0.15;
% Variables
A = y(1);
B = y(2);
P = y(3);
% Output
Ddv_div = [-K*A*B;-YB*(K*A*B);YP*(K*A*B)];
end
However, let me point out that if the initial concentration of one of the two reactants is zero (like in your case) you won't observe any change in the concentration of any of the compounds, since the reaction cannot occur.
Categorías
Más información sobre Chemistry en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!