ODE45has long runtime and graph will not plot
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
1 commentaire
Réponses (1)
Alan Stevens
le 2 Fév 2024
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
0 commentaires
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!