liapunov.m

计算Lyapunov指数,不知道正确与否

%这个程序计算Liapunov指数
%计算方法，取两个相邻很近的点作为初始点，然后每隔一个周期计算这两点之间的距离，不过，每次的初始点必须保证两个初始点之间的距离相等，而且新点与两个旧点
%必须在同一条直线上，具体参考混沌动力学
canshu;
em=1;
W=0.325;
%for W=1:-0.005:0.005
    calnum=1;
    period=2*pi/W;
    step=period/2000;
    e=1;
    k=50;
    tspan=0*period:step:k*period;          %先扔掉k个周期的数据
    x1=[0,0,0,0,0,0,0,0]';
    x1=shuzhidy(tspan,x1);
    x2=x1+1e-5*ones(8,1);
    d0=sqrt((x1-x2)'*(x1-x2));
    tn=1;
    step=step/100; 
    k=1;
    tspan=(k-1)*100*step:step:k*100*step;
    x11=shuzhidy(tspan,x1);
    x22=shuzhidy(tspan,x2);
    d1=sqrt((x11-x22)'*(x11-x22));
    e0=log(d1/d0);
    namuda0=e0/(tn*step);
    while(e>1e-5)
        k=k+1;
        tn=tn+1;
        tspan=(k-1)*100*step:step:k*100*step;
        %先更新计算点，使其到x1的距离仍然是d0
        x1=x11;
        x2=x1+(d0/sqrt(8))*(x22-x11)./abs(x22-x11);
        d0=sqrt((x1-x2)'*(x1-x2));
        x11=shuzhidy(tspan,x1);
        x22=shuzhidy(tspan,x2);
        d1=sqrt((x11-x22)'*(x11-x22));
        e0=e0+log(d1/d0);
        namuda1=e0/(tn*step);
        display([W,tn,namuda1,log(d1/d0)]);
        e=abs(namuda1-namuda0);
        namuda0=namuda1;
        calnum=calnum+1;
        if(calnum>1000)
            break;
        end
    end
    ew(em)=namuda1;
    save ew ew;
    em=em+1;
%end
%W=1:-0.0025:0.005;
%plot(W,ew);