绘制Duffing振子的分叉图的程序

绘制Duffing振子的分叉图的程序
这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。

因为我都不知道道理在哪里。

但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。

别照搬和迷恋别人的程序！
% % %%%% 绘制Duffing振子的庞加莱截面图的程序
% % buchang:已知激励下步长数值的大小, % % tend程序仿真达到150个激励周期的总时间,
% clear;clc
% global m c k1 k3 F0 omega
%
%
m=1;c=0.1;k1=0;k3=1;omega=1;F0=12 % x0=[3;4];
%
tstart=0;Tbushu=600;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*150; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% [maxvalue,indices]=max(abs(TData)) %
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% figure,plot(dis,velo,'b.','markersize',5);
% %%%% 绘制Duffing振子的分叉图的程序
% clear;clc
% global m c k1 k3 F0 omega;
% m=1;k1=0;k3=1;omega=1;F0=12;
% range=[0.01:0.01:1];
% YY=[];k=0;
% for c=range
% k=k+1;
% y0=[3,4];
% tspan=[0:0.01:200];
% [t,Y]=ode45('dafin3',tspan,y0);
% count=find(t>100);
% Y=Y(count,:);
% % 画x的分岔图。

% j=1;
% n=length(Y(:,1));
% for i=2:n-1
% if Y(i-1,1)+eps<Y(i,1) & Y(i,1)>Y(i+1,1)+eps %简单的取出局部最大值。

%
YY(k,j)=Y(i,1);
%使最大值计数个数自动增加
% j=j+1;
% end
% end
% if j>1
%
plot(c,YY(k,[1:j-1]),'b.','markersize',5); % end
% hold on;
% index(k)=j-1;
% end
% xlabel('c');
% ylabel('x max');
% title('dafin bifurcation diagram');
% % % 绘制分岔图的程序
% clear,clc
% global m c k1 k3 F0 omega
%
%
m=1;c=0.1;k1=0;k3=1;omega=1;F0=12;
% ccanshu=0.01:0.01:1;
% for k=1:100
% c=ccanshu(k)
% x0=[3;4];
% tspan=[0:0.01*2*pi:500];
% [t,y]=ode45('dafin3',tspan,x0); % dis=zeros(50,1);
% velo=zeros(50,1);
% for i=1:50
% dis(i,1)=y(100*(i+20),1);
% velo(i,1)=y(100*(i+20),2); % end
% Dismatrix(k,:)=dis';
% end
%
figure,plot(ccanshu,Dismatrix,'b.','marker size',3); %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% %%%%%
% %线性参数k1的变化产生的分岔图
% clear;clc
% global m c k1 k3 F0 omega
% m=1;c=0.1;k3=1;omega=1;F0=12; % kcanshu=0.01:0.01:2;
% for k=1:200
% k1=kcanshu(k)
% x0=[3;4];
% tspan=[0:0.01*2*pi:500];
% [t,y]=ode45('dafin3',tspan,x0); % dis=zeros(50,1);
% velo=zeros(50,1);
% for i=1:50
% dis(i,1)=y(100*(i+20),1);
% velo(i,1)=y(100*(i+20),2); % end
% Dismatrix(k,:)=dis';
% end
%
plot(kcanshu,Dismatrix,'b.','markersize',5 );
% title('参数变化下的分岔图')
% xlabel('线性刚度参数k1的变化')
% ylabel('X值')
% %非线性参数k3的变化产生的分岔图
% clear;clc
% global m c k1 k3 F0 omega
% m=1;c=0.1;k1=0;omega=1;F0=12; % kcanshu=0.01:0.01:2;
% for k=1:200
% k3=kcanshu(k)
% x0=[3;4];
% tspan=[0:0.01*2*pi:500];
% [t,y]=ode45('dafin3',tspan,x0); % dis=zeros(50,1);
% velo=zeros(50,1);
% for i=1:50
% dis(i,1)=y(100*(i+20),1);
% velo(i,1)=y(100*(i+20),2); % end
% Dismatrix(k,:)=dis';
% end
%
plot(kcanshu,Dismatrix,'b.','markersize',5 );
% title('参数变化下的分岔图')
% xlabel('非线性参数k3的变化')
% ylabel('X值')
% %激励参数F0变化产生的分岔图
% clear;clc
% global m c k1 k3 F0 omega
% m=1;c=0.1;k1=0;k3=1;omega=1; % F0canshu=0.1:0.1:20;
% for k=1:200
% F0=F0canshu(k)
% x0=[3;4];
% tspan=[0:0.01*2*pi:500];
% [t,y]=ode45('dafin3',tspan,x0); % dis=zeros(50,1);
% velo=zeros(50,1);
% for i=1:50
% dis(i,1)=y(100*(i+20),1); % velo(i,1)=y(100*(i+20),2); % end
% Dismatrix(k,:)=dis';
% end
%
plot(F0canshu,Dismatrix,'b.','markersize', 5);
% title('参数变化下的分岔图')
% xlabel('激励参数F0的变化')
% ylabel('X值')
% % %激励频率omega变化产生的分岔图% clear;clc
% global m c k1 k3 F0 omega
% m=1;c=0.1;k1=0;k3=1;F0=12;
% omegacanshu=0.1:0.1:10;
% for k=1:100
% omega=omegacanshu(k)
% x0=[3;4];
%
tspan=[0:0.01*2*pi/omega:500];
% [t,y]=ode45('dafin3',tspan,x0); % dis=zeros(50,1);
% velo=zeros(50,1);
% for i=1:50
%
dis(i,1)=y(round(100*omega*(i+20)),1);
%
velo(i,1)=y(round(100*omega*(i+20)),2 );
% end
% Dismatrix(k,:)=dis';
% end
%
plot(omegacanshu,Dismatrix,'b.','marker size',5);
% title('参数变化下的分岔图')
% xlabel('激励频率omega的变化')
% ylabel('X值')
% clear;clc
% global m c k1 k3 F0 omega
%
n=3,rhs_ext_fcn=@dafin_ext2,fcn_integ rator=@ode45,tstart=0,stept=0.5,tend= 200,
% ystart=[3 4 0],ioutp=10,
% m=1;c=0.1;k1=0;F0=12;k3=1;
% omegacanshu=0.1:0.1:10;
% for k=1:100
% omega = omegacanshu(1,k),lyapunovzhishu(k,:)=l yapunovfun(n,rhs_ext_fcn,fcn_integrator ,tstart,stept,tend,ystart,ioutp)
% end
%
figure,plot(omegacanshu,lyapunovzhishu ),
% title('Lyapunov 动力学指数');
% xlabel('激励频率omega变化'); ylabel('Lyapunov 指数');
% % % 绘制分岔图的程序
% clear;clc
% global m c k1 k3 F0 omega
%
%
m=1;c=0.1;k1=0;k3=1;omega=1;F0=12; % ccanshu=0.01:0.01:1;
% for k=1:100
% c=ccanshu(k)
% x0=[3;4];
%
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% if k==1
%
[maxvalue,indices]=max(abs(TData));
% end
%
pointnumber=round((tend-2*pi/omega*
40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% Dismatrix(k,:)=dis';
% end
%
plot(ccanshu,Dismatrix,'b.','markersize',3 );
% % % 绘制分岔图的程序
% clear,clc
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2; % ccanshu=0.01:0.01:1;
% for k=1:100
% c=ccanshu(k)
% x0=[2;0];
%
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% if k==1
%
[maxvalue,indices]=max(abs(TData));
% end
%
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% Dismatrix(k,:)=dis';
% end
%
figure,plot(ccanshu,Dismatrix,'b.','marker size',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20); % xlabel('随阻尼参数c变化','fontsize',20);
% % % 绘制分岔图的程序
% clear,clc
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;
% k1canshu=-1:0.01:0.99;
% for k=1:200
% k1=k1canshu(k)
% x0=[2;0];
%
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% if k==1
%
[maxvalue,indices]=max(abs(TData));
% end
%
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% Dismatrix(k,:)=dis';
% end
%
figure,plot(k1canshu,Dismatrix,'b.','mark ersize',3);
% axis([-1,1,-1,4])
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20); % xlabel('随线性刚度参数k1的变化','fontsize',20);
% % % 绘制分岔图的程序
% clear,clc
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2; % k3canshu=0.01:0.01:1;
% for k=1:100
% k3=k3canshu(k)
% x0=[2;0];
%
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% if k==1
%
[maxvalue,indices]=max(abs(TData)); % end
%
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% Dismatrix(k,:)=dis';
% end
%
figure,plot(k3canshu,Dismatrix,'b.','mark ersize',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20); % xlabel('随非线性刚度参数k3的变化
','fontsize',20);
% % % 绘制分岔图的程序
% clear,clc
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2; % F0canshu=0.1:0.1:10;
% for k=1:100
% F0=F0canshu(k)
% x0=[2;0];
%
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200; % tspan=[tstart:buchang:tend];
% [t,y]=ode45('dafin3',tspan,x0);
%
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
% Y=y(count,:);
%
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
% if k==1
%
[maxvalue,indices]=max(abs(TData)); % end
%
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
% dis=zeros(pointnumber,1);
% velo=zeros(pointnumber,1);
% for i=1:pointnumber
%
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
%
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% end
% Dismatrix(k,:)=dis';
% end
%
figure,plot(F0canshu,Dismatrix,'b.','mark ersize',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20); % xlabel('随外界激励幅值F0的变化','fontsize',20);
% %激励频率omega变化产生的分岔图clear;clc
global m c k1 k3 F0 omega
m=1;c=0.1;k1=0;k3=1;F0=12; omegacanshu=0.1:0.1:10;
for k=1:100
omega=omegacanshu(k)
x0=[3;4];
tstart=0;Tbushu=100;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*200;
tspan=[tstart:buchang:tend];
[t,y]=ode45('dafin3',tspan,x0);
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬态响应的影响
Y=y(count,:);
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);
if k==1
[maxvalue,indices]=max(abs(TData));
end
pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;
dis=zeros(pointnumber,1);
velo=zeros(pointnumber,1);
for i=1:pointnumber
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
end
Dismatrix(k,:)=dis';
end
figure,plot(omegacanshu,Dismatrix,'b.','
markersize',3);
set(gca,'fontsize',20);
title('随参数变化的分岔图','fontsize',20); xlabel('随外界激励频率omega的变化','fontsize',20);
% %%%%%%%%%%%%%%%% %%%%%
% clear,clc
%
n=3,rhs_ext_fcn=@dafin_ext2,fcn_integ rator=@ode45,tstart=0,stept=0.5,tend= 200,
% ystart=[2 0 0],ioutp=10,
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2; % kcanshu=-1:0.01:0.99;
% for k=1:200
%
k1=kcanshu(1,k),lyapunovzhishu(k,:)=ly apunovfun(n,rhs_ext_fcn,fcn_integrator,
% end
% figure,plot(kcanshu,lyapunovzhishu), % title('Lyapunov 动力学指数');
% xlabel('线性刚度参数k1,角频率为2rad/s'); ylabel('Lyapunov 指数');
%
%
% clear,clc
%
n=3,rhs_ext_fcn=@dafin_ext2,fcn_integ rator=@ode45,tstart=0,stept=0.5,tend= 200,
% ystart=[2 0 0],ioutp=10,
% global m c k1 k3 F0 omega
%
m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2; % kcanshu=[0.1:0.1:10,10:10:1000];
% for k=1:200
%
k3=kcanshu(1,k),lyapunovzhishu(k,:)=ly apunovfun(n,rhs_ext_fcn,fcn_integrator,
% end
%
figure,plot(kcanshu,lyapunovzhishu),set( gca,'fontsize',20);
% title('Lyapunov 动力学指数','fontsize',20);
% xlabel('非线性刚度参数k3,角频率为2rad/s','fontsize',20); ylabel('Lyapunov 指数','fontsize',20);。