duffing.m

来自「duffing分岔图」· M 代码 · 共 81 行

function dx=duffing(t,X)

global F wd;

r=0.168;

x=X(1);

y=X(2);

psi=X(3);

dx=zeros(3,1);

dx(1)=y;

dx(2)=-r*y+1/2*x*(1-x^2)+F*sin(psi);

dx(3)=wd;







clear;

global F wd;

wd=1.0;

range=[0.10:0.0001:0.3];

period=2*pi/wd; % 

k=0;

YY2=[];

step=2*pi/100;  %步长。

for F=range

    y0=[0 0 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);