dde23 solver gets stuck in an infinite loop (I think?)
7 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Jan
am 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 Kommentare
Torsten
am 18 Jul. 2021
Bearbeitet: Torsten
am 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;
...
Siehe auch
Kategorien
Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!