绘制Duffing振子的分叉图的程序Word格式文档下载.docx

绘制Duffing振子的分叉图的程序Word格式文档下载.docx

《绘制Duffing振子的分叉图的程序Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《绘制Duffing振子的分叉图的程序Word格式文档下载.docx（9页珍藏版）》请在冰点文库上搜索。

%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振子的分叉图的程序

%globalmck1k3F0omega;

F0=12;

%range=[0.01:

0.01:

1];

%YY=[];

k=0;

%forc=range

%k=k+1;

%y0=[3,4];

%tspan=[0:

200];

%[t,Y]=ode45（'

tspan,y0）;

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;

%ifj>

1

%plot（c,YY（k,[1:

j-1]）,'

%holdon;

%index（k）=j-1;

%xlabel（'

c'

%ylabel（'

xmax'

%title（'

dafinbifurcationdiagram'

%%%绘制分岔图的程序

%clear,clc

%ccanshu=0.01:

1;

%fork=1:

100

%c=ccanshu（k）

0.01*2*pi:

500];

%dis=zeros（50,1）;

%velo=zeros（50,1）;

50

%dis（i,1）=y（100*（i+20）,1）;

%velo（i,1）=y（100*（i+20）,2）;

%Dismatrix（k,:

）=dis'

;

%figure,plot（ccanshu,Dismatrix,'

3）;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%线性参数k1的变化产生的分岔图

%kcanshu=0.01:

2;

200

%k1=kcanshu（k）

%plot（kcanshu,Dismatrix,'

参数变化下的分岔图'

）

线性刚度参数k1的变化'

X值'

%%非线性参数k3的变化产生的分岔图

%k3=kcanshu（k）

非线性参数k3的变化'

%%激励参数F0变化产生的分岔图

%F0canshu=0.1:

0.1:

20;

%F0=F0canshu（k）

%plot（F0canshu,Dismatrix,'

激励参数F0的变化'

%%%激励频率omega变化产生的分岔图

%omegacanshu=0.1:

10;

%omega=omegacanshu（k）

0.01*2*pi/omega:

%dis（i,1）=y（round（100*omega*（i+20））,1）;

%velo（i,1）=y（round（100*omega*（i+20））,2）;

%plot（omegacanshu,Dismatrix,'

激励频率omega的变化'

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

%ystart=[340],ioutp=10,

%omega=omegacanshu（1,k）,lyapunovzhishu（k,:

）=lyapunovfun（n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp）

%figure,plot（omegacanshu,lyapunovzhishu）,

Lyapunov动力学指数'

激励频率omega变化'

ylabel（'

Lyapunov指数'

Tbushu=100;

tend=（2*pi/omega）*200;

%ifk==1

%[maxvalue,indices]=max（abs（TData））;

%plot（ccanshu,Dismatrix,'

c=0.4;

k1=-1;

F0=3;

omega=2;

%x0=[2;

0];

%set（gca,'

fontsize'

20）;

随参数变化的分岔图'

随阻尼参数c变化'

%k1canshu=-1:

0.99;

%k1=k1canshu（k）

%figure,plot（k1canshu,Dismatrix,'

%axis（[-1,1,-1,4]）

随线性刚度参数k1的变化'

%k3canshu=0.01:

%k3=k3canshu（k）

%figure,plot（k3canshu,Dismatrix,'

随非线性刚度参数k3的变化'

%figure,plot（F0canshu,Dismatrix,'

随外界激励幅值F0的变化'

%%激励频率omega变化产生的分岔图

clear;

globalmck1k3F0omega

m=1;

omegacanshu=0.1:

fork=1:

omega=omegacanshu（k）

x0=[3;

tstart=0;

tspan=[tstart:

[t,y]=ode45（'

count=find（t>

Y=y（count,:

TData=Y（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:

dis（i,1）=Y（Tbushu*（i-1）+indices,1）;

velo（i,1）=Y（Tbushu*（i-1）+indices,2）;

Dismatrix（k,:

end

figure,plot（omegacanshu,Dismatrix,'

set（gca,'

title（'

xlabel（'

随外界激励频率omega的变化'

%%%%%%%%%%%%%%%%%%%%%%

%ys