计算李雅普诺夫指数，循环出错，请指点

huazi071783

Y.rar (74.62 KB, 下载次数: 18)

2011-7-5 21:23 上传

点击文件名下载附件
数据Y

我要计算杜芬混沌吸引子的李雅普诺夫指数，循环老是出错，请指点，下面是我的程序
a=2;R=10;h=100;delt_t=0.01
load('Y.mat')                           %Y是5001*2的杜芬混沌时间序列，plot(Y(:,1),Y(:,1))就得到相空间曲线
lambda=lyapunov(delt_t,Y,h,a,R)
%下面是计算李氏指数的函数

function lambda=lyapunov(delt_t,Y,h,a,R)
M=length(Y);                                 %Y的第1列是横坐标，第2列是纵坐标
date_N1=a/delt_t;                          %限制不是同一条轨道
date_N2=R/delt_t;                          %限制不是离的很远的轨道，减少计算量
for i=1:M-h
    for j=1:M-h
        num=0;
        if (abs(j-i)>date_N1)&(abs(j-i)<date_N2)      %寻找相空间中每个点的最近距离点，并记下该点下标
            num=num+1;
            index(num)=j;
            d_s(num)=((Y(j,1)-Y(i,1))^2+(Y(j,2)-Y(i,2))^2)^0.5;   % 求最短距离   
        else
            continue
        end
    end
        [d_s,ind]=sort(d_s);                           %求最近的临近点
     for k=1:h
            turn=1;
            jj=index(ind(turn));
           
            d1=((Y(i,:)-Y(i+turn,:))*(Y(i,:)-Y(i+turn,:))')^0.5;
            d2=((Y(jj+turn,:)-Y(jj,:))*(Y(jj+turn,:)-Y(jj,:))')^0.5;
            cosine=(Y(i+turn,:)-Y(i,:))*(Y(jj+turn,:)-Y(jj,:))'/(d1*d2);         %求夹角
            ind1=0;
            if acos(cosine)<45*pi/180      %限制夹角小于45°
                   ind1=ind1+1;
                    D1(ind1)=((Y(i,:)-Y(jj,:))*(Y(i,:)-Y(jj,:))')^0.5;
                    D2(ind1)=((Y(i+k,:)-Y(jj+k,:))*(Y(i+k,:)-Y(jj+k,:))')^0.5;
                    lambda1(ind1)=log(D2(ind1)/D1(ind1))/(h*delt_t);
            else
            turn=turn+1;   
                continue   
            end
      %lambda(i)=max(lambda1);  %如果把这一行放这里，在循环里面就能计算，但不是我要的结果，我想循环后求最大值
     end
     lambda(i)=max(lambda1);                               %这句报错，
end

下面是报错，前面不是有lambda1吗？为什说没有定义呢？请指教
??? Undefined function or variable "lambda1".
Error in ==> lyapunov at 36
        lambda(i)=max(lambda1);

huazi071783

下面是得到的杜芬混沌吸引子，就是计算它的李雅普诺夫指数

huazi071783

怎么没人帮我看看啊，急啊！就是一个循环的问题，我水平有限不知道原因！

huazi071783

怎么没人帮我看看啊，急啊！就是一个循环的问题，我水平有限不知道原因！先

ChaChing

help下continue!
若acos(cosine)<45*pi/180不成立时, 是不会执行循环里面内其他指令的!
若那句在循环外头, 就会执行, 这时没有lambda1

huazi071783

回复 5 # ChaChing 的帖子

解决了，就是你说的问题，谢谢

ChaChing

回复 6 # huazi071783 的帖子

专业不同, 帮助有限!
仅建议若解决问题, 与大伙分享经验!

轩小轩

楼主怎么解决的这个问题，把lambda1(ind1)=log(D2(ind1)/D1(ind1))/(h*delt_t);这句放在什么地方啊？