Lyapunov指数的计算方法

【总结】Lyapunov指数的计算方法非线性理论

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

1.关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方

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

（1）定义法

—对H堆连续动力系統z = 在—OBJ孙味“为中心.|拆（心0）||为丰笹啟存

«堆的球面*施著时间的演化，在t时討该球而0P变形为M继的椭球厨・设滾椭域面的第/ 个坐标轴方向的半轴长対卩兀|，则诙系统第i个指数対*

此即连续系统Lyapunov揩较飽定冥・

弼计尊时・取|处（心0）［为岀W为常数），以孔为球心・欧几里篇范敢为山的正衮矢量集仙测，…叮为初始球.由非线性徴分方崔“尸㈤可以分别计算出点細血创、引他、r引址经过时间t后淺化的轨迹・役其终了点分别为珊、砒、f 仙则令石f 陶一心■处严=甩-和，r 亦耳国=略一報#则可得新的矢重棄叶禺巴…后畀｝・

由于各牛妥臺在演化过程中舌焙向着是大的UapurOT IS数方何靠掘，因此需要通过Schimdt IE交化不断地讨新矢量逬行置换.SP Wolf to文章中提出的GSR^法.表述如下：

播着以他为球心，疤数対（I的正奁矢臺料创巴叫叫…伽严;为祈球继续进行演出设演化

至N步时，得到矢董慕冈㈤出巴…耳僅牛且N足够大,这可以得到Lyapunov 扌鐵的近似计算公式三

实际计算时，取为1・

定义法求解Lyapunov 指数JPG

关于定义法求解的程序，和matlab板块的连续系统LE求解程序”差不多。以

Rossie啄统为例

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的三个列向量为相互正交的单位向量

丫 = [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);

Ip = lp+log(abs(mod));

%三个Lyapunov指数

Lyapu nov1(i) = lp(1)/(tstart);

Lyapu nov2(i) = lp(2)/(tstart);

Lyapu nov3(i) = lp(3)/(tstart); y(4:12) = y0';

end

%作Lyapunov指数谱图

i = 1:iteratetimes;

plot(i,Lyap uno v1,i,Lyap uno v2,i,Lyap unov3)

程序中用到的ThreeGS程序如下：

%G-S正交化

function A = ThreeGS(V) % V 为3*3 向量

v1 = V(:,1);