绘制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)); % 去掉前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% 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% 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(t100);% Y=Y(count,:);% % 画x的分岔图。% j=1;% n=length(Y(:,1));% for i=2:n-1% if Y(i-1,1)+epsY(i,1) Y(i,1)Y(i+1,1)+eps %简单的取出局部最大值。% YY(k,j)=Y(i,1); %使最大值计数个数自动增加% j=j+1;% end% end% if j1% plot(c,YY(k,[1:j-1]),b.,markersize,5);% end% hold on;% index(k)=j-1;% end% xlabel(c);% ylabel(x max);% title(dafin bifurcation diagram);% % % 绘制分岔图的程序% 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% 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% dis(i,1)=y(100*(i+20),1);% velo(i,1)=y(100*(i+20),2);% end% Dismatrix(k,:)=dis;% end% figure,plot(ccanshu,Dism