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

75 views (last 30 days)
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.

Answers (1)

William Rose
William Rose on 10 Nov 2022
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.
  3 Comments
Rebeka
Rebeka on 14 Nov 2022
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')

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!