dde23 solver gets stuck in an infinite loop (I think?)
4 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
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
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;
...
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!