绘制Duffing振子的分叉图的程序.docx
《绘制Duffing振子的分叉图的程序.docx》由会员分享，可在线阅读，更多相关《绘制Duffing振子的分叉图的程序.docx（9页珍藏版）》请在冰豆网上搜索。
绘制Duffing振子的分叉图的程序
绘制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）+epsY（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,'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
升级会员 