Bifurcation diagram calcullation for ODE-system.
x_label - the paramter of the x label
rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
This function must include RHS of ODE-system coupled with
variational equation (n items of linearized systems, see Example).
fcn_integrator - handle of ODE integrator function, for example: @ode45
tstart - start values of independent value (time t)
stept - step on t-variable for Gram-Schmidt renormalization procedure.
tend - finish value of time
ystart - start point of trajectory of ODE system.
y_dest - intercept quadrant
cut_start - remove the previous transient state
X - x label
Y - the ordinate of the bifurcation diagram
clc;clear all;close all;
tic;
W1 = [];
W2 = [];
W3 = [];
W4 = [];
d =-1:0.01:1;
N = length(d);
matlabpool local 4
parfor i = 1:N
% for i = 1:N
fprintf('i = %d\n',i);
[x,y]= bifurcation(d(i),@(t,y)JM1(t,y),@ode23,0,1e-2,1000,[1 0 0 d(i)],3,5000);
W1 = [W1;x];
W2 = [W2;y];
[x,y]= bifurcation(d(i),@(t,y)JM1(t,y),@ode23,0,1e-2,1000,[-1 0 0 d(i)],3,5000);
W3 = [W3;x];
W4 = [W4;y];
end
save bifur
matlabpool close
toc;
figure(1);
plot(W1,W2,'b.','markersize',0.5); hold on
plot(W3,W4,'r.','markersize',0.5);
function dy=JM1(t,y)
dy=zeros(3,1);
a=0.6;c=2;a1=1;b1=0.1;b=1.3;
dy(1)=abs(y(2))-b;
dy(2)=(a1-b1*abs(y(4)))*y(3);
dy(3)=abs(y(1))-y(2)-a*y(3)-c;
dy(4)=y(3);