论坛 › 数学与非线性科学 › Duffing系统分岔程序

liliangbiao 发表于 2011-5-30 18:13

Duffing系统分岔程序

这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。因为我都不知道道理在哪里。但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。别照搬和迷恋别人的程序！
% % %%%% 绘制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=;
% tstart=0;Tbushu=600;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*150;
% tspan=;
% =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);
% =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=;
% YY=[];k=0;
% for c=range
%   k=k+1;
%   y0=;
%   tspan=;
%   =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,),'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=;
%   tspan=;
%   =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=;
%      tspan=;
%      =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=;
%      tspan=;
%      =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=;
%      tspan=;
%      =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=;
%      tspan=;
%      =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=,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=;
%   tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
%   tspan=;
%   =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
%         =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=;
%   tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
%   tspan=;
%   =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
%         =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=;
%   tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
%   tspan=;
%   =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
%         =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=;
%   tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
%   tspan=;
%   =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
%         =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=;
%   tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
%   tspan=;
%   =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
%         =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=;
    tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;
    tspan=;
    =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
      =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=,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=,ioutp=10,
% global m c k1 k3 F0 omega
% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;
% kcanshu=;
% 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 发表于 2011-6-1 10:25

感谢楼主！1！！

chenshumeionly 发表于 2011-9-25 23:08

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

liliangbiao 发表于 2011-9-26 15:21

本帖最后由 liliangbiao 于 2011-9-26 15:21 编辑

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

菜肴全 发表于 2012-4-2 13:09

感谢楼主啊

paopaotaiqiang 发表于 2012-4-20 21:51

感谢，学习中

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

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

ppsqtjn 发表于 2012-10-3 04:03

顶！！！！！！！











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

lovelycoco 发表于 2012-12-11 17:13

学习中，多谢楼主分享！

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

轩小轩 发表于 2012-9-7 15:10 static/image/common/back.gif
楼主的那个lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,iout ...
这个应该是楼主自己编写的李指数计算程序，只能联系楼主，给楼主发短消息看看

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

gghhjj 发表于 2012-12-13 04:49 static/image/common/back.gif
这个应该是楼主自己编写的李指数计算程序，只能联系楼主，给楼主发短消息看看

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

gghhjj 发表于 2012-12-17 09:21

轩小轩 发表于 2012-12-16 23:12 static/image/common/back.gif
哦，楼主不好联系啊，你有李雅普诺夫程序吗

搜索论坛，有好多

x589 发表于 2014-12-24 23:23

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

x589 发表于 2014-12-24 23:24

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

查看完整版本: Duffing系统分岔程序