How to plot bifurcation with Delay Differential equations?

I want to draw the bifurcation diagram for the model.

All parameters are positve constant.

The value of parameters are as:

A1 = 0.8463, A2 = 0.6891, K = 1.2708, beta1 = 0.4110, beta2 = 0.1421,

The diagram are vary tau from 68 to 72 in steps of 0.001. For inital conditions X(0) = 0.26 and Y(0) = 0.58.

Please ansers me for Matlab code to plot the bifurcation diagrams.

[Alan Stevens]

What is the definition of the bifurcation point here?

[Kitipol Jankaew]

That is the critical point causes the stability of an equilibrium point to change. Here bifurcation point is . It depends on time delay τ.

[Alan Stevens]

How do you decide that has happened here?

[Kitipol Jankaew]

I had use dde23 to plot limit cycle. But now i want plot Hopf bifurcation diagram, I tried to use Runge-Kutta 4th order method to approximate the solution of system and plot the diagram, but i don't how to use this method with time delay.

Here Matlab code that I used.

clc;

clear all;

close all;

%constant

A1 = 0.8463;

A2 = 0.6891;

K = 4.2708;

b1 = 0.4110;

b2 = 0.1421;

%function

fx =@(t,x,y) x*(1-x/K)-A1*x*y/(1+x)-b1*x;

fy =@(t,x,y) A2*x*y/(1+y)-b2*y;

%initial

t(1)= 30;

x(1)= 0.9;

y(1)= 12;

h=0.01;

tfinal=1000;

N=ceil((tfinal-t(1))/h);

for j=1:N

t(j+1)=t(j)+h;

%tempx(j+1)=tempx(j)+h;

%tempy(j+1)=tempy(j)+h;

k1x=fx(t(j), x(j), y(j));

k1y=fy(t(j), x(j), y(j));

%k1y=fy(t(j), x(j), y(j), tempx(j), tempy(j));

k2x=fx(t(j)+h/2, x(j)+h/2*k1x, y(j)+h/2*k1y);

k2y=fy(t(j)+h/2, x(j)+h/2*k1x, y(j)+h/2*k1y);

%k2y=fy(t(j)+h/2, x(j)+h/2*k1x, y(j)+h/2*k1y, tempx(j)+h/2*k1x, tempy(j)+h/2*k1y);

k3x=fx(t(j)+h/2, x(j)+h/2*k2x, y(j)+h/2*k2y);

k3y=fy(t(j)+h/2, x(j)+h/2*k2x, y(j)+h/2*k2y);

%k3y=fy(t(j)+h/2, x(j)+h/2*k2x, y(j)+h/2*k2y, tempx(j)+h/2*k2x, tempy(j)+h/2*k2y);

k4x=fx(t(j)+h, x(j)+h*k3x, y(j)+h*k3y);

k4y=fy(t(j)+h, x(j)+h*k3x, y(j)+h*k3y);

%k4y=fy(t(j)+h, x(j)+h*k3x, y(j)+h*k3y, tempx(j)+h*k3x, tempy(j)+h*k3y);

x(j+1)=x(j)+h/6*(k1x+2*k2x+2*k3x+k4x);

y(j+1)=y(j)+h/6*(k1y+2*k2x+2*k3y+k4y);

end

%plot solution

figure;

plot(t,x,'.','markersize',3);

xlabel('t');

ylabel('x');

xlim([190 400])

figure;

plot(t,y,'.','markersize',3);

xlabel('t');

ylabel('y');

figure;

plot(x,y,'.','markersize',3);

xlabel('x');

ylabel('y');

[Kitipol Jankaew]

And I want get diagram look like this

[kaushik dehingia]

Can anyone share the Bifurcation diagram code for a delayed system? I t will be very helpful for me.

[kaushik dehingia]

Can anyone share me the bifurcation code?

Answer 1

[Alan Stevens]

How about the following for your loop (it assumed you have defined tau earlier in the file):

for j=1:N+1
    if t(j)<=tau
        xd = x(1); 
        yd = y(1);
    else
       d = ceil((t(j)-tau)/h);
       xd = x(d);
       yd = y(d);
    end
       
        
    t(j+1)=t(j)+h;
    %tempx(j+1)=tempx(j)+h;
    %tempy(j+1)=tempy(j)+h;
      
    k1x=fx(t(j),    x(j),   y(j));
    k1y=fy(t(j),    xd,   yd, y(j));
    %k1y=fy(t(j),    x(j),   y(j),   tempx(j),  tempy(j));
    k2x=fx(t(j)+h/2,    x(j)+h/2*k1x,   y(j)+h/2*k1y);
    k2y=fy(t(j)+h/2,    xd+h/2*k1x,   yd+h/2*k1y, y(j)+h/2*k1y);
    %k2y=fy(t(j)+h/2,    x(j)+h/2*k1x,   y(j)+h/2*k1y,  tempx(j)+h/2*k1x,    tempy(j)+h/2*k1y);
    
    k3x=fx(t(j)+h/2,    x(j)+h/2*k2x,   y(j)+h/2*k2y);
    k3y=fy(t(j)+h/2,    xd+h/2*k2x,   yd+h/2*k2y, y(j)+h/2*k2y); 
    %k3y=fy(t(j)+h/2,    x(j)+h/2*k2x,   y(j)+h/2*k2y,  tempx(j)+h/2*k2x,    tempy(j)+h/2*k2y);
    
    k4x=fx(t(j)+h,  x(j)+h*k3x, y(j)+h*k3y);
    k4y=fy(t(j)+h,  xd+h*k3x, yd+h*k3y, y(j)+h*k3y); 
    %k4y=fy(t(j)+h,  x(j)+h*k3x, y(j)+h*k3y,  tempx(j)+h*k3x,    tempy(j)+h*k3y);
    
    x(j+1)=x(j)+h/6*(k1x+2*k2x+2*k3x+k4x);
    y(j+1)=y(j)+h/6*(k1y+2*k2x+2*k3y+k4y);
 end

[Alan Stevens]

And fy woud have to change to

fy =@(t,xd,yd,y) A2*xd*yd/(1+yd)-b2*y;

and set

t(1) = 0;

Though I don't see why using an RK4 like this would give any significant advantage over dde23.

[Kitipol Jankaew]

Thank you very much. I have used dde23 to plot limit cycle but don't know how to use it to plot the bifurcation diagram like above diagram. Here Matlab code:

function exam5f % Code for plot graph of Stability and Bifurcation in Holling type2 predator-prey model

% with Alle effect and time delay

clear;

clc;

A1 = 0.7519;

A2 = 0.2287;

K = 6.4187;

b1 = 0.4110;

b2 = 0.1421;

function v = exam5d(t,y,Z)

ylag1 = Z(:,1);

%ylag2 = Z(:,1);

v = zeros(2,1);

v(1) = y(1)*(1-y(1)/K)-A1*y(1)*y(2)/(1+y(1))-b1*y(1);

v(2) = A2*ylag1(1)*ylag1(2)/(1+ylag1(1))-b2*y(2);

end

sol = dde23(@exam5d, 10.3340, [2, 1.2], [0,10000]);

figure;

%plot(sol.y(2,:),sol.y(1,:),'LineWidth',1);

plot(sol.y(1,:),sol.y(2,:),'LineWidth',1);

%title('Figure 1. Mealybugs and Green Lacewings.')

xlabel('X');

ylabel('Y');

figure;

%plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'LineWidth',1);

plot(sol.x,sol.y(1,:))%,'.','LineWidth',1);

title('Figure 1.Number of Mealybugs')

xlabel('t');

ylabel('X(t)');

xlim([1000 4000])

figure;

plot(sol.x,sol.y(2,:),'LineWidth',1);

%plot(sol.y(1,:),sol.y(2,:));

%title('Figure 1. Number of Green Lacewings.')

%xlabel('time t');

xlabel('t');

ylabel('Y(t)');

%ylabel('G');

%xlim([1000 7000])

end

And my advisor give the pascal program code which use Runge-Kutta method to plot the diagram. I tried to use that program but unsuccessful. So I change to use Matlab because Matlab has many examples code.

[Alan Stevens]

Can you upload the Pascal listing?

[Kitipol Jankaew]

I use Pascal XE. Here the code:

Program PPMWD;

Uses Crt,Dos;

Var

a1,a2,b1,b2,k,l: Double;

k1 : Double;

tau: Double;

T,X,Y,TempX,TempY,h :Double;

t2: Integer;

Xmax,Ymax : Double;

tpX,tpY : Double;

create,createX,createY,createK : Text;

Nstep,i,q: Longint;

DataX,DataY: file Of Double;

c : Integer;

file_name: String[25];

{------------------------------------------------------------------}

Function F1(T,X,Y:Double): Double;

Begin

F1 := X*(1-X/k1)-(a1*X*Y)/(1+X)-b1*X;

End;

Function F2(T,Y,TempX,TempY:Double): Double;

Begin

F2 := a2*TempX*TempY/(1+TempX)-b2*Y;

End;

{--------------------------------------------------------------------}

Procedure RungeKuttaSix;

Var

F11,F21,F31,F41,F51,F61: Double;

F12,F22,F32,F42,F52,F62: Double;

DummyT,DummyX,DummyY,DummydX,DummydY: Double;

Xmax,Ymax: Double;

tpX,tpY: Double;

create,createX,createY,createK : Text;

I : Longint;

c : Integer;

Begin{RungeKuttaSix}

F11 := h*F1(T,X,Y);

F12 := h*F2(T,Y,TempX,TempY);

DummyT := T+h/3;

DummyX := X+h/3;

DummyY := Y+h*F12/3;

DummydX := TempX+h*F11/3;

DummydY := TempY+h*F12/3;

F21 := h*F1(DummyT,DummyX,DummyY);

F22 := h*F2(DummyT,DummyY,DummydX,DummydY);

DummyT := T+h*2/3;

DummyX := X+h*F21*2/3;

DummyY := Y-F12/3+F22;

DummydX := TempX+h*F21*2/3;

DummydY := TempY-F12/3+F22;

F31 := h*F1(DummyT,DummyX,DummyY);

F32 := h*F2(DummyT,DummyY,DummydX,DummydY);

DummyT := T+h;

DummyX := X+h;

DummyY := Y+F12-F22+F32;

DummydX := TempX+h;

DummydY := TempY+F12-F22+F32;

F41 := h*F1(DummyT,DummyX,DummyY);

F42 := h*F2(DummyT,DummyY,DummydX,DummydY);

k := (F11+3*F21+3*F31+F41)/8;

l := (F21+3*F22+3*F32+F42)/8;

X := X+h*k;

Y := Y+h*l;

{TempX := TempX+h*k;

TempY := TempY+h*l; }

T := T+h;

End; {Procedure RungeKuttaSix}

{--------------------------------------------------------------------}

Begin {Main Program}

Clrscr;

a1 := 0.8463;

a2 := 0.6891;

k1 := 4.2708;

b1 := 0.4110;

b2 := 0.1421;

tau := 80.3340;

h := 0.001;

Nstep := 1000;

q := 0;

{initial value}

X := 1.6; Y := 1.2; T := 0;

Writeln('------WAIT FOR A HALF OF MINUTE------');

Writeln('Record ----> ');

Assign(create,'C:\Users\user\Desktop\test\test1\bidp.TXT');

Rewrite(create);

Assign(DataX,'C:\Users\user\Desktop\test\test1\XData2.dat'); Rewrite(DataX);

Write(DataX,X); Close(DataX); Writeln(create,T,'',X,'',Y,'');

Assign(DataY,'C:\Users\user\Desktop\test\test1\YData2.dat'); Rewrite(DataY);

Write(DataY,Y); Close(DataY); Writeln(create,T,'',X,'',Y,'');

Assign(createX,'C:\Users\user\Desktop\test\test1\tau-Xmax.TXT'); Rewrite(createX);

Assign(createK,'C:\Users\user\Desktop\test\test1\tau-Tmax.TXT'); Rewrite(createK);

Assign(createY,'C:\Users\user\Desktop\test\test1\tau-Ymax.TXT'); Rewrite(createY);

Writeln('RUNNING !!!');

Writeln('DO NOT TURN OFF COMPUTER.');

Writeln;

t2 := Round(tau/h);

While (tau<85) Do

Begin

c := 1;

{X := 0.2;}

tpX := X; Xmax := X;

Gotoxy(20,16);

Writeln(tau:5:20,' ',Xmax:5:20);

{Y := 0.4;}

tpY := Y; Ymax := Y;

Gotoxy(20,18);

Writeln(tau:5:20,' ',Ymax:5:20);

For i:=1 To Nstep Do

Begin {for}

RungeKuttaSix;

Gotoxy(20,8);

Writeln(T:5:20,' ',X:5:20,' ',Y:5:20);

Writeln(create,T:5:20,' ',X:5:20,' ',Y:5:20);

Gotoxy(20,4);

Write(q);

q := q+1;

If ((i>= 500) And (c=12))Then

Begin

If ((tpX<Xmax) And (Xmax>X)) Then

Gotoxy(20,24);

Writeln(T:5:20,' ',Xmax:5:20,' ');

Writeln(createX,tau,' ',Xmax);

{Writeln(createK,tau);}

If ((tpY<Ymax) And (Ymax>Y)) Then

Gotoxy(20,24);

Writeln(T:5:20,' ',Ymax:5:20,' ');

Writeln(createY,tau,' ',Ymax);

Writeln(createK,Xmax,' ',Ymax);

End;

If (c=12) Then c := 1;

If (c=4) Then

Begin

tpX := Xmax;

tpY := Ymax;

End;

If (c=8) Then

Begin

Xmax := X;

Ymax := Y;

End;

c := c+1;

End; {for}

tau := tau+0.01;

End; { while }

Close(createX);

Close(createY);

Close(createK);

Close(DataX);

Close(DataY);

Close(create);

Writeln('-----COMPLETED-----');

Writeln('-----PRESS SPACEBAR TO CONTINUE-----');

delay(100);

Repeat

Until Readkey=#32;

End.{Main Program}

[Alan Stevens]

Perhaps the following begins to do what you want:

tau = 0.01:0.05:6;
IC = [1.6 1.2];
for i = 1:numel(tau)
sol = dde23(@ddefn,tau(i),IC,[0 1000]);
x(i,:) = sol.y(1,end-50:end);
y(i,:) = sol.y(1,end-50:end);
end
subplot(2,1,1)
plot(tau,x,'k.')
subplot(2,1,2)
plot(tau,y,'r.')

function dxydt = ddefn(~,xy,Z)
        A1 = 0.8463; A2 = 0.6891; 
        K = 4.2708; 
        b1 = 0.4110; b2 = 0.1421;

        xd = Z(1);
        yd = Z(2);
        x = xy(1);
        y = xy(2);
        if x<10^-9
            x=0;
        end
        if y<10^-9
            y = 0;
        end
        dxydt = [x*(1-x/K)-A1*x*y/(1+x)-b1*x;
                 A2*xd*yd/(1+yd)-b2*y];
end

[Kitipol Jankaew]

Thank you very much. These are diagram which i want. Thank you very much again.

[Kitipol Jankaew]

I have equations. What are means of x(i,:) = sol.y(1,end-50:end);

y(i,:) = sol.y(1,end-50:end); ?

[Alan Stevens]

x(i,:) = sol.y(1,end-50:end);

This picks out the last 50 points to be plotted. The reason for selecting only the last few points is that you want to ensure you are plotting only points from the limit cycle, and not from the initial transient before the limit cycle is reached. Points from the transient part of the trace just muddy the picture.

The number 50 is arbitrary, you can change this.

[Tianyu Cheng]

I have a question, If there is no bifurcation, does the code works well? For exampe, I change some parameter and I know in this case there is no bifurcation, what is the figure look like?

And If I choose A1 as bifurcation parameter, How can I realise it?

Thank you very much

[Alan Stevens]

The code should still work ok if there is no bifurcation - only the graph will look much more boring! Try it and see.

[Tianyu Cheng]

Thank you very much. It works. And I have another question. This system is continous, thus the bifurcation figure should be smooth curve and not look like as the figure of discrete model. So how can I optimize this method?

Thank you very much.

[Alan Stevens]

There will always be gaps where the bifurcations occur! To get a smoother curve way from the bifurcations just use smaller intervals of tau.

[Houssem eddin Nezali]

con you help me for this on

dx/dt=hy-ax+yz

dy/dt=hx-by-xz

dz/dt=-dz+xy+w^2

dw/dt=xy+cw

x(0)=0.8 ,y(0)=0.3,z(0)=10.1,w(0)=4.5,a=4,b=-0.5,d=1,h=5,c=-2.2

bifuraction

[Alan Stevens]

doc ode45

[Houssem eddin Nezali]

ok

[Houssem eddin Nezali]

i need the the bifuraction fonction

[ibtissam benamara]

Hi

Thanks @Alan Stevens for your code is very helpful, but I have a question: why you use x(i,:) = sol.y(1, end-50:end); and y(i): ) = sol.y(1, end-50: (end);

I saw the same figure for prey and predator...? ?

[Alan Stevens]

To repeat my comment above:

x(i,:) = sol.y(1,end-50:end);

This picks out the last 50 points to be plotted. The reason for selecting only the last few points is that you want to ensure you are plotting only points from the limit cycle, and not from the initial transient before the limit cycle is reached. Points from the transient part of the trace just muddy the picture.

[ibtissam benamara]

this is not my question, i undersand why you use (1,end-50:end);

the question why you will take x(i,:)=y(i,:), consequetly, it will give the same figure for the two species, this not true;

[Akanksha Rajpal]

Hello @Alan Stevens

Your code really helped, but I was wondering if we can use similar coding if we want to extend the work to two delays? I tried that but it was showing error.

If you could help me regarding this and provide a code for this example only where the delay residing in X and Y are tau1 and tau2 respectively.

Answer 2

[Priya Verma]

In question, the denominator term is define in first delay variable term. Why are you all this term is defining in second delay term.

i. e. fy =@(t,x,y) A2*x*y/(1+y)-b2*y; in this denominator term is (1+y) .....?

A2*xd*yd/(1+yd)-b2*y; in this denominator term is (1+yd) .....?

please, explain...!

[Torsten]

Which code are you referring to where both of these lines appear ?

[Priya Verma]

in second dde equation ...why are you defining denomenatoinator term in second delay varuabiable in matlab code ?

[Priya Verma]

variable *

[Torsten]

Z(1) equals X(t-tau), Z(2) equals Y(t-tau) in the code.

Thus in MATLAB notation with xy(1) = X and xy(2) = Y:

dxydt(1) = xy(1)*(1-xy(1)/K)-A1*xy(1)*xy(2)/(1+xy(1))-b1*xy(1)
dxydt(2) = A2*Z(1)*Z(2)/(1+Z(1))-b2*xy(2)

[Priya Verma]

Yes, now correct 💯

[Priya Verma]

Thank you,

[Priya Verma]

May you provide MATLAB code to plot graph between two different lags (x-axis tau1 and y-axais tau2) in dde?

[Torsten]

I don't understand what you mean by "graph between two different lags". Your differential equations only have one lag, namely tau.

[Priya Verma]

I am taking about this model.

[Priya Verma]

In this model tau1 and tau2 two lags are given. So, how to plot graph between these two delays?

[Torsten]

dde23 gives solutions for X-,X+,Y,M and P as functions of t.

If you want to plot the solutions between tau1 and tau2 (assuming tau1 < tau2), restrict the plot to the interval [tau1 tau2] by setting xlim([tau1 tau2]).

[Priya Verma]

But, i don't understand, how to plot it.

[Torsten]

tau1 = 2;
tau2 = 4;
fun = @(t,y) y;
tspan = [0 5];
y0 = 1;
sol = ode45(fun,tspan,y0);
plot(sol.x,sol.y(1,:))
xlim([tau1 tau2])
grid on

[Priya Verma]

Have you taken tau2 on x-axis and tau1 on y-axis?

[Torsten]

tau1 and tau2 are just two numbers used in the delay differential equations (like tau1 = 1 and tau2 = 2). What do you want to plot there ?

[Priya Verma]

I want to plot domain of stability region with respect to tau1 and tau2 (i.e. on x-axis tau1 and on y-axis tau2) for the above model.

[Torsten]

I have no experience with stability regions for delay differential equations with respect to the delay vector. How do you determine this region numerically ?

[Priya Verma]

How to plot this type of graph for dde ?

[Torsten]

Looks like a plot of a solution variable at a certain time (I don't know which time) if the delay tau2 is varied from 0 to 200.

[Priya Verma]

Is there any code, package, etc to fit the parameter values of dde?

[Priya Verma]

or find tau value according to model?