《-Lyapunov指数的计算方法[互联网+]》由会员分享，可在线阅读，更多相关《-Lyapunov指数的计算方法[互联网+]（12页珍藏版）》请在金锄头文库上搜索。
1、【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有79本书！下面以吕金虎混沌时间序列分析及其应用、马军海复杂非线性系统的重构技术为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法 连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。（1）定义法定义法求解Lyapunov指数.JPG 关于定义法
2、求解的程序，和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(
3、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 0x-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
4、);% 初始化输入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
5、= 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)
6、= 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)
7、 = 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*
8、a2)*a2;A = a1,a2,a3;计算得到的Rossler系统的LE为0.0632310.092635-9.8924Wolf文章中计算得到的Rossler系统的LE为0.09 0 -9.77需要注意的是定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。正交化程序可以根据上面的扩展到N*N向量，这里就不加以说明了，对matlab用户来说应该还是比较简单的！（2）Jacobian方法 通过资料检索，发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法。基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jaco
9、bian矩阵特征值的乘积，最后计算出LE和分数维。经过计算也证明了这种方法精度较高，对目前常见的混沌系统，如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高，而且程序编写有一定的规范，个人很推荐使用。（虽然我自己要做的系统并不适用 ）LET工具箱可以在网络上找到，这里就不列出了！关于LET工具箱如果有问题，欢迎加入本帖讨论！Jacobian法求解Lyapunov指数.JPG 对离散动力系统，或者说是非线性时间序列，往往不需要计算出所有的Lyapunov指数，通常只需计算出其最大的Lyapunov指数即可。“1983年，格里波基证明了只要最大Lyapunov指数大
10、于零，就可以肯定混沌的存在”。目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种：（1）由定义法延伸的Nicolis方法（2）Jacobian方法（3）Wolf方法（4）P范数方法（5）小数据量方法其中以Wolf方法和小数据量方法应用最为广泛，也最为普遍。下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍。（1）Nicolis方法 这种方法和连续系统的定义方法类似，而且目前应用很有限制，因此只对其理论进行介绍，编程应用方面就省略了 Nicolis方法求最大LE.JPG（2）Wolf方法 Wolf方法求最大LE.JPGWolf方法的Matlab程序如下：functi
11、on lambda_1=lyapunov_wolf(data,N,m,tau,P)%该函数用来计算时间序列的最大Lyapunov 指数-Wolf 方法%m: 嵌入维数 ！ 一般选大于等于10%tau:时间延迟 ！一般选与周期相当，如我选2000%data:时间序列 ！可以选1000；%N:时间序列长度 满足公式：M=N-(m-1)*tau=24000-(10-1)*1000=5000%P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|IJ|P的相点中搜寻 ! P=周期=600%lambda_1: 返回最大lyapunov指数值min_point=1; %
12、&要求最少搜索到的点数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 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 ; %若在
