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

4 views (last 30 days)
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
Torsten
Torsten on 18 Jul 2021
Edited: Torsten on 18 Jul 2021
You could combine both (continuous information about the actual integration time and stopping after a certain number of time steps) by leaving the function ddefunc as is (without any changes concerning t) and changing the events function by inserting three new lines
...
count = count + 1;
if mod(count,1000) == 0
T
end
VALUE = 0;
...
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.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!