"Busy" when running this code with ODE45
显示 更早的评论
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 个评论
Les Beckham
2024-2-13
What is your question? FYI: it is normal for the Matlab status bar to show "Busy" when executing code.
回答（2 个）
Torsten
2024-2-13
0 个投票
Use ode15s instead of ode45 - you have a stiff system of differential equations.
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
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2024-2-13
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
