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

实用标准文案文档大全这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。因为我都不知道道理在哪里。但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。别照搬和迷恋别人的程序！

%%%%%%绘制Duffing振子的庞加莱截面图的程序

%%buchang:已知激励下步长数值的大小,

%%tend程序仿真达到150个激励周期的总时间,

%clear;clc

%globalmck1k3F0omega

%

%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);

%fori=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

%globalmck1k3F0omega;

%m=1;k1=0;k3=1;omega=1;F0=12;

%range=[0.01:0.01:1];

%YY=[];k=0;

%forc=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));

%fori=2:n-1

%ifY(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

%ifj>1

%plot(c,YY(k,[1:j-1]),'b.','markersize',5);

%end

%holdon;

%index(k)=j-1;

%end

%xlabel('c');

%ylabel('xmax');

%title('dafinbifurcationdiagram');

%%%绘制分岔图的程序

%clear,clc

%globalmck1k3F0omega

%

%m=1;c=0.1;k1=0;k3=1;omega=1;F0=12;

%ccanshu=0.01:0.01:1;

%fork=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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.1;k3=1;omega=1;F0=12;

%kcanshu=0.01:0.01:2;

%fork=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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.1;k1=0;omega=1;F0=12;

%kcanshu=0.01:0.01:2;

%fork=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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.1;k1=0;k3=1;omega=1;

%F0canshu=0.1:0.1:20;

%fork=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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.1;k1=0;k3=1;F0=12;

%omegacanshu=0.1:0.1:10;

%fork=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);

%fori=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

%globalmck1k3F0omega

%n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@ode45,tstart=0,stept=0.5,tend=200,

%ystart=[340],ioutp=10,

%m=1;c=0.1;k1=0;F0=12;k3=1;

%omegacanshu=0.1:0.1:10;

%fork=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

%globalmck1k3F0omega

%

%m=1;c=0.1;k1=0;k3=1;omega=1;F0=12;

%ccanshu=0.01:0.01:1;

%fork=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);

%ifk==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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;

%ccanshu=0.01:0.01:1;

%fork=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);

%ifk==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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;

%k1canshu=-1:0.01:0.99;

%fork=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);

%ifk==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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;

%k3canshu=0.01:0.01:1;

%fork=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);

%ifk==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);

%fori=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

%globalmck1k3F0omega

%m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;

%F0canshu=0.1:0.1:10;

%fork=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);

%ifk==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);

%fori=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

globalmck1k3F0omega

m=1;c=0.1;k1=0;k3=1;F0=12;

omegacanshu=0.1:0.1:10;

fork=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);

ifk==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);

fori=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