最大Lyapunov指数计算问题请教！

wannete

%计算Duffing振子的lyapunov指数
clear all;
clc;
global X V
yinit=[1,1];
orthmatrix=[1 0;
            0 1];
y=zeros(6,1);
%初始化输入
IC(1:2)=yinit;
IC(3:6)=orthmatrix;
InitiaTime=0;  %时间初始值
TimeStep=1e-2;  %时间步长
wholetimes=100000;  %总的循环次数
UpdateStepNum=10;%每次演化的步数
iteratetimes =wholetimes/UpdateStepNum; %演化的次数
DiscardItr=200;
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitiaTime;
T1=InitiaTime;
T2=T1+Iteration;
Tspan=[T1:TimeStep:T2];
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
mod=zeros(1,2);
d=2;
Sum=zeros(1,d);
n=0;
k=0;
xData=[];
yData=[];

for i=1:iteratetimes
    n=n+1
    [t,X]=ode45('twofun',Tspan,IC,options);
    %取积分得到的最后一个时刻的值
    [rX,cX]=size(X);
    L=cX-d*d;
   
    for i=1:d
        m1=L+1+(i-1)*d;
        m2=m1+d-1;
        A(:,i)=(X(rX,m1:m2))';
    end
   y0=TwoGS(A);
   
    %取三个向量的模
    mod(1)=sqrt(y0(:,1)'*y0(:,1));
    mod(2)=sqrt(y0(:,2)'*y0(:,2));
    y0(:,1)=y0(:,1)/mod(1);
    y0(:,2)=y0(:,2)/mod(2);
    if T2>DiscardTime
        Q0=y0;
    else
        Q0=eye(d);
    end
   
    if (T2>DiscardTime)
        k=k+1;
        T=k*Iteration;
        TT=n*Iteration+InitiaTime;
        Sum=Sum+log(abs(mod));
        lambda=Sum/T;
        
        Lambda=fliplr(sort(lambda));
        [rxD,cxD]=size(xData);
         [ryD,cyD]=size(yData);
         
         xData=[xData;TT];
         yData=[yData;lambda];
    end
    ic=X(rX,1:L);
    T1=T1+Iteration;
    T2=T2+Iteration;
    Tspan=[T1:TimeStep:T2];
   
    IC=[ic(:);y0(:)];
end
plot(xData,yData(:,1),xData,yData(:,2))


function dX=twofun(t,X)
k=0.5;
B=2.804;
wd=500*2*pi;
x=X(1);
y=X(2);
Y=[X(3) X(5);
   X(4) X(6)];
dX=zeros(6,1);
dX(1)=wd*y;
dX(2)=wd*(-k*y-x^3+B*cos(wd*t));
J=[ 0,  wd*1;
wd*(-3*x^2),   wd*(-k)];
dX(3:6)=J*Y;





function A=TwoGS(V)  %V为3*3的向量
global a1 a2
v1=V(:,1);
v2=V(:,2);
a1=zeros(2,1);
a2=zeros(2,1);
a1=v1;
a2=v2-((a1'*v2)/(a1'*a1))*a1;
A=[a1,a2];
大家有谁懂的，分析一下这个程序 啊，也是出自本论坛，转载一下！自己用了一下，不是很明白！
需要计算一段时间内的最大指数，最后只得两个值！

[ 本帖最后由 wannete 于 2009-7-23 19:44 编辑 ]

octopussheng

先弄清一个概念，所需要求解的LE，就是计算时间终了时刻的LE，根据你这个程序，最后得到的两个值就是所要的结果。

这个程序就是GSR（GS正交化）方法求解的程序。相关理论课参考Wolf的那篇文章。

wannete

如果要改成计算一段时间的指数，比如就是Iteration时间内的最大指数，怎么改改啊，这样计算太慢了！

octopussheng

这个程序就是计算一段时间的LE的。

如果只要最大LE的话，可以根据最大LE的定义，编个简单程序即可求解。

xw263

图像上得到是两条曲线,也就是两个LE指数随积分时间变化的曲线.但并不意味着只有两个LE值.
实际上程序运行得到的LE指数是一个k*3的矩阵.而矩阵最后一行的三个数就是你所要的三个最终的LE指数.

程序运行完后,只要在命令窗口输入:yData 再按回车就可以得到了.如果想得到最大LE指数.只要输入
Lambda(1),即可.
程序时间应该并不长吧,大概1分半钟.

偶认为 plot(xData,yData(:,1),xData,yData(:,2))应改为:plot(xData,yData(:,1),xData,yData(:,2),xData,yData(:,3))
这样才完整一些.画出的是三条曲线.

xw263

如果想缩短运行时间,恐怕只能减少循环次数了.也就是减小wholetimes的值.

但是注意也不能太小的,因为总要除去程序所产生暂态数据点.比如这个程序中就要除去前面的200次循环.

另外,octopussheng总结的一篇贴子中,改进了这个程序,去掉了一些不必要的语句和循环,可能运行时间会短些.
你搜下,试试吧.

水平有限,还望多多指教.

octopussheng

说的还是很有道理的，呵呵。希望楼主能有所收获。

wannete

谢谢大家的指教！非常感谢！