Lyapunov指数的计算方法

余年寄山水
779次浏览
2021年02月20日 06:29
最佳经验
本文由作者推荐

-

2021年2月20日发(作者：同喜同乐)


【总结】


Lyapunov


指数的计算方法



非线性理论



近 期为了把计算


LE


的一些问题弄清楚，看了有

< br>7


〜


9


本书！下面以吕金虎《混 沌



时间


序列分析及其应用》、马军海 《复杂非线性系统的重构技术》为主线，把目前



已有的


LE


计算方法做一个汇总！



1.


关于连续系统


Lyapunov

< br>指数的计算方法



连续系统


LE


的计算方法主要有定义



方

< p>
法、


Jacobian


方法、

QR


分解方法、奇异值分解方法，或者通过求解系统的微分



方


程，得到微分方程解的时间序列，然后利用时间序列（即离散 系统）的



LE


求解



方法


来计算得到。关于连续系统



LE


的计算，主要以定义方法、


Jacob ian


方法做主



要介绍


内容。



（


1


）定义法



—



对


H


堆连


续动力系統


z =

< br>在


—


OBJ


孙

< br>


味“为


中心


.|


拆（心


0


）


||

< p>
为丰笹啟存



«


堆的球面


*


施著时间的演化，在


t


时


討


该球而


0P

< p>
变形


为


M


继的椭球厨・设 滾椭域面的第


/


个坐


标轴方向的半< /p>


轴长対


卩


兀


|< /p>


，


则诙


系统第


i


个


指


数対


*< /p>




此即连续系统


Lyapunov


揩较


飽定冥


・



弼计尊时・取


|


处


（心


0


）［


为岀


W


为常


数），以孔


为球心・欧几里篇范敢为山的正衮



矢


量集仙测


，…叮


为初始球


.< /p>


由非线性徴分方崔


“尸


㈤


可以


分别计


算出点細



血


创、



引他 、


r


引址


经过时间

t


后淺化的轨


迹・役


其终了


点分别为珊


、


砒、


f


仙



则令



石


f


陶一心■处

严


=


甩


-


和，


r


亦耳国


=

略一報


#


则可


得新的矢重棄



叶禺巴


…后畀｝


・



由于各牛妥臺在演


化过


程


中舌焙


向着是大的


Uap urOT IS


数方何


靠掘，因此需要通过



Schimdt IE


交化不断地讨新矢量逬行置换

< p>
.


SP Wolf to


文章中

< br>提出的


GSR^


法


.

< p>
表述如


下：




播着以他为球


心，


疤数対


（


I


的


正奁矢


臺料创


巴叫叫…伽严


;


为

< p>
祈球继续进行演



出



设演化


1




至


N


步时，得到矢董慕


冈


㈤


出巴…耳僅牛


且


N


足够大


,


这可以得到


Lyapunov


扌鐵的


近似计


算


公式三




实际计算时，取为


1


・


定义法求解



Lyapunov


指数



JPG


关于定义法求解的程序，和



matlab


板块的



连续系统


LE


求解程序”差不多。以

< br>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;


2





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

< br>程序如下：



%G-S


正交化



function A = ThreeGS(V) % V


为



3*3


向量



v1 = V(:,1);


3




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.8924


Wolf


文章中计算得到的



Rossler


系统的


LE

为



--------


0.09 0 -9.77


需要注意的是一一定义法求解的精度有限，对有些系统的计算往往出现 计果和理论



值有偏


差的现象。



正交化程序可以根据上面的扩展到



N*N


向量，这里就不加以说明了，对



matlab


用



户来


说应该还是比较简单的！



(2


)


Jacobian


方法



通过资料检索，发现论坛中用的较多的



LET


工具箱的算法原理就是


Jacobian


方



法。基本原理就是首先求解出连续系统微分方程 的近似解，然后对系统的



Jacobian



矩阵进行


QR


分解，计算


Jacobian


矩阵特征值的乘积，最后计算出



LE


和分数维。


经过计算也证明了这种方法精度较高，对目前常见的混沌系统，如



Lore nz


、


< br>Henon


、


Duffing


等 的


Lyapunov


指数的计算精度都很高，而且程序编写有一 定的规



范，


个人很推荐使用。（虽然 我自己要做的系统并不适用



罢


@


）




4




LET


工具箱可以在网络上找到，这里就不列出了！关于



LET


工具箱如果有问题，



欢迎加


入本帖讨论！



Jacobian


法求解



Lyapunov


指数



JPG


对离散动力系统，或者说是非线性时间序列，往往


不需要计算出所有的



Lyapunov< /p>


指数，通常只需计算出其最大的



Lyap u nov


指数即可。“


1 98


年，格里波



基证明了只要最大< /p>


Lyapunov


指数大于零，就可以肯定混沌的存在”。



目前常用的计算混沌序列最大



Lyapu nov


指数的方法主要有一下几种：



（


1


）由定义法延伸的



Nicolis


方法



（


2


）


Jacobian< /p>


方法



（



3


）



Wolf


方法


（


4


）


P


—范数方法



（



5


）小数据量方法



其中以



Wolf


方法和小数据量方法应用最为广泛，也最为普遍。


< p>
下面对


Nicolis


方法、

Wolf


方法以及小数据量方法作



介绍。



（



1


）


Nicolis


方法



这种方法和连续系统的定义方法 类似，而且目前应用很有限制，因此只对其理论



进行介


绍，编程应用方面就省略了



Nicolis


方法求最大





2


）


Wolf


方法



Wolf


方法求最大




Wolf


方法的



Matlab


程序如下：



function lambda_1=lyapunov_wolf(data,N,m,tau,P) %


该函数用来计算时间序列的最大



Lyapunov


指数


--Wolf


方法



！



一般选大于等于



10



% m:


嵌入维数



% tau


时间延迟



！



一般选与周期相当，如我选



2000


！可以选



1000


；



% data


时间序列



% N:


时间序列长度



满足公式：


M=N-


（


m- 1


）


*tau=24000-


（


10-1


）


*1000=5000


% P


：


时间序列的平均周期


，


选择演化相点距当前点的位置差，即若当前相点为


I


，则



演化相点只能在


|I


—



J|>P


的相点中搜寻



！


P=


周期


=600



% lambda_1:


返回最大


l yapunov


指数值



min_point=1


%&&


要求最少搜索到的点数



MAX_CISHU=5 %&&


最大增加搜索范围次数



%FLYINGHAWK


%


求最大、最小和平均相点距离



%


最大相点距


max_d = 0; min_d = 1.0e+100;



离



%


最小相 点距


5




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


点时，对


max_eps


的放宽幅度

< p>


%


若在


min_eps


〜


max_eps


中找不到演化


相



%


演化相点与当前相点距离的最小限



%&&


演化相点与当前相点距离的最大限



%


从



P+1


〜



M-1


个相点中找与第一个相点最近的相点位置



(Loc_DK)


及其最短距



离



DK


DK = 1.0e+100;


%


第



i


个相点到其最近距离点的距离



Loc_DK = 2;


%


第



i


个相点对应的最近距离点的下标



%


限制短暂分离，从


for i = (P+1):(M-1) d


点


P+1


开始搜索



= 0;


for k = 1 : m



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)); end


d = sqrt(d);


if (d < DK) & (d > min_eps)


DK = d;


Loc_DK = i;


end


end


%

< br>以下计算各相点对应的李氏数保存到


lmd()


数组中< /p>



% i


为相点序号，从


1


到


(M-1)


，也是


i-1


点的演化点；


Loc_DK

< p>
为相点


i-1


对应最



短



距离的相点位置，



DK


为其对应的最短距离



%


Loc_DK+1


为



Loc_DK


的演化点，



DK1


为



i


点到



Loc_DK+1


点的距离，称为演化距离



6



-


-


-


-


-


-


-


-