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

4 ビュー (過去 30 日間)
Sofia Tsangaridou
Sofia Tsangaridou 2021 年 7 月 17 日
コメント済み: Sofia Tsangaridou 2021 年 7 月 18 日
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

回答 (1 件)

Jan
Jan 2021 年 7 月 17 日
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 件のコメント
Torsten
Torsten 2021 年 7 月 18 日
編集済み: Torsten 2021 年 7 月 18 日
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 2021 年 7 月 18 日
Alright, I'll give that a try. Thanks for all your help. I really appreciate it :)

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeLoops and Conditional Statements についてさらに検索

製品


リリース

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by