%這個(gè)程序計(jì)算Liapunov指數(shù)
%計(jì)算方法，取兩個(gè)相鄰很近的點(diǎn)作為初始點(diǎn)，然后每隔一個(gè)周期計(jì)算這兩點(diǎn)之間的距離，不過，每次的初始點(diǎn)必須保證兩個(gè)初始點(diǎn)之間的距離相等，而且新點(diǎn)與兩個(gè)舊點(diǎn)
%必須在同一條直線上，具體參考混沌動(dòng)力學(xué)
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個(gè)周期的數(shù)據(jù)
    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;
        %先更新計(jì)算點(diǎn)，使其到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);