关于计算杜芬混沌吸引子的最大李雅普诺夫指数的问题

huazi071783

首先说明我是做土木结构损伤检测的，最近才研究用混沌振子来检测损伤，碰到就算混沌吸引子的李氏指数问题。刚开始以为挺简单，没想到这里水这么深。看了很多文献，参考一篇杨绍清的一篇《一种最大李雅普诺夫指数估计的稳健算法》，



我按他的方法编好了matlab程序，但是不知道对不对，水平有限。有几个问题不懂，请教
1，演化步长选取后，当i=1时第一次迭代计算完后，第二次迭代是取i=2还是取i=1+h？总的迭代次数是总长度M还是多少呢？
2，公式一中为什么是ip+jh，而不是ip+h呢？
下面是我的程序，不知道理解的正不正确

huazi071783

这是我计算得到的吸引子
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
        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
        [d_s,ind]=sort(d_s);                           %求最近的临近点
        for k=1:h
            d1=((Y(i,:)-Y(i+k,:))*(Y(i,:)-Y(i+k,:))')^0.5;
            jj=index(ind(1));
            d2=((Y(jj+k,:)-Y(jj,:))*(Y(jj+k,:)-Y(jj,:))')^0.5;
            cosine=(Y(i+k,:)-Y(i,:))*(Y(jj+k,:)-Y(jj,:))'/(d1*d2);         %求夹角
             ind1=0;
            if acos(cosine)<80*pi/180      %限制夹角小于80°
                   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);
                    lambda(i)=max(lambda1);
              break
              else
               continue   
            end
        end
    end
end

meiyongyuandeze

结构损伤也可以用混沌来判断啊，楼主可以不可以给点参考文献，很感兴趣呢！

huazi071783

回复 3 # meiyongyuandeze 的帖子

就是不知道行不行啊，现在在试着看

huazi071783

算是算出来了，但是结果明显不对，应该是绝大部分大于零的，变化不大的曲线才对。下面是我算的结果，请做非线性的高手来指点江山！！

huazi071783

哪位能帮帮忙，指点指点啊！！在线静候

huazi071783

我把程序从新调试好了，但是得出来的结果还是不对啊，曲线怎么这么乱七八糟的，不像上述文献那么平缓，而且小于零的很？求平均后得到的指数为0.8058。求指点！谢

huazi071783

调整后的程序，但结果不理想，望指点
function [lambda,LE]=lyapunov(delt_t,Y,h,a,R,bata)
% lambda=lyapunov(0.01,Y,100,2,20);
% load('Y.mat');
% delt_t=0.01;
% a=2;h=100;R=10;
% bata=40;
M=length(Y);                               %Y的第1列是横坐标，第2列是纵坐标
date_N1=a/delt_t;                          %限制不是同一条轨迹
date_N2=R/delt_t;                          %限制不是离的很远的轨迹，减少计算量
for i=1:M-h
        num=0;
    for j=1:M-h
        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);                           %求最近的临近点
    turn=1;
     for k=1:h
%             turn=1;
            jj=index(ind(turn));
           

            d1=((Y(i,:)-Y(jj,:))*(Y(i,:)-Y(jj,:))')^0.5;
            d2=((Y(i+turn,:)-Y(jj+turn,:))*(Y(i+turn,:)-Y(jj+turn,:))')^0.5;
            cosine=(Y(jj,:)-Y(i,:))*(Y(jj+turn,:)-Y(i+turn,:))'/(d1*d2);         %求夹角
            ind1=0;
            if acos(cosine)<bata*pi/180      %限制夹角小于bata
                    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))/(k*delt_t);
            else
            turn=turn+1;   
                continue   
            end
     end
     lambda(i)=max(lambda1);
%        lambda(i)=mean(lambda1);
end
LE=mean(lambda);

cqupenghao

回复 8 # huazi071783 的帖子

楼主好像是在用wolf法求解最大Lyapunov指数，如果这种方法不行，楼主可以尝试另一种方法：小数据量法。
另外我附上用LET工具箱求解的duffing equation的Lyapunov指数情况。

小影子

您好！请问您算出最大李雅普诺夫指数来了吗？能否讨论一下