第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
绘制Duffing振子的分叉图的程序.doc
这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。因为我都
不知道道理在哪里。但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更
为精确的结果。别照搬和迷恋别人的程序！
% % %%%% 绘制 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];
% 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));
应的影响
% 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);
% for i=1:pointnumber
%
%
% end
% figure,plot(dis,velo,'b.','markersize',5);
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
% 去掉前40个周期的激励时间以消除瞬态响
% %%%% 绘制 Duffing 振子的分叉图的程序
% clear;clc
% global m c k1 k3 F0 omega;
% m=1;k1=0;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(t>100);
%
%
Y=Y(count,:);
% %
%
画 x 的分岔图。
j=1;
if Y(i-1,1)+epsY(i+1,1)+eps %简单的取出局部最大值。
%使最大值计数个数自动增
plot(c,YY(k,[1:j-1]),'b.','markersize',5);
j=j+1;
end
end
if j>1
YY(k,j)=Y(i,1);
n=length(Y(:,1));
for i=2:n-1
%
%
%
%
加
%
%
%
%
%
%
%
%
% end
% xlabel('c');
% ylabel('x max');
% title('dafin bifurcation diagram');
end
hold on;
index(k)=j-1;
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);
for i=1:50
% % % 绘制分岔图的程序
% 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
%
%
%
%
%
%
%
%
%
%
%
% 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;omega=1;F0=12;
dis(i,1)=y(100*(i+20),1);
velo(i,1)=y(100*(i+20),2);
end
Dismatrix(k,:)=dis';
end
Dismatrix(k,:)=dis';
dis(i,1)=y(100*(i+20),1);
velo(i,1)=y(100*(i+20),2);
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
% kcanshu=0.01:0.01:2;
% for k=1:200
%
%
%
%
%
%
%
%
%
%
%
% end
% plot(kcanshu,Dismatrix,'b.','markersize',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
%
%
%
%
%
%
%
%
%
%
%
% 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;
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:50
dis(i,1)=y(100*(i+20),1);
velo(i,1)=y(100*(i+20),2);
end
Dismatrix(k,:)=dis';
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
% F0canshu=0.1:0.1:20;
% for k=1:200
%
%
%
%
%
%
%
%
%
%
%
% end
% plot(F0canshu,Dismatrix,'b.','markersize',5);
% title('参数变化下的分岔图')
% xlabel('激励参数 F0的变化')
% ylabel('X 值')
dis(i,1)=y(100*(i+20),1);
velo(i,1)=y(100*(i+20),2);
end
Dismatrix(k,:)=dis';
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
% % %激励频率 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
%
%
%
%
%
%
%
%
%
%
%
% end
% plot(omegacanshu,Dismatrix,'b.','markersize',5);
% title('参数变化下的分岔图')
% xlabel('激励频率 omega 的变化')
% ylabel('X 值')
end
Dismatrix(k,:)=dis';
dis(i,1)=y(round(100*omega*(i+20)),1);
velo(i,1)=y(round(100*omega*(i+20)),2);
% 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,
% ystart=[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
%
=
omegacanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,st
ept,tend,ystart,ioutp)
% end
% figure,plot(omegacanshu,lyapunovzhishu),
% title('Lyapunov 动力学指数');
% xlabel('激励频率 omega 变化'); ylabel('Lyapunov 指数');
omega
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个周期的激励时间以消除瞬
% % % 绘制分岔图的程序
% 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
%
%
%
%
%
%
态响应的影响
%
%
%
%
%
%
%
%
%
%
%
%
%
% end
% plot(ccanshu,Dismatrix,'b.','markersize',3);
end
Dismatrix(k,:)=dis';
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
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
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);
if k==1
[maxvalue,indices]=max(abs(TData));
% % % 绘制分岔图的程序
% 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
%
%
%
%
%
%
态响应的影响
%
%
%
%
%
%
%
%
%
%
%
%
%
% end
% figure,plot(ccanshu,Dismatrix,'b.','markersize',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20);
% xlabel('随阻尼参数 c 变化','fontsize',20);
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
end
Dismatrix(k,:)=dis';
end
pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;
dis=zeros(pointnumber,1);
velo=zeros(pointnumber,1);
for i=1:pointnumber
% % % 绘制分岔图的程序
% 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:tend];
[t,y]=ode45('dafin3',tspan,x0);
count=find(t>(2*pi/omega*40));
% 去掉前40个周期的激励时间以消除瞬
end
pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;
dis=zeros(pointnumber,1);
velo=zeros(pointnumber,1);
for i=1:pointnumber
Y=y(count,:);
TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);
if k==1
[maxvalue,indices]=max(abs(TData));
%
态响应的影响
%
%
%
%
%
%
%
%
%
%
%
%
%
% end
% figure,plot(k1canshu,Dismatrix,'b.','markersize',3);
% axis([-1,1,-1,4])
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20);
% xlabel('随线性刚度参数 k1的变化','fontsize',20);
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
end
Dismatrix(k,:)=dis';
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个周期的激励时间以消除瞬
% % % 绘制分岔图的程序
% 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
%
%
%
%
%
%
态响应的影响
%
%
%
%
%
%
%
%
%
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
end
Dismatrix(k,:)=dis';
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
%
%
%
%
% 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;
% F0canshu=0.1:0.1:10;
% for k=1:100
%
%
%
%
%
%
态响应的影响
%
%
%
%
%
%
%
%
%
%
%
%
%
% end
% figure,plot(F0canshu,Dismatrix,'b.','markersize',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20);
% xlabel('随外界激励幅值 F0的变化','fontsize',20);
dis(i,1)=Y(Tbushu*(i-1)+indices,1);
velo(i,1)=Y(Tbushu*(i-1)+indices,2);
end
Dismatrix(k,:)=dis';
% %激励频率 omega 变化产生的分岔图
clear;clc
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);
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
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
