Lyapunov指数的计算方法

《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 = -yz，% 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

3、（9), X（12）；% 输出向量的初始化，必不可少dX = zeros（12，1）;% Rossler吸引子dX（1) = -yz；dX(2) = x+a*y;dX（3) = b+z*(xc）；% Rossler吸引子的Jacobi矩阵Jaco = 0 -1 -1； 1 a0； z 0x-c；dX(4:12) = JacoY；求解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（

4、1:3） = yinit；y(4：12） = orthmatrix；tstart = 0； % 时间初始值tstep = 1e3; % 时间步长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(iteratetime

5、s，1）；for i=1：iteratetimes tspan = tstart：tstep：(tstart + tstepsteps); T，Y = ode45(Rossler_ly， tspan, y）； % 取积分得到的最后一个时刻的值 y = Y(size（Y,1),:）； 重新定义起始时刻 tstart = tstart + tstepsteps; 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)）; m

6、od（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 =

7、 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-（(a1v3)/(a1a1)）*a1（(a2v3)/（a2*a2)）a2;A = a1,a2，a3；计算得到的Ro

8、ssler系统的LE为-0。0632310。0926359。8924Wolf文章中计算得到的Rossler系统的LE为-0。09 0 -9.77需要注意的是定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象.正交化程序可以根据上面的扩展到N*N向量,这里就不加以说明了，对matlab用户来说应该还是比较简单的！（2)Jacobian方法通过资料检索，发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法。基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维。经过

9、计算也证明了这种方法精度较高,对目前常见的混沌系统，如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高，而且程序编写有一定的规范，个人很推荐使用。（虽然我自己要做的系统并不适用)LET工具箱可以在网络上找到，这里就不列出了!关于LET工具箱如果有问题,欢迎加入本帖讨论！Jacobian法求解Lyapunov指数.JPG对离散动力系统，或者说是非线性时间序列，往往不需要计算出所有的Lyapunov指数，通常只需计算出其最大的Lyapunov指数即可。“1983年,格里波基证明了只要最大Lyapunov指数大于零，就可以肯定混沌的存在”。目前常用的计算混沌序列最大L

10、yapunov指数的方法主要有一下几种:(1）由定义法延伸的Nicolis方法（2）Jacobian方法（3)Wolf方法(4)P范数方法(5)小数据量方法其中以Wolf方法和小数据量方法应用最为广泛，也最为普遍。下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍。（1)Nicolis方法这种方法和连续系统的定义方法类似，而且目前应用很有限制，因此只对其理论进行介绍，编程应用方面就省略了Nicolis方法求最大LE.JPG（2）Wolf方法Wolf方法求最大LE。JPGWolf方法的Matlab程序如下:function lambda_1=lyapunov_wolf(data,N

11、,m,tau，P)该函数用来计算时间序列的最大Lyapunov 指数-Wolf 方法m： 嵌入维数! 一般选大于等于10%tau:时间延迟 ！一般选与周期相当，如我选2000%data：时间序列！可以选1000；N：时间序列长度满足公式：M=N（m-1）tau=24000-(101)1000=5000%P：时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I,则演化相点只能在IJ|P的相点中搜寻 ! P=周期=600%lambda_1:返回最大lyapunov指数值min_point=1； &要求最少搜索到的点数MAX_CISHU=5 ；%最大增加搜索范围次数FLYINGHAWK

12、% 求最大、最小和平均相点距离 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 : (M1) 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 min_ep

13、s） 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 : (M1） % 计算演化距离 DK1 = 0； for k = 1 : m DK1 = DK1 + （Y

14、(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对应的李氏指数，见完整程序.。%以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小 po

15、int_num = 0； % &在指定距离范围内找到的候选相点的个数 cos_sita = 0； &夹角余弦的比较初值 要求一定是锐角 zjfwcs=0;%&增加范围次数 while （point_num = 0） * 搜索相点 for j = 1 : （M-1) if abs(ji） =(P1) &候选点距当前点太近，跳过！ 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 max_eps

16、） %&不在距离范围，跳过! 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

17、 +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 ： (M1) if abs(iii） = （P1） &候选点距当前点太近，跳过！ 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

18、) 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为嵌入

19、空间维数% tau为时间延迟% data为输入时间序列% N为时间序列长度 X为输出,是m*n维矩阵M=N-(m-1)tau；相空间中点的个数for j=1:M %相空间重构 for i=1：m X(i,j)=data(i1)*tau+j）; endend这里声明一下，这些程序并非我自己编写的，均是转载，其使用我已经验证过，绝对可以运行！（3）小数据量方法说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧，呵呵！小数据量方法求最大Lyapunov指数。JPG上面两种方法，即Wolf方法和小数据量方法,在计算LE之前,都要求对时间序列进行重构相空间，重构相空间的优良对于最大LE的计算精度影

20、响非常大!因此重构相空间的几个参数的确定就非常重要。（1）时间延迟主要推荐两种方法-自相关函数法、CC方法自相关函数法对一个混沌时间序列,可以先写出其自相关函数，然后作出自相关函数关于时间t的函数图像。根据数值试验结果,当自相关函数下降到初始值的11/e时,所得的时间t即为重构相空间的时间延迟。CC方法-可以同时计算出时间延迟和时间窗口，个人推荐使用这种方法!（2)平均周期平均周期的计算可以采用FFT方法。在matlab帮助中有一个太阳黑子的例子，现摘录如下：load sunspot.dat 装载数据文件year = sunspot（：，1）； %读取年份信息wolfer = sunspot(

21、:，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

22、； 年份（周期）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提出的GP算法计算出序列的关联维数d，然后利用嵌入维数m=2d+1，选取合适的嵌入维数.GP算法程序如下：function ln_r，ln_C=G_P(data,N,tau,min_m,m

23、ax_m，ss)% the function is used to calculate correlation dimention with GP algorithm% 计算关联维数的GP算法 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

24、of r r的步长skyhawkfor 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:M1 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 ma

25、x 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)；

26、lnr fprintf(d/d/d/dn,k,ss,m，max_m）; end plot(ln_r（m,:），ln_C（m，:）)； hold on;endfid=fopen(lnr。txt,w);fprintf(fid，6.2f %6。2fn，ln_r）;fclose（fid）;fid = fopen（lnC。txt，w)；fprintf(fid,%6.2f %6.2fn,ln_C);fclose(fid);程序中的correlation_integral函数如下:function C_I=correlation_integral(X,M，r)the function is used to

27、calculate correlation integralC_I:the value of the correlation integralX：the reconstituted state space,M is a m*M matrixm:the embedding dementionM:M is the number of embedded points in mdimensional sapce%r:the radius of the Heaviside function，sigma/2r2sigma%calculate the sum of all the values of Hea

28、visideskyhawksum_H=0；for i=1:M fprintf（d/%dn，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； endendC_I=2*sum_H/（M(M1）;%the value of correlation integra

29、l无水1324发布于2007-0803 20:35：01感谢octopussheng的总结！这是现有的方法的一个总结,不知道你对这些方法有些什么体会，或者说他们的局限，现在还有作这方面的改进的吗?先总结总结，具体的比较会逐步贴出来的，而且这也需要大家的一起努力啊！光是我在这里说，也只是我自己的一点体会,希望用过这些方法的一起来参与这个总结工作!以上的各种方法在实际应用的时候要根据具体情况来选择。一般地，如果已知系统方程(当然系统不能太过复杂)时，则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些！而如果系统方程比较复杂(如超维系统）、或者为一时间序列时,则推荐采样Wolf

30、方法、小数据量方法。Wolf方法的特点是时间序列无噪声，空间中小向量的演变高度非线性，而Jacobian方法则是噪声大,空间中小向量的演变接近线性。小数据量方法的优点在于：（1）对小数据组的计算可靠；（2）计算量较小，比wolf方法快很多;（3）编程、操作较为容易。而关于时间延迟、嵌入维数、平均周期的确定，还是推荐CC方法和GP算法,结果更为可靠一些!你这已经做得很好了。这个具体的应用我也是做得很少的,还有请注意一下matlab版中的一个帖子, 早就注意了，我这里面一些代码就是从那边套用的！试过了吗？都没有问题吧计算起来，原来论坛里面一些程序不能运行,套用，那你就要注意，在这里重写就要有我们的特色。不能够是简单的重请详细读帖子，上面有说明,当然不是简单的重复啦，我这里是理论结合实践,呵呵！程序都测试能算的,不过对不同的系统，精度不敢保证的太高,呵呵，只要能运行就可以参照自己来编写精度可以根据自己需要来选择呵呵,确实,保证能够使用！希望这样总结一下对大家有所帮助吧！非常有价值。谢谢！