Duffing系统分岔程序

liliangbiao

这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。因为我都不知道道理在哪里。但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。别照搬和迷恋别人的程序！
% % %%%% 绘制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(t>100);
%     Y=Y(count,:);
% %     画x的分岔图。
%     j=1;
%     n=length(Y(:,1));
%     for i=2:n-1
%         if Y(i-1,1)+eps<Y(i,1) & Y(i,1)>Y(i+1,1)+eps  %简单的取出局部最大值。
%             YY(k,j)=Y(i,1);                           %使最大值计数个数自动增加
%             j=j+1;
%         end
%     end
%     if j>1
%         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,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;
%  kcanshu=0.01:0.01:2;
%  for k=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);
%      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
%  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
%     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';
%  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;
%  F0canshu=0.1:0.1:20;
%  for k=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);
%      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
%  plot(F0canshu,Dismatrix,'b.','markersize',5);
% title('参数变化下的分岔图')
% xlabel('激励参数F0的变化')
% ylabel('X值')

% % %激励频率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
%     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
%          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
% 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
% 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
% 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];
%     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
%         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
% 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
%     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));
%     end
%     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
%     Dismatrix(k,:)=dis';
% end
% figure,plot(ccanshu,Dismatrix,'b.','markersize',3);
% set(gca,'fontsize',20);
% title('随参数变化的分岔图','fontsize',20);
% xlabel('随阻尼参数c变化','fontsize',20);

% % % 绘制分岔图的程序
% 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个周期的激励时间以消除瞬态响应的影响
%     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
%         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
% 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
%     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);
%     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
%         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
% 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
%     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
%         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
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
    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);
    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
        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,
% ystart=[2 0 0],ioutp=10,
% global m c k1 k3 F0 omega
% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;
% kcanshu=-1:0.01:0.99;
% for k=1:200
% k1=kcanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp)
% end
% figure,plot(kcanshu,lyapunovzhishu),
% title('Lyapunov 动力学指数');
% xlabel('线性刚度参数k1,角频率为2rad/s'); ylabel('Lyapunov 指数');
%
%
% clear,clc
% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@ode45,tstart=0,stept=0.5,tend=200,
% ystart=[2 0 0],ioutp=10,
% global m c k1 k3 F0 omega
% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;
% kcanshu=[0.1:0.1:10,10:10:1000];
% for k=1:200
% k3=kcanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp)
% end
% figure,plot(kcanshu,lyapunovzhishu),set(gca,'fontsize',20);
% title('Lyapunov 动力学指数','fontsize',20);
% xlabel('非线性刚度参数k3,角频率为2rad/s','fontsize',20); ylabel('Lyapunov 指数','fontsize',20);

hmlz

感谢楼主！1！！

chenshumeionly

楼主漂亮的分岔图都是用这样的程序做的么？

liliangbiao

不是的，但是思想几乎都一样，但是细节处理不同。只要读懂了我最初给的思想，就能编写。不需要向别人求程序。

菜肴全

感谢楼主啊

paopaotaiqiang

感谢，学习中

轩小轩

楼主的那个lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp)
这个函数有吗，可不可以传上来学习一下，谢谢！

ppsqtjn

顶！！！！！！！











lagunamoon  汗巾 電子防潮箱 男子漢  wikipps.hk

lovelycoco

学习中，多谢楼主分享！

gghhjj

轩小轩 发表于 2012-9-7 15:10

                               
登录/注册后可看大图

楼主的那个lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,iout ...

这个应该是楼主自己编写的李指数计算程序，只能联系楼主，给楼主发短消息看看

轩小轩

gghhjj 发表于 2012-12-13 04:49

                               
登录/注册后可看大图

这个应该是楼主自己编写的李指数计算程序，只能联系楼主，给楼主发短消息看看

哦，楼主不好联系啊，你有李雅普诺夫程序吗

gghhjj

轩小轩 发表于 2012-12-16 23:12

                               
登录/注册后可看大图

哦，楼主不好联系啊，你有李雅普诺夫程序吗

搜索论坛，有好多

x589

我也想要，恩能够共享下不  

x589

贴出来吧，楼主  这么好的东西  咱急需哦