MATLAB Answers

dde23 solver gets stuck in an infinite loop (I think?)

1 view (last 30 days)
Sofia Tsangaridou
Sofia Tsangaridou on 17 Jul 2021
Commented: Sofia Tsangaridou on 18 Jul 2021
So, I'm trying to solve a system of 6 delay differential equations using the dde23 solver but I think my system is too complicated? Whenever I try to run it, it just doesn't do anything and when I forcefully stop it (using Ctrl+C) I'm made aware that it gets stuck either on line 414 or 415 of the dde23 code. Is there a way around this?
This is my ddefunc (saved as a separate script):
function dXdt = ddefunc(~,X,Z)
D0 = 0.21829;
betaP = 8.6;
Vm = (4*pi*(21)^3)/24;
kP = 0.0072419;
Qstar = 1.1;
gamaP = 0.05;
alphaP = 212.95;
bP = 308;
nP = 2;
Tprod = 61.6;
gamaT = 0.01;
alphaT = 0.00075*144.87;
kS = 2/3;
kT = 0.003*3180;
nT = 2;
etammin = 0.38874;
etammax = 2.6828;
bm = 706;
etaemin = 0.41022;
etaemax = 0.69335;
be = 92.1;
lag1 = Z(:,1); % lag1 = taum
lag2 = Z(:,2); % lag2 = taue
lag3 = Z(:,3); % lag3 = taue + taum
dXdt = [D0/betaP*Vm*kP*Qstar*X(4)*X(5) - gamaP*X(1) - alphaP*X(1)^nP/(bP^nP + X(1)^nP); % X(1) = P
Tprod - gamaT*X(2) - alphaT*(X(6) + kS*betaP*X(1))*X(2)^nT/(kT^nT + X(2)^nT); % X(2) = T
X(3)*((etammax - etammin)*(X(2)/(bm +X(2)) - lag1(2)/(bm + lag1(2)))); % X(3) = phi1
X(4)*((etammax - etammin)*(lag2(2)/(bm + lag2(2)) - lag3(2)/(bm + lag3(2)))); % X(4) = phi2
X(5)*((etaemax - etaemin)*(X(2)/(be + X(2)) - lag2(2)/(be + lag2(2)))); % X(5) = psi
Vm*kP*Qstar*(X(3) - X(4)*X(5)) + X(6)*(etaemin + (etaemax - etaemin)*X(2)/(be + X(2)))]; % X(6) = Me
end

Answers (1)

Jan
Jan on 17 Jul 2021
dde23 cannot stuck in an infinite loop. But it is easy to provide a system, which has a pole, such that the step size gets tiny and the computing time can be very slow.
You can use an event function to stop the integration, when it was called a certain number of times:
function [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y,Z)
persistent count
if isempty(count)
count = 0;
end
count = count + 1;
VALUE = 0;
DIRECTION = 0;
ISTERMINAL = (count == 1e6);
end
Now display the trajectory you get. Does the function run in a pole?
  16 Comments
Sofia Tsangaridou
Sofia Tsangaridou on 18 Jul 2021
Alright, I'll give that a try. Thanks for all your help. I really appreciate it :)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!