繪制Duffing振子地分叉圖地程序

實用標準文案文檔大全這些程序思想有些可能不正確，有問題，自己改進，我不再負責對這些程序解釋。因為我都不知道道理在哪里。但是期望您能在程序的提示下，進一步的做改進或者改正，以期獲得更為精確的結(jié)果。別照搬和迷戀別人的程序！

%%%%%%繪制Duffing振子的龐加萊截面圖的程序

%%buchang:已知激勵下步長數(shù)值的大小,

%%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個周期的激勵時間以消除瞬態(tài)響應的影響

%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);%使最大值計數(shù)個數(shù)自動增加

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

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

%%線性參數(shù)k1的變化產(chǎn)生的分岔圖

%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('參數(shù)變化下的分岔圖')

%xlabel('線性剛度參數(shù)k1的變化')

%ylabel('X值')

%%非線性參數(shù)k3的變化產(chǎn)生的分岔圖

%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('參數(shù)變化下的分岔圖')

%xlabel('非線性參數(shù)k3的變化')

%ylabel('X值')

%%激勵參數(shù)F0變化產(chǎn)生的分岔圖

%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('參數(shù)變化下的分岔圖')

%xlabel('激勵參數(shù)F0的變化')

%ylabel('X值')

%%%激勵頻率omega變化產(chǎn)生的分岔圖

%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('參數(shù)變化下的分岔圖')

%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動力學指數(shù)');

%xlabel('激勵頻率omega變化');ylabel('Lyapunov指數(shù)');

%%%繪制分岔圖的程序

%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個周期的激勵時間以消除瞬態(tài)響應的影響

%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個周期的激勵時間以消除瞬態(tài)響應的影響

%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('隨參數(shù)變化的分岔圖','fontsize',20);

%xlabel('隨阻尼參數(shù)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個周期的激勵時間以消除瞬態(tài)響應的影響

%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('隨參數(shù)變化的分岔圖','fontsize',20);

%xlabel('隨線性剛度參數(shù)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個周期的激勵時間以消除瞬態(tài)響應的影響

%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('隨參數(shù)變化的分岔圖','fontsize',20);

%xlabel('隨非線性剛度參數(shù)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個周期的激勵時間以消除瞬態(tài)響應的影響

%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('隨參數(shù)變化的分岔圖','fontsize',20);

%xlabel('隨外界激勵幅值F0的變化','fontsize',20);

%%激勵頻率omega變化產(chǎn)生的分岔圖

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個周期的激勵時間以消除瞬態(tài)響應的影響

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