Modelling and Simulation: Regarding DDEs

Specific Domains Modelling & Simulations
,
1

I was facing a problem regarding delay. For low delay e.g. tau = 0.101 the code was fine as tau increases slightly as tau = 0.104 or more it shows some error not exactly error but warning and shows no promising result. While in Matlab those were fine there were no problem even for higher delay upto tau = 0.15! Why this is happening and how to overcome it!

Here is the .jl file:

function delay_model(du,u,h,p,t)

a,b,c,tau = p

hist2 = h(p, t-tau)[2]

du[1] = -a*u[1] + b*(hist2)*u[3]

du[2] = a*u[1] - b*(hist2)*u[3]  -c*u[2]^2 

du[3] = c*u[2]^2

end

h(p, t) = ones(3)

tau = 0.11

lags = [tau]

a = 0.04; b = 10; c = 3

p = (a,b,c,tau)

tspan = (0,20)

u0 = [1.0,0.9,1.0]

using DelayDiffEq

prob = DDEProblem(delay_model,u0,h,tspan,p; constant_lags=lags)

alg = MethodOfSteps(Tsit5())

sol = solve(prob,alg)

using Plots; plot(sol)

########

Here are the corresponding errors:
for tau = 0.104:
Warning: Interrupted. Larger maxiters is needed.

for tau = 0.108:
Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

Here are the attached file regarding various warnings:

delayerror
delayerror1241Ă—327 37.1 KB

2

Hello and welcome to the community! Your question needs a better title and tags for the right people to find it and help you. What packages do you need help with? Is it a question about modeling and simulation etc?

3

Yes! It is related to modeling and simulation

4

I can’t try it right now, but while you wait for an expert maybe try a stiff solver. Rosenbrock23 or Rodas5 for a start. What solver were you using in matlab?

5

In Matlab I used dde23.

6

This post is a great starting point for debugging issues like this: PSA: How to help yourself debug differential equation solving issues

The documentation suggests BS3 may be close to dde23: DDE Solvers · DifferentialEquations.jl

1 Like
7

What was the MATLAB code? How did you verify that the rhs was giving exactly the same derivatives?

8

Let me share corresponding MATLAB code:
%%Function
function dydt = Delay(t,y,Z);

global a b c

lagy = Z(:,1);

dydt=[-ay(1) + blagy(2)y(3);
ay(1) - blagy(2)y(3) - cy(2)^2;
cy(2)^2];

%%%%%

%%Output
clc;
clear all;
global a b c

a = 0.04; b = 10; c = 3;

delay = 0.1;
Exphst = [1.0 0.9 1.0];
sol = dde23(@Delay,[delay],Exphst,[0, 20]);
plot(sol.x,sol.y)

Now let see the corresponding results when tau = 0.1:

MATLAB result:
Julia_Delay_0.1

9

Julia Result:

Matlab_Delay_0.1
Matlab_Delay_0.1840Ă—630 6.7 KB

10

Now if I choose delay upto tau = 0.10519 Matblab shows the result as

Matlab_Delay_0.10519
Matlab_Delay_0.10519840Ă—630 15.6 KB

11

But for Julia it shows some warning like:
Warning: Interrupted. Larger maxiters is needed.

while it shows almost similar result for tau = 0.102918!! see the figure
Julia_Delay_0.102918

So the value of the delay is inconsistent for different coding language!! But it should not be happened right?

I would also like to share one thing that the system is unbounded regarding high delay. For example for the maximum time point 20, in MATLAB if tau crosses 0.10519 i.e., if tau = 0.105195 it shows the solution beyond the integration level. But in case of Julia same thing happens for lower delay. From which it is sure that the system is unbounded regarding high delay.

12

Did you try lower tolerances? Did you try increasing maxiters there?

13

I have tried with lower tolerance, found no improvement!! For “increasing of maxiters” I have no idea how to handle this.

14

I set a reminder for next week in case someone like @devmotion doesn’t find the time to take a look.

15

Thank you and will be waiting.

16

How are the history and initial conditions set in Matlab? Is Exphst the history? How do you set the initial conditions? If you set h(p, t) = [1.0,0.9,1.0] instead of h(p, t) = ones(3) (i.e., make initial conditions and history the same) what happens?

17

Exphst is the history and initial conditions (considered same in MATLAB).

Following @ANasc, initial conditions and history were set same in fact I omitted u0 and set h(p, t) = [1.0 0.9 1.0] and found the result was fine at tau = 0.105 (improved!). Then I set h(p, t) = [1.0 0.1 1.0], surprisingly the result turns even better at high value of tau (=0.14)!! But in this case I found another issue, if time is set tspan = (0,230) then it is fine, once it set as tspan=(0,240) an warning appear as

Warning: Instability detected. Aborting

and simulations stops accordingly!

18

what does that instability look like and why?

19

Setting h(p, t) = [1.0 0.1 1.0] and for tau = 0.14, it shows better result upto time scale tspan = (0,200)
See the fig specially time scale 0-200:
Julia_Delay_ts_200

but if time scale set up as tspan = (0,300) it shows that warning:
┌ Warning: Instability detected. Aborting

See the fig and see time scale was given (0,300) but if runs upto 240 approx:
Julia_Delay_ts_300

20

Can you open an issue on DelayDiffEq.jl?