mumianhua 发表于 2014-5-27 21:12

Lyapunov指数程序————大家讨论一哈

本帖最后由 牛小贱 于 2014-5-28 12:41 编辑


大家好,希望大家一起探讨一下Lyapunov指数!!!
我们都知道,Lyapunov指数是一个衡量系统动力学特性的重要指标,它表征系统在相空间中相邻轨道间收敛
或者是发散的平均指数率,对于系统是否存在动力学混沌,可以从最大的Lyapunov指数是否大于0来直观的判断。
在我研究的非线性转子动力学,遇到的问题就是Lyapunov指数不会求解,希望哪位大神能够给予Lyapunov指数相关程序

1. 关于连续系统Lyapunov指数的计算方法
    连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列(即离散系统)的LE求解方法来计算得到。
    关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容。
(1)定义法
关于定义法求解的程序,和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+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-c);
% Rossler吸引子的Jacobi矩阵
Jaco = [0 -1 -1;
      1 a0;
      z 0x-c];
dX(4:12) = Jaco*Y;

求解LE代码:
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = ;
orthmatrix = [1 0 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 = 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:tstepundefinedtstart + tstep*steps);   
    = ode45('Rossler_ly', tspan, y);
    % 取积分得到的最后一个时刻的值
    y = Y(size(Y,1),undefined;
    % 重新定义起始时刻
    tstart = tstart + tstep*steps;
    y0 = [y(4) y(7) y(10);
          y(5) y(8) y(11);
          y(6) y(9) y(12)];
    %正交化
    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(:,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;
plot(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))*a2;
A = ;

计算得到的Rossler系统的LE为――――0.0632310.092635-9.8924
Wolf文章中计算得到的Rossler系统的LE为――――0.09   0   -9.77


需要注意的是――定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象。

正交化程序可以根据上面的扩展到N*N向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的!

(2)Jacobian方法

通过资料检索,发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法。
基本原理就是首先求解出连续系统微分方程的近似解,然后对系统的Jacobian矩阵进行QR分解,计算Jacobian矩阵特征值的乘积,最后计算出LE和分数维。
经过计算也证明了这种方法精度较高,对目前常见的混沌系统,如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高,而且程序编写有一定的规范,个人很推荐使用。(虽然我自己要做的系统并不适用)

LET工具箱可以在网络上找到,这里就不列出了!关于LET工具箱如果有问题,欢迎加入本帖讨论!

对离散动力系统,或者说是非线性时间序列,往往不需要计算出所有的Lyapunov指数,通常只需计算出其最大的Lyapunov指数即可。“1983年,格里波基证明了只要最大Lyapunov指数大于零,就可以肯定混沌的存在”。


目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种:
(1)由定义法延伸的Nicolis方法
(2)Jacobian方法
(3)Wolf方法
(4)P-范数方法
(5)小数据量方法
其中以Wolf方法和小数据量方法应用最为广泛,也最为普遍。
下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍。


(1)Nicolis方法

这种方法和连续系统的定义方法类似,而且目前应用很有限制,因此只对其理论进行介绍,编程应用方面就省略了

(2)Wolf方法
Wolf方法的Matlab程序如下:
function lambda_1=lyapunov_wolf(data,N,m,tau,P)
%该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
%m: 嵌入维数
%tau:时间延迟
%data:时间序列
%N:时间序列长度
%undefined:时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为I,则演化相点只能在|I-J|>undefined的相点中搜寻
%lambda_1:返回最大lyapunov指数值
min_point=1; %&&要求最少搜索到的点数
MAX_CISHU=5 ;%&&最大增加搜索范围次数
%FLYINGHAWK
%   求最大、最小和平均相点距离
    max_d = 0;                                       %最大相点距离
    min_d = 1.0e+100;                                  %最小相点距离
    avg_dd = 0;
    Y=reconstitution(data,N,m,tau);                  %相空间重构
    M=N-(m-1)*tau;                                     %重构相空间中相点的个数
    for i = 1 : (M-1)
      for j = i+1 : M
            d = 0;
            for 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
               max_d = d;
            end
            if min_d > d
               min_d = d;
            end
            avg_dd = avg_dd + d;
      end
    end
    avg_d = 2*avg_dd/(M*(M-1));                %平均相点距离
   
    dlt_eps = (avg_d - min_d) * 0.02 ;         %若在min_eps~max_eps中找不到演化相点时,对max_eps的放宽幅度
    min_eps = min_d + dlt_eps / 2 ;            %演化相点与当前相点距离的最小限
    max_eps = min_d + 2 * dlt_eps;         %&&演化相点与当前相点距离的最大限
   
%   从P+1~M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK
    DK = 1.0e+100;                           %第i个相点到其最近距离点的距离
    Loc_DK = 2;                              %第i个相点对应的最近距离点的下标
    for i = (P+1)undefinedM-1)                        %限制短暂分离,从点P+1开始搜索
      d = 0;
      for k = 1 : m
            d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));
      end
      d = sqrt(d);
      if (d < DK) & (d > min_eps)
         DK = d;
         Loc_DK = i;
      end
    end
%   以下计算各相点对应的李氏数保存到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 ;                              % 存放前i个log2(DK1/DK)的累计和
    for i = 2 : (M-1)                        % 计算演化距离      
      DK1 = 0;
      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 ;                  % 保存原最近位置相点
      old_DK=DK;
%   计算前i个log2(DK1/DK)的累计和以及保存i点的李氏指数
      if (DK1 ~= 0)&( DK ~= 0)
         sum_lmd = sum_lmd + log(DK1/DK) /log(2);
      end
      lmd(i-1) = sum_lmd/(i-1);
%   以下寻找i点的最短距离:要求距离在指定距离范围内尽量短,与DK1的角度最小
      point_num = 0; % &&在指定距离范围内找到的候选相点的个数
      cos_sita = 0; %&&夹角余弦的比较初值 ――要求一定是锐角
      zjfwcs=0   ;%&&增加范围次数
         while (point_num == 0)
         % * 搜索相点
            for j = 1 : (M-1)
                if abs(j-i) <=(P-1)      %&&候选点距当前点太近,跳过!
                   continue;   
                end
               
                %*计算候选点与当前点的距离
                dnew = 0;
                for k = 1 : m
                   dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
                end
                dnew = sqrt(dnew);
               
                if (dnew < min_eps)|( dnew > max_eps )   %&&不在距离范围,跳过!
                  continue;            
                end
                              
                %*计算夹角余弦及比较
                DOT = 0;
                for k = 1 : m
                  DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));
                end
                CTH = DOT/(dnew*DK1);
               
                if acos(CTH) > (3.14151926/4)      %&&不是小于45度的角,跳过!
                  continue;
                end
               
                if CTH > cos_sita   %&&新夹角小于过去已找到的相点的夹角,保留
                  cos_sita = CTH;
                  Loc_DK = j;
                  DK = dnew;
                end
                point_num = point_num +1;
               
            end      
      
            if point_num <= min_point
               max_eps = max_eps + dlt_eps;
               zjfwcs =zjfwcs +1;
               if zjfwcs > MAX_CISHU    %&&超过最大放宽次数,改找最近的点
                   DK = 1.0e+100;
                   for ii = 1 : (M-1)
                      if abs(i-ii) <= (P-1)      %&&候选点距当前点太近,跳过!
                     continue;   
                      end
                      d = 0;
                      for k = 1 : m
                        d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));
                      end
                      d = sqrt(d);
      
                      if (d < DK) & (d > min_eps)
                         DK = d;
                         Loc_DK = ii;
                      end
                   end
                   break;
               end
               point_num = 0          ;   %&&扩大距离范围后重新搜索
               cos_sita = 0;
            end
      end
   end
%取平均得到最大李雅普诺夫指数
lambda_1=sum(lmd)/length(lmd);


程序中用到的reconstitution函数如下:
function X=reconstitution(data,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);
    end
end


这里声明一下,这些程序并非我自己编写的,均是转载,其使用我已经验证过,绝对可以运行!

(3)小数据量方法说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧,呵呵!


*** 上面两种方法,即Wolf方法和小数据量方法,在计算LE之前,都要求对时间序列进行重构相空间,重构相空间的优良对于最大LE的计算精度影响非常大!因此重构相空间的几个参数的确定就非常重要。

(1)时间延迟
主要推荐两种方法――自相关函数法、C-C方法
自相关函数法――对一个混沌时间序列,可以先写出其自相关函数,然后作出自相关函数关于时间t的函数图像。根据数值试验结果,当自相关函数下降到初始值的1-1/e时,所得的时间t即为重构相空间的时间延迟。
C-C方法――可以同时计算出时间延迟和时间窗口,个人推荐使用这种方法!

(2)平均周期
平均周期的计算可以采用FFT方法。在matlab帮助中有一个太阳黑子的例子,现摘录如下:
load sunspot.dat             %装载数据文件
year = sunspot(:,1);       %读取年份信息
wolfer = sunspot(:,2);    %读取黑子活动数据
plot(year,wolfer)             %绘制原始数据图
title('Sunspot Data')

Y = fft(wolfer);               %快速FFT变换

N = length(Y);                %FFT变换后数据长度
Y(1) = [];                     %去掉Y的第一个数据,它是所有数据的和
power = abs(Y(1:N/2)).^2;%求功率谱
nyquist = 1/2;               
freq = (1:N/2)/(N/2)*nyquist;%求频率
plot(freq,power), grid on      %绘制功率谱图
xlabel('cycles/year')
title('Periodogram')

period = 1./freq;                     %年份(周期)
plot(period,power), axis(), grid on%绘制年份-功率谱曲线
ylabel('Power')
xlabel('Period(Years/Cycle)')

= max(power);       %求最高谱线所对应的年份下标
period(index)                           %由下标求出平均周期

(3)嵌入维数
目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的G-P算法计算出序列的关联维数d,然后利用嵌入维数m>=2d+1,选取合适的嵌入维数。
G―P算法程序如下:
function =G_P(data,N,tau,min_m,max_m,ss)
% the function is used to calculate correlation dimention with G-P algorithm
%    计算关联维数的G-P算法
% data:the time series                     时间序列
% N: the length of the time series         时间序列长度
% tau: the time delay                        时间延迟
% min_m:the least embedded dimention m       最小的嵌入维数
% max_m:the largest embedded dimention m   最大的嵌入维数
% ss:the stepsize of r                     r的步长
%skyhawk
for m=min_m:max_m
    Y=reconstitution(data,N,m,tau);%reconstitute state space
    M=N-(m-1)*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         
      end                              %points in state space计算状态空间中每两点之间的距离
    end
    max_d=max(max(d));%the max distance of all points   得到所有点之间的最大距离
    d(1,1)=max_d;
    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_integral(Y,M,r);%calculate the correlation integral
      ln_C(m,k)=log(C(k));%lnC(r)
      ln_r(m,k)=log(r);%lnr
      fprintf('%d/%d/%d/%d\n',k,ss,m,max_m);
    end
    plot(ln_r(m,undefined,ln_C(m,undefined);
    hold on;
end
fid=fopen('lnr.txt','w');
fprintf(fid,'%6.2f %6.2f\n',ln_r);
fclose(fid);
fid = fopen('lnC.txt','w');
fprintf(fid,'%6.2f %6.2f\n',ln_C);
fclose(fid);

程序中的correlation_integral函数如下:
function C_I=correlation_integral(X,M,r)
%the function is used to calculate correlation integral
%C_I:the value of the correlation integral
%X:the reconstituted state space,M is a m*M matrix
%m:the embedding demention
%M:M is the number of embedded points in m-dimensional sapce
%r:the radius of the Heaviside function,sigma/2<r<2sigma
%calculate the sum of all the values of Heaviside
%skyhawk
sum_H=0;
for i=1:M
%   fprintf('%d/%d\n',i,M);
    for j=i+1:M
      d=norm((X(:,i)-X(:,j)),inf);%calculat the distances of each two points in matris M with sup-norm
      sita=heaviside(r,d);%calculate the value of the heaviside function
      sum_H=sum_H+sita;
    end
end
C_I=2*sum_H/(M*(M-1));%the value of correlation integral


以上的各种方法在实际应用的时候要根据具体情况来选择。

一般地,如果已知系统方程(当然系统不能太过复杂)时,则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些!

而如果系统方程比较复杂(如超维系统)、或者为一时间序列时,则推荐采样Wolf方法、小数据量方法。

Wolf方法的特点是时间序列无噪声,空间中小向量的演变高度非线性,而Jacobian方法则是噪声大,空间中小向量的演变接近线性。

小数据量方法的优点在于:(1)对小数据组的计算可靠;(2)计算量较小,比wolf方法快很多;(3)编程、操作较为容易。


而关于时间延迟、嵌入维数、平均周期的确定,还是推荐C-C方法和G-P算法,结果更为可靠一些!







truleeee 发表于 2014-5-28 00:13

把那橙色的背景去掉吧,能把人眼睛累死

mumianhua 发表于 2014-5-28 10:13

谁有好的李雅诺谱指数程序,可以交流交流哈相互学习

liguangzhigong 发表于 2014-5-29 13:17

您好,我请教一个问题。对于实际的系统,一般都含有噪声,而且噪声强度较大,那么是否就意味着,只能够采用小数据方法呢

VibInfo 发表于 2014-5-29 14:43

mumianhua 发表于 2014-5-28 10:13
谁有好的李雅诺谱指数程序,可以交流交流哈相互学习

参考http://forum.vibunion.com/thread-24201-1-1.html
别人已经汇总了很多了

VibInfo 发表于 2014-5-29 14:47

liguangzhigong 发表于 2014-5-29 13:17
您好,我请教一个问题。对于实际的系统,一般都含有噪声,而且噪声强度较大,那么是否就意味着,只能够采用 ...

由于在实际问题中,往往不知道动力系统的运动方程,再加上使用GSR方法计算相空间各维上的体积元的长期增长速率比较困难,使得原始定义计算Lyapunov指数几乎不可能;
Wolf方法适用于时间序列无噪声,切空间中小向量的演变高度非线性;
Jocobian方法适用于时间序列噪声大,切空间中小向量的演变接近线性;
而P-范数的选取和计算过于复杂,实际操作起来较为困难;
小数据量方法只适用于计算混沌时间序列的最大Lyapunov指数。

iwecabt 发表于 2014-7-9 10:19

非常感谢~~~~~~~~~~~











http://www.bb520hk.com/   http://ppspps.org/   http://www.tzingchu.com/

jlzhoushuang 发表于 2014-7-28 18:45

{:{03}:}

窃以为洒脱 发表于 2014-8-30 09:42

感谢分享

shenyongjun 发表于 2014-8-30 22:28

这个已经很多了吧?
谁有较好的时滞非线性或者分数阶微分方程的Lyapunov指数方法和程序?

fangdongzi 发表于 2015-10-5 12:27

shenyongjun 发表于 2014-8-30 22:28
这个已经很多了吧?
谁有较好的时滞非线性或者分数阶微分方程的Lyapunov指数方法和程序?

楼主解决了嘛?我也是分数阶的李雅普洛夫问题

suffer 发表于 2015-10-10 10:44

shenyongjun 发表于 2014-8-30 22:28
这个已经很多了吧?
谁有较好的时滞非线性或者分数阶微分方程的Lyapunov指数方法和程序?

Lyapunov指数计算一般是基于时间序列的,和方程具体是什么形式的好像关系不大吧
方程求解时间序列后就可以采用这些李指数分析程序了
页: [1]
查看完整版本: Lyapunov指数程序————大家讨论一哈