Duffing系统的分岔图的Matlab实现

 振动资讯 V:~
ghGi"B:w.T#xK

function dx=duffing(t,X)振动资讯^_|/Jz'W#e-Xo G)j

global F wd;振动资讯N.v'b eFk!J

r=0.168;振动资讯#Xq%soug.@.r6b

x=X(1);

y=X(2);

psi=X(3);振动资讯"OM8v:C7S]

dx=zeros(3,1);

dx(1)=y;振动资讯%{BLJORjZt
Z

dx(2)=-r*y+1/2*x*(1-x^2)+F*sin(psi);振动资讯:`1vc)|2c jm

dx(3)=wd;

clear;振动资讯th&\:K[ R

global F wd;

wd=1.0;

range=[0.10:0.0001:0.3];振动资讯DBa.b#q v)?

period=2*pi/wd; %

k=0;振动资讯ca:T"l)d'n~`

YY2=[];振动资讯L0duZ,?3}(JR

step=2*pi/100;  %步长。振动资讯J2tXJi O g}

for F=range振动资讯b5jq{/U1_7VM

    y0=[0 0 0];振动资讯wIfU:s'{

   F

    k=k+1;

    % discard the first 60 periodic data;

    %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值

    tspan=[0:step:60*period];振动资讯\&x/GhL3r7`

    [t,Y]=ode45(@duffing,tspan,y0);振动资讯_:]N"O8tx0L'\

    y0=Y(end,:);

    j=1;

    for i=60:200

        tspan=[i*period:step:(i+1)*period];振动资讯K/K LcD1uK

        [t,Y]=ode45(@duffing,tspan,y0);

        YY1(k,j)=Y(end,1);   % get the omega data from every period end

        j=j+1;               %取出每一个周期内的第一个解的最后一个值。振动资讯mbc1W8O"O}r

        y0=Y(end,:);振动资讯'|F`n
U9rZ

    end振动资讯K:WY4W9l r}X,Q(^]

end振动资讯[m)vy C3N6V3}ryk

bifdata=YY1(:,end-51:end);振动资讯[b{]#e'V3o2]6M

plot(range,bifdata,'k.','markersize',1);

Duffing Bifurcation