ODE45has long runtime and graph will not plot

1 次查看（过去 30 天）
Yasmine
Yasmine 2024-1-31
回答： Alan Stevens 2024-2-2

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2024-2-2
Ran in:
Here's my attempt to make sense of your equations
tauspan = 0:100:10000;
% p0 = [0.5 0.5 0 0];
% dp1dt = -rf + rr
% dp2dt = 2*(-rf + rr)
% dp3dt = rf - rr
% dp4dt = 4*(rf - rr)
% d(2*p1-p2)/dt = 0 so 2*p1 - p2 = c2 (where c2 is a constant)
% d(p1 + p3)/dt = 0 so p1 + p3 = c3 (where c3 is a constant)
% d(4*p1+p4)/dt = 0 so 4*p1 + p4 = c4 (where c4 is a constant)
% Using initial conditions we have:
% 2*0.5 - 0.5 = c2 so p2 = 2*p1 - 0.5
% 0.5 + 0 = c3 so p3 = -p1 + 0.5
% 4*0.5 + 0 = c4 so p4 = -4*p1 + 2
kf = 1.32E10*exp(-236.7)/8.314;
% kr = 1.32E10*exp(-285.2)/8.314
% rf = kf*p1*(2*p1-0.5)^2/T
% rr = kr*(-p1+0.5)*(-4*p1+2)^2/T
% Scale the time base:
% tau = sf*t dp1dt = dp1dtau*dtau/dt = dp1dtau*sf
% sf*dp1dtau = -kf*p1*(2*p1-0.5)^2/T + kr*(-p1+0.5)*(-4*p1+2)^2/T
% Let kf/sf = 1 so
sf = kf;
% and kr/sf = exp(-285.2)/exp(-236.7) = exp(-48.5)
% dp1dtau = (-p1*(2*p1-0.5)^2 + exp(-48.5)*(-p1+0.5)*(-4*p1+2)^2)/T
p10 = 0.5;
T = 900:75:1200;
for i = 1:numel(T)
[tau, p1] = ode45(@(t,p) ODET(t,p,T(i)), tauspan, p10);
p2 = 2*p1-0.5;
p3 = -p1+0.5;
p4 = -4*p1+2;
t = tau/sf;
figure
hold on
Tlbl = ['T = ', int2str(T(i))];
subplot(2,2,1)
plot(t,p1),grid
title(Tlbl)
xlabel('t'), ylabel('p1')
subplot(2,2,2)
plot(t,p2),grid
title(Tlbl)
xlabel('t'), ylabel('p2')
subplot(2,2,3)
plot(t,p3),grid
title(Tlbl)
xlabel('t'), ylabel('p3')
subplot(2,2,4)
plot(t,p4),grid
title(Tlbl)
xlabel('t'), ylabel('p4')
hold off
end
function dpdtau = ODET(~,p,T)
dpdtau = (-p*(2*p-0.5)^2 + exp(-48.5)*(-p+0.5)*(-4*p+2)^4)/T;
end

类别

在 Help Center 和 File Exchange 中查找有关 Response Computation and Visualization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by