You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to plot bifurcation with Delay Differential equations?
56 views (last 30 days)
Show older comments
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 Comments
Kitipol Jankaew
on 16 Nov 2020
That is the critical point causes the stability of an equilibrium point to change. Here bifurcation point is . It depends on time delay τ.
Kitipol Jankaew
on 17 Nov 2020
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
on 17 Nov 2020
Edited: Kitipol Jankaew
on 17 Nov 2020
And I want get diagram look like this
kaushik dehingia
on 11 Feb 2021
Moved: Dyuman Joshi
on 15 Mar 2024
Can anyone share the Bifurcation diagram code for a delayed system? I t will be very helpful for me.
Accepted Answer
Alan Stevens
on 17 Nov 2020
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 Comments
Alan Stevens
on 17 Nov 2020
Edited: Alan Stevens
on 17 Nov 2020
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
on 18 Nov 2020
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
on 20 Nov 2020
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
on 20 Nov 2020
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
on 21 Nov 2020
Thank you very much. These are diagram which i want. Thank you very much again.
Kitipol Jankaew
on 23 Nov 2020
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
on 23 Nov 2020
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
on 9 Jan 2021
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
on 9 Jan 2021
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
on 10 Jan 2021
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
on 10 Jan 2021
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
on 2 Mar 2021
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
ibtissam benamara
on 2 Jun 2021
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
on 3 Jun 2021
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
on 20 Jun 2021
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
on 30 Jan 2022
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.
More Answers (1)
Priya Verma
on 15 Mar 2024
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 Comments
Priya Verma
on 15 Mar 2024
in second dde equation ...why are you defining denomenatoinator term in second delay varuabiable in matlab code ?
Torsten
on 15 Mar 2024
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
on 16 Mar 2024
May you provide MATLAB code to plot graph between two different lags (x-axis tau1 and y-axais tau2) in dde?
Torsten
on 16 Mar 2024
I don't understand what you mean by "graph between two different lags". Your differential equations only have one lag, namely tau.
Priya Verma
on 16 Mar 2024
In this model tau1 and tau2 two lags are given. So, how to plot graph between these two delays?
Torsten
on 16 Mar 2024
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]).
Torsten
on 17 Mar 2024
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
on 24 Mar 2024
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
on 24 Mar 2024
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
on 25 Mar 2024
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
on 26 Mar 2024 at 3:02
Is there any code, package, etc to fit the parameter values of dde?
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)