How to plot bifurcation with Delay Differential equations?

2020 11 14
2 个回答
回答已采纳
更新时间：2024 3 26
15 次查看（30 天）

0 个投票

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更早的评论

Alan Stevens
Alan Stevens 2020-11-14
What is the definition of the bifurcation point here?
Kitipol Jankaew
Kitipol Jankaew 2020-11-16
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
Alan Stevens 2020-11-16
How do you decide that has happened here?
Kitipol Jankaew
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');
Kitipol Jankaew
Kitipol Jankaew 2020-11-17
编辑：Kitipol Jankaew 2020-11-17
And I want get diagram look like this
kaushik dehingia
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
kaushik dehingia 2021-2-11
Can anyone share me the bifurcation code?

请先登录，再进行评论。

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

 采纳的回答

Alan Stevens
Alan Stevens 2020-11-17

0 个投票

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更早的评论

Alan Stevens
Alan Stevens 2020-11-17
编辑：Alan Stevens 2020-11-17
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
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.
Alan Stevens
Alan Stevens 2020-11-18
Can you upload the Pascal listing?
Kitipol Jankaew
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}
Alan Stevens
Alan Stevens 2020-11-20
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
Kitipol Jankaew 2020-11-21
Thank you very much. These are diagram which i want. Thank you very much again.
Kitipol Jankaew
Kitipol Jankaew 2020-11-23
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
Alan Stevens 2020-11-23
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
Tianyu Cheng 2021-1-9
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
Alan Stevens 2021-1-9
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
Tianyu Cheng 2021-1-10
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
Alan Stevens 2021-1-10
There will always be gaps where the bifurcations occur! To get a smoother curve way from the bifurcations just use smaller intervals of tau.
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
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
Alan Stevens 2021-6-3
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
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
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.

请先登录，再进行评论。

更多回答（1 个）

Priya Verma
Priya Verma 2024-3-15

0 个投票

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...!

21 个评论
显示 19更早的评论 隐藏 19更早的评论

Torsten
Torsten 2024-3-15
Which code are you referring to where both of these lines appear ?
Priya Verma
Priya Verma 2024-3-15
in second dde equation ...why are you defining denomenatoinator term in second delay varuabiable in matlab code ?
Torsten
Torsten 2024-3-15
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
Priya Verma 2024-3-16
May you provide MATLAB code to plot graph between two different lags (x-axis tau1 and y-axais tau2) in dde?
Torsten
Torsten 2024-3-16
I don't understand what you mean by "graph between two different lags". Your differential equations only have one lag, namely tau.
Priya Verma
Priya Verma 2024-3-16
In this model tau1 and tau2 two lags are given. So, how to plot graph between these two delays?
Torsten
Torsten 2024-3-16
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
Priya Verma 2024-3-17
But, i don't understand, how to plot it.
Torsten
Torsten 2024-3-17
编辑：Torsten 2024-3-17
Ran in:
Ran in:
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
Priya Verma 2024-3-17
Have you taken tau2 on x-axis and tau1 on y-axis?
Torsten
Torsten 2024-3-17
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
Priya Verma 2024-3-24
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
Torsten 2024-3-24
I have no experience with stability regions for delay differential equations with respect to the delay vector. How do you determine this region numerically ?
Torsten
Torsten 2024-3-25
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
Priya Verma 2024-3-26
Is there any code, package, etc to fit the parameter values of dde?
Priya Verma
Priya Verma 2024-3-26
or find tau value according to model?

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

产品

版本

R2019a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by