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

1、绘制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;% ts

2、tart=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-

3、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

4、;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(t100);% 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% en

5、d% if j1% 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:

6、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;

7、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,

8、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

9、: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:20

10、0% 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的变化)% ylab

11、el(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*

12、(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,% yst

13、art=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变化); yl

14、abel(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

15、=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

16、)=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;Tbus

17、hu=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(te

18、nd-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(随参

19、数变化的分岔图,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:ten

20、d;% 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(pointnu

21、mber,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);% % % 绘制分

22、岔图的程序% 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);

23、% 去掉前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

24、,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

25、;% 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):Tbush

26、u*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% figu

27、re,plot(F0canshu,Dismatrix,b.,markersize,3);% set(gca,fontsize,20);% title(随参数变化的分岔图,fontsize,20);% xlabel(随外界激励幅值F0的变化,fontsize,20);% %激励频率omega变化产生的分岔图clear;clcglobal m c k1 k3 F0 omegam=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;buc

28、hang=(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)/b

29、uchang/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;endfigure,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,% ys