# ODE45has long runtime and graph will not plot

2 views (last 30 days)
Yasmine on 31 Jan 2024
Answered: Alan Stevens on 2 Feb 2024
Walter Roberson on 31 Jan 2024
It is probably a stiff system; use ode15s instead of ode45

Alan Stevens on 2 Feb 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