绘制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/omega)/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):Tbushu*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)+epsY(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.','markersize',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.','markersize',5); % title('参数变化下的分岔图')

% xlabel('激励频率omega的变化') % ylabel('X值')

% clear;clc

% global m c k1 k3 F0 omega

% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@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,:)=lyapunovfun(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/omega)/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):Tbushu*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/omega)/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):Tbushu*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.','markersize',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/omega)/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):Tbushu*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.','markersize',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/omega)/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):Tbushu*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.','markersize',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/omega)/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):Tbushu*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.','markersize',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/omega)/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):Tbushu*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_integrator=@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,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp) % end

% figure,plot(kcanshu,lyapunovzhishu), % title('Lyapunov 动力学指数');

% xlabel('线性刚度参数k1,角频率为2rad/s'); ylabel('Lyapunov 指数');

% %

% clear,clc

% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@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,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp) % end

% figure,plot(kcanshu,lyapunovzhishu),set(gca,'fontsize',20); % title('Lyapunov 动力学指数','fontsize',20);

% xlabel('非线性刚度参数k3,角频率为2rad/s','fontsize',20); ylabel('Lyapunov 指数','fontsize',20);

本文来源：https://www.bwwdw.com/article/64oh.html