Lyapunov指数的计算方法

1、【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有79本书！下面以吕金虎混沌时间序列分析及其应用、马军海复杂非线性系统的重构技术为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法    连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。（1）定义法定义法求解Lyapuno

2、v指数.JPG 关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例Rossler系统微分方程定义程序function dX = Rossler_ly(t,X)%  Rossler吸引子，用来计算Lyapunov指数%        a=0.15,b=0.20,c=10.0%        dx/dt = -y-z,%        dy/dt = x

3、+ay,%        dz/dt = b+z(x-c),a = 0.15;b = 0.20;c = 10.0;x=X(1); y=X(2); z=X(3);% Y的三个列向量为相互正交的单位向量Y = X(4), X(7), X(10);    X(5), X(8), X(11);    X(6), X(9), X(12);% 输出向量的初始化，必不可少dX = zeros(12,1);% Rossler吸引子dX(1) = -y-z;dX(2) = x+a*y;dX(3) = b+z*(x

4、-c);% Rossler吸引子的Jacobi矩阵Jaco = 0 -1 -1;         1 a   0;         z 0  x-c;dX(4:12) = Jaco*Y;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = 1,1,1;orthmatrix = 1 0 0;            

5、0; 0 1 0;              0 0 1;a = 0.15;b = 0.20;c = 10.0;y = zeros(12,1);% 初始化输入y(1:3) = yinit;y(4:12) = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数 iteratetimes = wholetimes/steps; % 演化的次数mod

6、 = zeros(3,1);lp = zeros(3,1);% 初始化三个Lyapunov指数Lyapunov1 = zeros(iteratetimes,1);Lyapunov2 = zeros(iteratetimes,1);Lyapunov3 = zeros(iteratetimes,1);for i=1:iteratetimes    tspan = tstart:tstep:(tstart + tstep*steps);       T,Y = ode45('Rossler_ly', tspan, y);&#

7、160;   % 取积分得到的最后一个时刻的值    y = Y(size(Y,1),:);    % 重新定义起始时刻    tstart = tstart + tstep*steps;    y0 = y(4) y(7) y(10);          y(5) y(8) y(11);          y(6) y(9) y(12);    %正交化 

8、   y0 = ThreeGS(y0);     % 取三个向量的模    mod(1) = sqrt(y0(:,1)'*y0(:,1);    mod(2) = sqrt(y0(:,2)'*y0(:,2);    mod(3) = sqrt(y0(:,3)'*y0(:,3);    y0(:,1) = y0(:,1)/mod(1);    y0(:,2) = y0(:,2)/mod(2);    y0(:,3) = y0

9、(:,3)/mod(3);    lp = lp+log(abs(mod);    %三个Lyapunov指数    Lyapunov1(i) = lp(1)/(tstart);    Lyapunov2(i) = lp(2)/(tstart);    Lyapunov3(i) = lp(3)/(tstart);        y(4:12) = y0'end% 作Lyapunov指数谱图i = 1:iteratetimes;p

10、lot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGS(V)  % V 为3*3向量v1 = V(:,1);v2 = V(:,2);v3 = V(:,3);a1 = zeros(3,1);a2 = zeros(3,1);a3 = zeros(3,1);a1 = v1;a2 = v2-(a1'*v2)/(a1'*a1)*a1;a3 = v3-(a1'*v3)/(a1'*a1)*a1-(a2'*v3)/(a2'*a2

11、)*a2;A = a1,a2,a3;计算得到的Rossler系统的LE为  0.063231  0.092635  -9.8924Wolf文章中计算得到的Rossler系统的LE为0.09   0   -9.77需要注意的是定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。正交化程序可以根据上面的扩展到N*N向量，这里就不加以说明了，对matlab用户来说应该还是比较简单的！（2）Jacobian方法 通过资料检索，发现论坛中用的较多的LET工具箱的算法原理就是Jacobi

12、an方法。基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维。经过计算也证明了这种方法精度较高，对目前常见的混沌系统，如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高，而且程序编写有一定的规范，个人很推荐使用。（虽然我自己要做的系统并不适用 ）LET工具箱可以在网络上找到，这里就不列出了！关于LET工具箱如果有问题，欢迎加入本帖讨论！Jacobian法求解Lyapunov指数.JPG 对离散动力系统，或者说是非线性时间序列，往往不需要计算出所有的Lyapun

13、ov指数，通常只需计算出其最大的Lyapunov指数即可。“1983年，格里波基证明了只要最大Lyapunov指数大于零，就可以肯定混沌的存在”。目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种：（1）由定义法延伸的Nicolis方法（2）Jacobian方法（3）Wolf方法（4）P范数方法（5）小数据量方法其中以Wolf方法和小数据量方法应用最为广泛，也最为普遍。下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍。（1）Nicolis方法 这种方法和连续系统的定义方法类似，而且目前应用很有限制，因此只对其理论进行介绍，编程应用方面就省略了 Nicolis方法

14、求最大LE.JPG（2）Wolf方法 Wolf方法求最大LE.JPGWolf方法的Matlab程序如下：function lambda_1=lyapunov_wolf(data,N,m,tau,P)%  该函数用来计算时间序列的最大Lyapunov 指数-Wolf 方法%  m: 嵌入维数 ！ 一般选大于等于10%  tau:时间延迟 ！一般选与周期相当，如我选2000%  data:时间序列 ！可以选1000；%  N:时间序列长度 满足公式：M=N-(m-1)*tau=24000-(10-1

15、)*1000=5000%  P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|IJ|>P的相点中搜寻 ! P=周期=600%  lambda_1: 返回最大lyapunov指数值min_point=1   %&&要求最少搜索到的点数MAX_CISHU=5 ;  %&&最大增加搜索范围次数%FLYINGHAWK%   求最大、最小和平均相点距离    max_d = 0;   &#

16、160;                                       %最大相点距离    min_d = 1.0e+100;                     

17、60;         %最小相点距离    avg_dd = 0;    Y=reconstitution(data,N,m,tau);           %相空间重构 可将此段程序加到 整个程序中，在时间循环内，可以保存时间序列的地方。见完整程序。    M=N-(m-1)*tau;               

18、                     %重构相空间中相点的个数    for i = 1 : (M-1)        for j = i+1 : M            d = 0;            for

19、 k = 1 : m                d = d + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);            end            d = sqrt(d);            if max_d < d

20、               max_d = d;            end            if min_d > d               min_d = d;      &#

21、160;     end            avg_dd = avg_dd + d;        end    end    avg_d = 2*avg_dd/(M*(M-1);                %平均相点距离        dl

22、t_eps = (avg_d - min_d) * 0.02 ;         %若在min_epsmax_eps中找不到演化相点时，对max_eps的放宽幅度    min_eps = min_d + dlt_eps / 2 ;            %演化相点与当前相点距离的最小限    max_eps = min_d + 2 * dlt_eps       &

23、#160;     %&&演化相点与当前相点距离的最大限    %     从P+1M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK    DK = 1.0e+100;                             %第i个相点

24、到其最近距离点的距离    Loc_DK = 2;                                 %第i个相点对应的最近距离点的下标    for i = (P+1):(M-1)                

25、       %限制短暂分离，从点P+1开始搜索        d = 0;        for k = 1 : m            d = d + (Y(k,i)-Y(k,1)*(Y(k,i)-Y(k,1);        end      

26、;  d = sqrt(d);        if (d < DK) & (d > min_eps)            DK = d;           Loc_DK = i;        end    end%  以下计算各相点对应的李氏数

27、保存到lmd()数组中%  i 为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK为相点i-1对应最短 距离的相点位置，DK为其对应的最短距离% Loc_DK+1为Loc_DK的演化点，DK1为i点到Loc_DK+1点的距离，称为演化距离%  前i个log2（DK1/DK）的累计和用于求i点的lambda值    sum_lmd = 0 ;                       

28、0;      % 存放前i个log2（DK1/DK）的累计和    for i = 2 : (M-1)                           % 计算演化距离              DK1 = 0;     &

29、#160;  for k = 1 : m            DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1)*(Y(k,i)-Y(k,Loc_DK+1);        end        DK1 = sqrt(DK1);        old_Loc_DK = Loc_DK ;  &

30、#160;               % 保存原最近位置相点        old_DK=DK;%     计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数        if (DK1 = 0)&( DK = 0)          

31、0;sum_lmd = sum_lmd + log(DK1/DK) /log(2);        end        lmd(i-1) = sum_lmd/(i-1); 此处可以保存不同相点i对应的李氏指数，见完整程序。% 以下寻找i点的最短距离:要求距离在指定距离范围内尽量短,与DK1的角度最小        point_num = 0; % &&在指定距离范围内找到的候选相点的个数&

32、#160;       cos_sita = 0  %&&夹角余弦的比较初值 要求一定是锐角        zjfwcs=0  %&&增加范围次数         while (point_num = 0)           % * 搜索相点      &

33、#160;     for j = 1 : (M-1)                if abs(j-i) <=(P-1)      %&&候选点距当前点太近，跳过！                   continue;      

34、              end                                %*计算候选点与当前点的距离                dnew = 0;   

35、0;            for k = 1 : m                   dnew = dnew + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);                end          

36、     dnew = sqrt(dnew);                                if (dnew < min_eps)|( dnew > max_eps )   %&&不在距离范围，跳过！            

37、      continue;                             end                                  

38、;             %*计算夹角余弦及比较                DOT = 0;                for k = 1 : m                  

39、0; DOT = DOT+(Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,old_Loc_DK+1);                end                CTH = DOT/(dnew*DK1);                    

40、            if acos(CTH) > (3.14151926/4)      %&&不是小于45度的角，跳过！                  continue;                end   &#

41、160;                            if CTH > cos_sita    %&&新夹角小于过去已找到的相点的夹角，保留                    cos_sita = CTH;   &

42、#160;                Loc_DK = j;                    DK = dnew;                end          

43、;      point_num = point_num +1;                            end                            if p

44、oint_num <= min_point               max_eps = max_eps + dlt_eps;               zjfwcs =zjfwcs +1;               if zjfwcs > MAX_CISHU  

45、0; %&&超过最大放宽次数，改找最近的点                   DK = 1.0e+100;                   for ii = 1 : (M-1)                  

46、60;   if abs(i-ii) <= (P-1)       %&&候选点距当前点太近，跳过！                       continue;                       

47、;    end                      d = 0;                      for k = 1 : m                  &

48、#160;       d = d + (Y(k,i)-Y(k,ii)*(Y(k,i)-Y(k,ii);                      end                      d = sqrt(d);     

49、0;                        if (d < DK) & (d > min_eps)                          DK = d;          

50、;               Loc_DK = ii;                      end                   end           

51、0;       break;                end               point_num = 0 ;      %&&扩大距离范围后重新搜索               cos

52、_sita = 0;            end        end   end%取平均得到最大李雅普诺夫指数（此处只有一个值，若为正说明体系是混沌的，若为负则说明体系是周期性的确定性运动）lambda_1=sum(lmd)/length(lmd);程序中用到的reconstitution函数如下： 此段程序可直接放在时间循环内部，即可以保存时间序列的地方。见完整程序范例。function X=reconstitution(data

53、,N,m,tau)%该函数用来重构相空间% m 为嵌入空间维数% tau 为时间延迟% data 为输入时间序列% N 为时间序列长度% X 为输出,是m*n维矩阵M=N-(m-1)*tau; %相空间中点的个数for j=1:M            %相空间重构    for i=1:m        X(i,j)=data(i-1)*tau+j);    endend这里声明一下，这些程序并非我自己编写的，均是

54、转载，其使用我已经验证过，绝对可以运行！（3）小数据量方法 说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧，呵呵！ 小数据量方法求最大Lyapunov指数.JPG 上面两种方法，即Wolf方法和小数据量方法，在计算LE之前，都要求对时间序列进行重构相空间，重构相空间的优良对于最大LE的计算精度影响非常大！因此重构相空间的几个参数的确定就非常重要。（1）时间延迟主要推荐两种方法自相关函数法、CC方法自相关函数法对一个混沌时间序列，可以先写出其自相关函数，然后作出自相关函数关于时间t的函数图像。根据数值试验结果，当自相关函数下降到初始值的11/e时，所得的时间t即为重构相空间的时间延迟。

55、CC方法可以同时计算出时间延迟和时间窗口，个人推荐使用这种方法！（2）平均周期平均周期的计算可以采用FFT方法。在matlab帮助中有一个太阳黑子的例子，现摘录如下：load sunspot.dat             %装载数据文件year = sunspot(:,1);       %读取年份信息wolfer = sunspot(:,2);    %读取黑子活动数据plot(year,wolfer)      

56、;       %绘制原始数据图title('Sunspot Data')Y = fft(wolfer);               %快速FFT变换N = length(Y);                %FFT变换后数据长度Y(1) = ;           &

57、#160;            %去掉Y的第一个数据，它是所有数据的和power = abs(Y(1:N/2).2;   %求功率谱nyquist = 1/2;               freq = (1:N/2)/(N/2)*nyquist; %求频率plot(freq,power), grid on         %绘制功率谱

58、图xlabel('cycles/year')title('Periodogram')period = 1./freq;                      %年份（周期）plot(period,power), axis(0 40 0 2e7), grid on   绘制年份功率谱曲线ylabel('Power')xlabel('Period(Years/Cycle)'

59、)mp,index = max(power);       %求最高谱线所对应的年份下标period(index)                            %由下标求出平均周期（3）嵌入维数目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的GP算法计算出序列的关联维数d，然后利用嵌入维数m>=2d+1，选取合适的嵌入维数。GP算

60、法程序如下：function ln_r,ln_C=G_P(data,N,tau,min_m,max_m,ss)% the function is used to calculate correlation dimention with G-P algorithm%    计算关联维数的GP算法% data:the time series                        时间序列% N: the leng

61、th of the time series            时间序列长度% tau: the time delay                         时间延迟% min_m:the least embedded dimention m       最小的嵌入维数% max_m:the largest

62、 embedded dimention m     最大的嵌入维数% ss:the stepsize of r                        r的步长%skyhawkfor m=min_m:max_m    Y=reconstitution(data,N,m,tau);%reconstitute state space    M=N-(m-1

63、)*tau;%the number of points in state space    for i=1:M-1        for j=i+1:M            d(i,j)=max(abs(Y(:,i)-Y(:,j);%calculate the distance of each two              

64、0;    end                     %points in state space  计算状态空间中每两点之间的距离    end    max_d=max(max(d);%the max distance of all points   得到所有点之间的最大距离    d(1,1)=max_d; 

65、;   min_d=min(min(d);%the min distance of all points   得到所有点间的最短距离    delt=(max_d-min_d)/ss;%the stepsize of r             得到r的步长    for k=1:ss        r=min_d+k*delt;        C(k)=correlation