非自治系统分岔图绘制实例——有疑问的都到这里面来

octopussheng

近期关于分岔图计算的疑问贴比较多，这里以duffing系统为例，把非自治系统分岔图绘制的问题贴出来，有问题的都可以参考该例！
duffing系统：
function dx=duffing(t,X)
global F wd;
r=0.25;
x=X(1);
y=X(2);
psi=X(3);
dx=zeros(3,1);
dx(1)=y;
dx(2)=-r*y-x^3+F*cos(psi);
dx(3)=wd;

分岔图绘制程序：
function duffing_bifur_F
clear;
global F wd;
wd=1.0;
range=[0.1:0.001:0.35];
period=2*pi/wd;
k=0;
YY1=[];
step=2*pi/100;  %步长。
for F=range
    y0=[0 0.001 0];
    F
    k=k+1;
    % discard the first 60 periodic data;
    %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
    tspan=[0:step:60*period];
    [t,Y]=ode45(@duffing,tspan,y0);
    y0=Y(end,:);
    j=1;
    for i=60:200
        tspan=[i*period:step:(i+1)*period];
        [t,Y]=ode45(@duffing,tspan,y0);
        YY1(k,j)=Y(end,1);   % get the omega data from every period end
        j=j+1;               %取出每一个周期内的第一个解的最后一个值。
        y0=Y(end,:);
    end
end
bifdata=YY1(:,end-51:end);
plot(range,bifdata,'k.','markersize',1);


下面的图是其中计算得到的一个结果！

liliangbiao

你再换个非自治系统试验一下，看看这个程序是否合适？

octopussheng

合适的，我试过了！

ruichard101

有几处看不明白:loveliness:
yo是不是表示x(1),X(2),X(3)的初值呀?
YY1=[],  YY1(k  j),   YY1(:,end-51:end)又分别表示的是什么意思呢?

liliangbiao

这个程序 我原来贴在我的个人空间里面的，可是我没有做很多验证，不知道合适不？
还希望大家多多验证其他的系统试验一下！
回答上边的疑问：
1. yo是表示x(1),X(2),X(3)的初值！
2. YY1=[],  定义一个空数组，以便以后存储数值！
YY1(k  j),   表示取出每一个周期内的第一个解的最后一个值赋值到上边的这个空数组中！
YY1(:,end-51:end) 表示画积分200-60 后的200-60-51个数值的图像！

ruichard101

谢谢指点呵呵！我也试试

octopussheng

把你的东西弄过来了， 呵呵，不要见怪啊！

linyxff

你好，我也调试过这个程序做分岔图，这个duffing方程是一个外激方程，我现在要做一个参激方程的分岔图，即dx(2)=-r*y-x^3+F*x*cos(psi)，F后多一个x，分岔图就做不出来了，请问是什么原因啊？我是新手，可能理论就有些不对，希望大家帮帮忙！

学习者

楼主，用你贴出的程序计算得到的图，与你贴出的图完全不一样。附件我发不上去，不知如何将我的计算结果贴出来？盼望你的回复

无水1324

其本上差不多了，


附件你要点下面那个“发表新回复”之后就可以了

咕噜噜

原帖由 linyxff 于 2009-4-9 10:44 发表

                               
登录/注册后可看大图

你好，我也调试过这个程序做分岔图，这个duffing方程是一个外激方程，我现在要做一个参激方程的分岔图，即dx(2)=-r*y-x^3+F*x*cos(psi)，F后多一个x，分岔图就做不出来了，请问是什么原因啊？我是新手，可能理论就有 ...

不能吧，应该是问题不大啊，把你的程序贴上来瞅瞅

学习者

学习者

无水：上面是我拷贝楼主的程序，修改其中的range=[0.1:0.001:0.5]的计算图，为何与楼主贴出的图差距这么大？盼复

无水1324

确实差异好大哦，你Poincare截面取点没有问题吧？都是用1楼的程序嘛？

学习者

是的。我用该程序在不同的机子上算了4遍，结果相同