"Busy" when running this code with ODE45
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Here is the code. I have also tried running it using ODE23 instead of 45, so I don't think that's the issue (I could be wrong though, I will take any help)
tspan = [0 300];
c0 = [0.01; 0.2; 0; 55.5; 0; 0;];
[t,C] = ode45(@ca2p1_v1,tspan,c0);
ext = (c0(1) - C(:,1))/(-1);
hold on
subplot(2,1,1)
plot(t,ext)
legend('Extent of Epoxidation Reaction');
xlabel('time');
ylabel('Extent');
subplot(2,1,2)
plot(t,C(:,1),t,C(:,2),t,C(:,3))
legend('Concentrations vs. Time')
xlabel('time');
ylabel('concentration');
legend('C8H16','H2O2','C8H16O')
function [dCdt] = ca2p1_v1(t,y)
%assign species
C8H16 = y(1);
H2O2 = y(2);
C8H16O = y(3);
H2O = y(4);
O2 = y(5);
CO2 = y(6);
%define rate constants
kE = 2.96e5;
kD = 1.5;
kC = 1.46e11;
%define rate laws
rE = kE * C8H16 * H2O2;
rD = kD * H2O2;
rC = kC * C8H16O * O2;
%define ODE for each species
dC8H16_dt = -rE;
dH2O2_dt = -rE - 2*rD;
dC8H16O_dt = rE - rC;
dH2O_dt = rE + 2*rD + 8*rC;
dO2_dt = rD + 11.5*rC;
dCO2_dt = 8*rC;
%derivative vector
dCdt = [dC8H16_dt; dH2O2_dt; dC8H16O_dt; dH2O_dt; dO2_dt; dCO2_dt];
end
1 Kommentar
Les Beckham
am 13 Feb. 2024
What is your question? FYI: it is normal for the Matlab status bar to show "Busy" when executing code.
Antworten (2)
Torsten
am 13 Feb. 2024
Use ode15s instead of ode45 - you have a stiff system of differential equations.
0 Kommentare
Walter Roberson
am 13 Feb. 2024
Verschoben: Walter Roberson
am 13 Feb. 2024
If you use ode15s and set the tspan to 15.01 or less than integration succeeds and shows most concentrations on the order of 10^-3 . But set the upper bound to 15.02 and you end earlier than that with concentrations going up to 10^8
If you use ode23s then it runs to completion.
tspan = [0 300];
c0 = [0.01; 0.2; 0; 55.5; 0; 0;];
[t,C] = ode23s(@ca2p1_v1,tspan,c0);
ext = (c0(1) - C(:,1))/(-1);
hold on
tiledlayout(4,2);
nexttile();
plot(t,ext)
title('Extent of Epoxidation Reaction');
xlabel('time');
ylabel('Extent');
nexttile();
plot(t,C(:,1));
title('C8H16')
xlabel('time');
ylabel('concentration');
nexttile();
plot(t,C(:,2));
title('H2O2')
xlabel('time');
ylabel('concentration');
nexttile();
plot(t,C(:,3));
title('C8H160')
xlabel('time');
ylabel('concentration');
nexttile();
plot(t,C(:,4));
title('H20')
xlabel('time');
ylabel('concentration');
nexttile();
plot(t,C(:,5));
title('O2')
xlabel('time');
ylabel('concentration');
nexttile();
plot(t,C(:,6));
title('CO2')
xlabel('time');
ylabel('concentration');
function [dCdt] = ca2p1_v1(t,y)
%assign species
C8H16 = y(1);
H2O2 = y(2);
C8H16O = y(3);
H2O = y(4);
O2 = y(5);
CO2 = y(6);
%define rate constants
kE = 2.96e5;
kD = 1.5;
kC = 1.46e11;
%define rate laws
rE = kE * C8H16 * H2O2;
rD = kD * H2O2;
rC = kC * C8H16O * O2;
%define ODE for each species
dC8H16_dt = -rE;
dH2O2_dt = -rE - 2*rD;
dC8H16O_dt = rE - rC;
dH2O_dt = rE + 2*rD + 8*rC;
dO2_dt = rD + 11.5*rC;
dCO2_dt = 8*rC;
%derivative vector
dCdt = [dC8H16_dt; dH2O2_dt; dC8H16O_dt; dH2O_dt; dO2_dt; dCO2_dt];
end
0 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!