关于连续系统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 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:tsteptstart + 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.24

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:时间序列长度

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

% 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