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:iteratetimestspan = 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);v2 = V(:,2);v3 = V(:,3);al = zeros(3,1);a2 = zeros(3,l);a3 = zeros(3,1);al = 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.8924Wolf文章中计算得到的Rossler系统的LE为-------- 0.09 0 -9.77需要注意的是一一定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。

正交化程序可以根据上面的扩展到N*N向量，这里就不加以说明了，对matlab用户来说应该还是比较简单的！(2) Jacobian 方法通过资料检索，发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法。

基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian 矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维。

经过计算也证明了这种方法精度较高，对目前常见的混沌系统，如Lore nz、Henon、Duffing等的Lyapunov指数的计算精度都很高，而且程序编写有一定的规范，个人很推荐使用。

（虽然我自己要做的系统并不适用罢@）LET 工具箱可以在网络上找到，这里就不列出了！关于 LET 工具箱如果有问题， 欢迎加入本帖讨论！Jacobian 法求解 Lyapunov 指数 JPG 对离散动力系统，或者说是非线性时间序列，往往不需要计算出所有的Lyapunov 指数，通常只需计算出其最大的 Lyap u nov 指数即可。

“198年，格里波 基证明了只要最大Lyapunov 指数大于零，就可以肯定混沌的存在”。

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

下面对Nicolis 方法、Wolf 方法以及小数据量方法作 介绍。

（ 1）Nicolis 方法这种方法和连续系统的定义方法类似，而且目前应用很有限制，因此只对其理论 进行介绍，编程应用方面就省略了 Nicolis 方法求最大 LE.JPG2）Wolf 方法 Wolf 方法求最大 LE.JPGWolf 方法的 Matlab 程序如下：function lambda_1=lyapunov_wolf(data,N,m,tau,P) % 该函数用来计算时间序列的最大 Lyapunov 指数--Wolf 方法！ 一般选大于等于 10 ！ 一般选与周期相当，如我选 2000 ！可以选 1000；% N:时间序列长度 满足公式：M=N-（m-1）*tau=24000-（10-1）*1000=5000% P ：时间序列的平均周期，选择演化相点距当前点的位置差，即若当前相点为 I ，则 演化相点只能在|I — J|>P 的相点中搜寻 ！ P=周期=600% lambda_1: 返回最大lyapunov 指数值min_point=1 ; %&& 要求最少搜索到的点数MAX_CISHU=5 ; %&& 最大增加搜索范围次数%FLYINGHAWK% 求最大、最小和平均相点距离max_d = 0; min_d = 1.0e+100;% m: 嵌入维数% tau 时间延迟 % data 时间序列 %最大相点距离%最小相点距avg_dd = 0;Y=reconstitution(data,N,m,tau); %相空间重构 可将此段程序加到整个程序中，在时间循环内，可以保存时间序列的地方。

见完整程序M=N-(m-1)*tau; %重构相空间中相点的个数for i = 1 : (M-1)for j = i+1 : Md = 0;for k = 1 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d < dmax_d = d;endif min_d > d min_d = d;endavg_dd = avg_dd + d;endendavg_d = 2*avg_dd/(M*(M-1)); dlt_eps = (avg_d - min_d) * 0.02 ;点时，对max_eps 的放宽幅度 min_eps = min_d + dlt_eps / 2 ;max_eps = min_d + 2 * dlt_eps ;d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1)); endd = sqrt(d);if (d < DK) & (d > min_eps)DK = d;Loc_DK = i;endend%以下计算各相点对应的李氏数保存到lmd()数组中% i 为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK 为相点i-1对应最 短 距离的相点位置， DK 为其对应的最短距离% Loc_DK+1 为 Loc_DK 的演化点， DK1 为 i 点到 Loc_DK+1 点的距离，称为演化距离%平均相点距离%若在min_eps 〜max_eps 中找不到演化相离 DKDK = 1.0e+100; Loc_DK = 2; for i = (P+1):(M-1) d = 0;for k = 1 : m %第 i 个相点到其最近距离点的距离 %第 i 个相点对应的最近距离点的下标 %限制短暂分离，从点P+1开始搜索 %演化相点与当前相点距离的最小限%&&演化相点与当前相点距离的最大限% 从 P+1 〜 M-1 个相点中找与第一个相点最近的相点位置 (Loc_DK) 及其最短距%前i个Iog2 ( DK1/DK )的累计和用于求i点的lambda值sum」md = 0 ; %存放前i个Iog2( DK1/DK )的累计和for i = 2 : (M-1) % 计算演化距离DK1 = 0;for k = 1 : mDK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));endDK1 = sqrt(DK1); oId_Loc_DK = Loc_DK ; % 保存原最近位置相点oId_DK=DK;% 计算前i个Iog2(DK1/DK )的累计和以及保存i点的李氏指数if (DK1 ~= 0)&( DK ~= 0)sum_Imd = sum_Imd + Iog(DK1/DK) /Iog(2);endImd(i-1) = sum_Imd/(i-1); 此处可以保存不同相点i 对应的李氏指数，见完整程序。