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

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

Answer 1

[Jan]

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?

[Walter Roberson]

If you do the above then before each run

clear EVENTS

where EVENTS is the name of the function.

[Sofia Tsangaridou]

Thank you both for your insight.

@Jan would I need to know where the poles are to do that? Also, where and how do I call that EVENTS function? I copied and pasted your code and I ran into the same issue I had before (regardless of whether or not I add 'clear EVENTS'). My code just takes too long to run and it always gets stuck on line 415 in dde23. I'm probably doing something wrong though. For context, I called it by adding 'EVENTS(T,Y,Z)' above its definition.

[Torsten]

Include the t argument in the definition of the function ddefunc and tell us up to which value t proceeds.

[Sofia Tsangaridou]

@Torsten t = linspace(0,100). Is that what you mean?

[Torsten]

Replace

function dXdt = ddefunc(~,X,Z)

by

function dXdt = ddefunc(t,X,Z)

t

and tell us which values for t are reported.

[Sofia Tsangaridou]

@Torsten I have done that and it doesn't make a difference? What I mean is, the code runs too slowly so I'm not getting an output.

[Torsten]

Did you add the line

t

as the first line of the function ddefunc ?

Then each time ddefunc is called, the actual time t should be printed to your console.

Please tell us up to which time dde23 integrates and if it gets stuck at some time t > 0.

[Sofia Tsangaridou]

@Torsten I've just added that line and when I forcefully stop the code it still doesn't print any t values. It just shows this: "Operation terminated by user during dde23 (line 414)".

[Torsten]

Maybe you stopped too early :-)

Can you show us the complete code you are using at the moment?

[Sofia Tsangaridou]

@Torsten, hope this will make it make more sense. Thanks for being willing to help.

% Time-delay Values

taum = 8.09;

taue = 5;

% Implementation of the dde23 solver to solve system of DDEs

tspan = [0 100];

lags = [taum taue (taue + taum)]; % constant time-delays of the system

sol = dde23(@ddefunc,lags,@Xhist,tspan,@EVENTS);

%Plot

plot(sol.x,sol.y);

xlabel('time t');

ylabel('Y')

function [VALUE,ISTERMINAL,DIRECTION] = EVENTS(t,X,Z)

persistent count

if isempty(count)

count = 0;

end

count = count + 1;

VALUE = 0;

DIRECTION = 0;

ISTERMINAL = (count == 1e6);

end

% Definion of DDE system

function dXdt = ddefunc(t,X,Z)

t;

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

%Definition of the history functions

function h = Xhist(~)

Pstar = 24*(10)^9;

Tstar = 100;

P0 = Pstar;

Me0 = 1;

T0 = Tstar + 10;

h = [P0;

T0;

1;

1;

1;

Me0];

end

[Torsten]

Please just include the line

t

not

t;

[Sofia Tsangaridou]

@Torsten I let the code run for around half an hour and then I stopped it cause it was just taking too long. t got to 2.7826 and when I forcefully stopped it, it was still on line 415 of the dde23 code.

[Torsten]

Then I suggest you integrate up to t =2.78 and check whether the curves obtained so far are reasonable. If no, you have to check your equations. If yes, you have to be a little patient.

You can reduce output of the actual time by inserting

persistent count

if isempty(count)

count = 0;

end

count = count + 1;

if mod(count,1000) == 0

count = 0;

t

end

at the start of function ddefunc.

[Sofia Tsangaridou]

@Torsten thank you! I'll try that. Do you think that I should still keep the EVENTS function in as well?

[Torsten]

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]

Alright, I'll give that a try. Thanks for all your help. I really appreciate it :)