How to plot bifurcation with Delay Differential equations?

18 次查看（过去 30 天）

显示更早的评论

Kitipol Jankaew 2020-11-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/647578-how-to-plot-bifurcation-with-delay-differential-equations

评论： Priya Verma 2024-3-26

采纳的回答： Alan Stevens

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Kitipol Jankaew 2020-11-17

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');

kaushik dehingia 2021-2-11

移动：Dyuman Joshi 2024-3-15

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

kaushik dehingia 2021-2-11

Can anyone share me the bifurcation code?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alan Stevens 2020-11-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/647578-how-to-plot-bifurcation-with-delay-differential-equations#answer_546928

在 MATLAB Online 中打开

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

20 个评论
显示 18更早的评论隐藏 18更早的评论

Kitipol Jankaew 2020-11-18

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.

Kitipol Jankaew 2020-11-20

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}

ibtissam benamara 2021-6-20

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 2022-1-30

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.

请先登录，再进行评论。