How do I plot a bifurcation diagram with known dynamical equations of motion (> 3 dimension)?

I have equation of motions for a system.The system is given by,

d[y]/dt=[f(x,t,parameters)] where, [y] is say 5 dimensional. I want to see how the solution changes with parameter. How do I write a code for that? The system is non-linear. Earlier I tried to plot bifurcation tree diagram for the system by storing variable after each period,by changing the parameter gradually to observe period doubling bifurcation, but in some cases I find some discontinuities that I don't quite understand. Any reference on how to correctly plot bifurcation diagram would be appreciated.

Answer 1

[William Rose]

@Rebeka,

Let's suppose that y is 1x5 and that the parameter is A, as in your diagram. Let's also assume that the equaitons of the systyem are too complicated to solve analytically. Therefore we will have to rely on numerical simulation. Make vector A with size Nx1, containing all the A values you want to have on your bifurcaiton diagram. Simulate the system with A=A(1). After a suitable simulation duraiton, evalute y(:) to see if it has reached a steady state (1). If y has reached a steady state, save the steady state values of y(1:5) as one row of a Nx5 array, for example yss1() (for steady state solution 1). We will use array yss1 for plotting later. Repeat with all N values of A. If you find no steady state solution for a particular value of A, save a row of NaNs in yss1(), for that particular vlaue of A. After you have done all N values of A, plot yss1(:,1) versus A. If there are multiple steady state solutions , you will have to find them by repeating the process above, but with different initial conditions. Your knowledge of the system should help you in this regard. On this second pass, in which you start from different inital points than what you used before, save the steady state solutions in yss2(). When done, add the plot of yss2(:,1) to the plot.

1. To see if the system has reached a steady state, check y(:,1), y(:,2), ..., y(:,5) for multiple time points at the end of the simulation, to see if the vector y has stabilized.

I am attaching a recent publication (here, but may be paywalled), in which we published experimental and theoretical bifurcation diagrams for a complex system. Figures 4-9 are bifurcation diagrams in which the parameter on the horizontal axis is the frequency of external stimulation. This was a study of the repsonse of a non-linear optical system to external laser modulation. We did not do the bifurcaiton analysis as described above, because "steady state" meant a phase-locked solution, which was not a steady state in the traditional sense.

[William Rose]

@Rebeka,

My answer above is just a suggesition. I apologize if it sounded like "you have to do it like this". You know your system better than anyone. If the bifurcaiton diagram in your posting was generated by you, feel free to share the equations or code that made it. I bet some of the smart and helpful people on this site will be able to assist you.

[Rebeka]

I agree with you. What if the stable state is a periodic orbit rather than a fixed point? As in my case I am expecting to observe how the periodic solution changes with the change of a parameter ('A'). Also, in my code, with some knowledge of the dynamics I am storing the last dataset after every driving period for 50 iteration at a fixed parameter (A). Then I am gradually changing the parameter ('A'). I am attching here a plot, and the code I have used to derive it with. It would be of great help if someone can look into it to see if I am going about the problem correctly. If the result is correct then what are the types of bifurcation there?

p.s Is there any other way to see the bifurcation of solutions

clc
clear all
mu=1; %mass ratio(m_2/m_1)
a1=1; %length ratio(l_1/l_2)
b_1=0.1;
b_2=b_1/mu;
tic;
w=1.779;
g=1;
A=0.5;
p=15000;
t_f=p;
step=0.01;
tt=t_f/step;
k_s=1; %natural frequency of internal pendulum (sqrt(g/l))
f1=@(a,b,c,d,e) mu*sin(b-a)*(a1*cos(b-a)*c^2+d^2)/(a1*(1+mu*(sin(b-a))^2))+(mu*sin(2*b-a)-(2+mu)*sin(a))*k_s^2/(2*(1+mu*(sin(b-a))^2))+A*((2+mu)*sin(a)-mu*sin(2*b-a))*sin(e)*k_s^2/(2*(1+mu*(sin(b-a))^2))+b_1*A*(sin(2*b-a)-3*sin(a))*cos(e)*k_s^2/(w*(1+mu*(sin(b-a))^2))-2*b_1*(2-(cos(b-a))^2)*c/(1+mu*(sin(b-a))^2);
f2=@(a,b,c,d,e) -sin(b-a)*((1+mu)*a1*c^2+mu*cos(b-a)*d^2)/(1+mu*(sin(b-a))^2)-a1*(1+mu)*(sin(b-2*a)+sin(b))*(k_s^2)/(2*(1+mu*(sin(b-a))^2))+A*a1*(1+mu)*(sin(b-2*a)+sin(b))*sin(e)*k_s^2/(2*(1+mu*(sin(b-a))^2))-2*b_2*A*a1*(sin(b)+mu*sin(b-2*a))*cos(e)*k_s^2/(w*(1+mu*(sin(b-a))^2))+2*cos(b-a)*(2*b_1*a1-(1+mu)*a1*b_2)*c/(1+mu*(sin(b-a))^2)+(2*b_1*(cos(b-a)^2)-2*b_2*(1+mu))*d/(1+mu*(sin(b-a))^2);

t=linspace(0,t_f,tt)
f=@(t,x) [x(2);f1(x(1),x(3),x(2),x(4),x(5));x(4);f2(x(1),x(3),x(2),x(4),x(5));w];
[t,x]=ode45(f,t,[pi/180 0 pi/180 0 0 ]);
% [x(1)=theta_1, x(2)=d(theta_1)/dt, x(3)=theta_2, x(4)=d(theta_2)/dt, w] 
y1=x(740000:tt,1);
y2=x(740000:tt,2);
y3=x(740000:tt,3);
y4=x(740000:tt,4);
time=t(740000:tt);
[value,index]=max(y1);
%[value,index]=max(y4);
t_new=time(index);
tic;
  count=1;
for A= 0.5:0.005:3
 p=15000;
 t_f=p;
 int=[pi/180 0  pi/180 0 0 ];
 tt=1000;
for i=1:50
%step=0.1;
%tt=t_f/step;
k_s=1; %natural frequency of internal pendulum (sqrt(g/l))
f1=@(a,b,c,d,e) mu*sin(b-a)*(a1*cos(b-a)*c^2+d^2)/(a1*(1+mu*(sin(b-a))^2))+(mu*sin(2*b-a)-(2+mu)*sin(a))*k_s^2/(2*(1+mu*(sin(b-a))^2))+A*((2+mu)*sin(a)-mu*sin(2*b-a))*sin(e)*k_s^2/(2*(1+mu*(sin(b-a))^2))+b_1*A*(sin(2*b-a)-3*sin(a))*cos(e)*k_s^2/(w*(1+mu*(sin(b-a))^2))-2*b_1*(2-(cos(b-a))^2)*c/(1+mu*(sin(b-a))^2);
f2=@(a,b,c,d,e) -sin(b-a)*((1+mu)*a1*c^2+mu*cos(b-a)*d^2)/(1+mu*(sin(b-a))^2)-a1*(1+mu)*(sin(b-2*a)+sin(b))*(k_s^2)/(2*(1+mu*(sin(b-a))^2))+A*a1*(1+mu)*(sin(b-2*a)+sin(b))*sin(e)*k_s^2/(2*(1+mu*(sin(b-a))^2))-2*b_2*A*a1*(sin(b)+mu*sin(b-2*a))*cos(e)*k_s^2/(w*(1+mu*(sin(b-a))^2))+2*cos(b-a)*(2*b_1*a1-(1+mu)*a1*b_2)*c/(1+mu*(sin(b-a))^2)+(2*b_1*(cos(b-a)^2)-2*b_2*(1+mu))*d/(1+mu*(sin(b-a))^2);


t=linspace(0,t_f,tt)
f=@(t,x) [x(2);f1(x(1),x(3),x(2),x(4),x(5));x(4);f2(x(1),x(3),x(2),x(4),x(5));w];
[t,x]=ode45(f,t,int);
% [x(1)=theta_1, x(2)=d(theta_1)/dt, x(3)=] 
%y1=x(19900:tt,1);
%y2=x(19900:tt,2);
%y3=x(19900:tt,3);
%tim=t(19900:tt);


%plot(y1,y2)
A1(count)=A;
T1(count)=x(tt,1);
T2(count)=x(tt,2);
T3(count)=x(tt,3);
T4(count)=x(tt,4);
count=count+1;
%plot(A1,T,'.b')
%hold on
int=[x(tt,1) x(tt,2) x(tt,3) x(tt,4) x(tt,5)];
t_f=2*pi/w;
tt=20;
end
end
% filename=('bifurcation_sh1_b02_fixed_inc.mat')
% save(filename)
toc;
plot(A1,T3,'.r')   
    

[William Rose]

@Rebeka, I'm sorry but I do not have more to deevote to theis problem. I hope others will assist you.

Answer 2

[Lazaros Moysis]

The bifurcation diagram depicts changes in the dynamical behavior of the system, as one parameter changes. Changing a parameter of the system will change the dynamical behavior of the system. The bifurcation diagrams depicts exactly these changes. So by looking at it, we can understand for which parameter values the system is periodic with period-n, or chaotic. We do this usually by computing how many intersections the system's trajectory will have with a specific plane. If it is periodic, it will have a fixed amount of intersections, if it is chaotic it will have inumerable amount of intersections. The diagram depicts these points of intersection, so we can observe the behavior of the system, as well as how it transitions in and out of chaos as a parameter changes. The following two videos explain the theory behind the diagram, as well as the Matlab code to generate it. They will help you immensly in understanding the diagram.

https://www.youtube.com/watch?v=OXTS-d9HdB4

https://www.youtube.com/watch?v=O506G8DtmQk