日志

.1. 关于连续系统Lyapunov指数的计算方法

已有 1124 次阅读2011-5-4 16:20

【总结】Lyapunov指数的计算方法

近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！

1. 关于连续系统Lyapunov指数的计算方法
连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法

2007-6-1 00:14

nedu-fly

助理工程师

精华 0
积分 23
帖子 38
水位 76
技术分 0

关于定义法求解的程序，和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 a  0;
      z 0  x-c];
dX(4:12) = Jaco*Y;

求解LE代码：
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = [1,1,1];
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:tstep

tstart + tstep*steps);
[T,Y] = ode45('Rossler_ly', tspan, y);
% 取积分得到的最后一个时刻的值
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)];
%正交化
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 = [a1,a2,a3];

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

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

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

2007-6-1 16:20

leslie0805

助理工程师

精华 0
积分 24
帖子 40
水位 80
技术分 0

（2）Jacobian方法

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

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

2007-6-2 06:46

wl21c

助理工程师

精华 0
积分 35
帖子 59
水位 118
技术分 0

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

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

（1）Nicolis方法

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

2007-6-2 15:08

zhouxiangnj

助理工程师

精华 0
积分 51
帖子 85
水位 170
技术分 0

（2）Wolf方法

2007-6-3 13:43

sdfdedff

助理工程师

精华 0
积分 78
帖子 130
水位 260
技术分 0

Wolf方法的Matlab程序如下：
function lambda_1=lyapunov_wolf(data,N,m,tau,P)
%  该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
%  m: 嵌入维数
%  tau:时间延迟
%  data:时间序列
%  N:时间序列长度
%

:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|I－J|>

的相点中搜寻
%  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)

M-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

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

2007-6-3 22:22

gushouli

助理工程师

精华 0
积分 11
帖子 18
水位 35
技术分 0

（3）小数据量方法

说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧，呵呵！

2007-6-4 13:05

huweeson

助理工程师

精华 0
积分 11
帖子 19
水位 37
技术分 0

上面两种方法，即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([0 40 0 2e7]), grid on  ％绘制年份－功率谱曲线
ylabel('Power')
xlabel('Period(Years/Cycle)')

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

（3）嵌入维数
目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的G－P算法计算出序列的关联维数d，然后利用嵌入维数m>=2d+1，选取合适的嵌入维数。
G―P算法程序如下：
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
% 计算关联维数的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,

,ln_C(m,

);
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

收藏分享邀请举报

wangjianbao11的个人空间 http://home.vibunion.com/?166599 [收藏] [复制] [分享] [RSS]

日志

.1. 关于连续系统Lyapunov指数的计算方法

全部作者的其他最新日志

评论 (0 个评论)

wangjianbao11