-Lyapunov指數(shù)的計(jì)算方法

1、【總結(jié)】Lyapunov指數(shù)的計(jì)算方法非線性理論近期為了把計(jì)算LE的一些問題弄清楚，看了有79本書！下面以呂金虎混沌時(shí)間序列分析及其應(yīng)用、馬軍海復(fù)雜非線性系統(tǒng)的重構(gòu)技術(shù)為主線，把目前已有的LE計(jì)算方法做一個(gè)匯總！1.關(guān)于連續(xù)系統(tǒng)Lyapunov指數(shù)的計(jì)算方法連續(xù)系統(tǒng)LE的計(jì)算方法主要有定義方法、Jacobian方法、QR分解方法、奇異值分解方法，或者通過求解系統(tǒng)的微分方程，得到微分方程解的時(shí)間序列，然后利用時(shí)間序列（即離散系統(tǒng)）的LE求解方法來計(jì)算得到。關(guān)于連續(xù)系統(tǒng)LE的計(jì)算，主要以定義方法、Jacobian方法做主要介紹內(nèi)容。（1）定義法對(duì)用維連塊動(dòng)力系統(tǒng),二網(wǎng)力在時(shí)刻以孫為卬心,區(qū)為半徑膿

2、一個(gè)«維的球面.隨著時(shí)間的演化，在t時(shí)刻讀球面即變形為M縫的橢球面.設(shè)該橢球面的第i個(gè)坐標(biāo)軸方向的半軸長(zhǎng)為氏（與，則該系統(tǒng)第i個(gè)L河皿3指數(shù)為*% = hm-inI。t網(wǎng)（wM|阿有|此即連續(xù)系統(tǒng)LyaMnw指孩的定義.實(shí)際計(jì)向眇取忸武飛,0）|為d5為常敷），以中為球心，歐兒里痣范敬為d的正交矢量集（.,叫£為初始球,由非線性微分方程讓尸可以分別計(jì)菖出點(diǎn)Q町+如孫理小、足+%經(jīng)過時(shí)間才后演化的軌跡,設(shè)具終了專分別為了如、際八、杭，則令則國(guó)=6一.加fJ"=W3歷L則可得新的矢量集A4NJ由于各個(gè)矢量在演化過程中吉噲向著最大的UapurOT指數(shù)方向31撥，因此需要

3、通過Sthmidt正交化不斷地對(duì)新矢量進(jìn)行置換,即Wulf的文章中提出的GSR方法.表述如下：必就，齒嚴(yán)一<曲鏟及耀>蠟，壯；時(shí)刃料叫vcn5Kmk4入出支口）>以<>>以尸一附出叫排著以砒為球心，范毅為（I的正交矢境惠心？叫叱了;為新球繼續(xù)進(jìn)行演化，設(shè)演化至N步時(shí),得到矢量罪化必乜匕巧，且N足夠大,這可以得到Lyapunov才徽拊近做計(jì)皙公式.卜亨以各n4.4*一竽+$£貼.叫痛計(jì)亶時(shí),可以將d取為1定義法求解Lyapunov指數(shù).JPG關(guān)于定義法求解的程序，和matlab板塊的連續(xù)系統(tǒng)LE求解程序”差不多。以Rossler系統(tǒng)為例Rossler

4、系統(tǒng)微分方程定義程序functiondX=Rossler_ly(t,X)%Rossler吸引子，用索計(jì)算Lyapunov指數(shù)%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的三個(gè)列向量為相互正交的單位向量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

5、)=b+z*(x-c);%Rossler吸引子的Jacobi矩陣Jaco=0-1-1;1a0;z0x-c;dX(4:12)=Jaco*Y;求解LE代碼：%計(jì)算Rossler吸引子的Lyapunov指數(shù)clear;yinit=1,1,1;orthmatrix=100;010;001;a=0.15;b=0.20;c=10.0;y=zeros(12,1);%初始化輸入y(1:3)=yinit;y(4:12)=orthmatrix;tstart=0;%時(shí)間初始值tstep=1e-3;%時(shí)間步長(zhǎng)wholetimes=1e5;%總的循環(huán)次數(shù)steps=10;%每次演化的步數(shù)iteratetimes=who

6、letimes/steps;%演化的次數(shù)mod=zeros(3,1);lp=zeros(3,1);%初始化三個(gè)Lyapunov指數(shù)Lyapunov1=zeros(iteratetimes,1);Lyapunov2=zeros(iteratetimes,1);Lyapunov3=zeros(iteratetimes,1);fori=1:iteratetimestspan=tstart:tstep:(tstart+tstep*steps);T,Y=ode45('Rossler_ly',tspan,y);%取積分得到的最后一個(gè)時(shí)刻的值y=Y(size(Y,1),:);%重新定義起始時(shí)

7、刻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);%取三個(gè)向量的模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);%三個(gè)Lyapunov指數(shù)L

8、yapunovl(i)=lp(1)/(tstart);Lyapunov2(i)=lp(2)/(tstart);Lyapunov3(i)=lp(3)/(tstart);y(4:12)=y0'end%作Lyapunov指數(shù)譜圖i=1:iteratetimes;plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)程序中用到的ThreeGS程序如下：%G-S正交化functionA=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

9、=v1;a2=v2-(a1'*v2)/(a1'*a1)*a1;a3=v3-(aT*v3)/(aT*a1)*a1-(a2'*v3)/(a2'*a2)*a2;A=a1,a2,a3;計(jì)算得到的Rossler系統(tǒng)的LE為0.0632310.092635-9.8924Wolf文章中計(jì)算得到的Rossler系統(tǒng)的LE為0.090-9.77需要注意的是一一定義法求解的精度有限，對(duì)有些系統(tǒng)的計(jì)算往往出現(xiàn)計(jì)果和理論值有偏差的現(xiàn)象。正交化程序可以根據(jù)上面的擴(kuò)展到N*N向量，這里就不加以說明了，對(duì)matlab用戶來說應(yīng)該還是比較簡(jiǎn)單的！Jacobian方法通過資料檢索，發(fā)現(xiàn)論壇中用的

10、較多的LET工具箱的算法原理就是Jacobian方法。基本原理就是首先求解出連續(xù)系統(tǒng)微分方程的近似解，然后對(duì)系統(tǒng)的Jacobian矩陣進(jìn)行QR分解，計(jì)算Jacobian矩陣特征值的乘積，最后計(jì)算出LE和分?jǐn)?shù)維。經(jīng)過計(jì)算也證明了這種方法精度較高，對(duì)目前常見的混沌系統(tǒng)，如Lorenz、Henon、Duffing等的Lyapunov指數(shù)的計(jì)算精度都很高，而且程序編寫有一定的規(guī)范，個(gè)人很推薦使用。（雖然我自己要做的系統(tǒng)并不適用于孑）LET工具箱可以在網(wǎng)絡(luò)上找到，這里就不列出了！關(guān)于LET工具箱如果有問題，歡迎加入本帖討論！Jacobian法求解Lyapunov指數(shù).JPG對(duì)離散動(dòng)力系統(tǒng)，或者說是非線性

11、時(shí)間序列，往往不需要計(jì)算出所有的Lyapunov指數(shù)，通常只需計(jì)算出其最大的Lyapunov指數(shù)即可?！?983,格里波基證明了只要最大Lyapunov指數(shù)大于零，就可以肯定混沌的存在”。目前常用的計(jì)算混沌序列最大Lyapunov指數(shù)的方法主要有一下幾種：(1)由定義法延伸的Nicolis方法(2) Jacobian方法(3) 3)Wolf方法(4) P范數(shù)方法(5)小數(shù)據(jù)量方法其中以Wolf方法和小數(shù)據(jù)量方法應(yīng)用最為廣泛，也最為普遍。下面對(duì)Nicolis方法、Wolf方法以及小數(shù)據(jù)量方法作一一介紹。1) 1)Nicolis方法這種方法和連續(xù)系統(tǒng)的定義方法類似，而且目前應(yīng)用很有限制，因此只對(duì)其

12、理論進(jìn)行介紹，編程應(yīng)用方面就省略了Nicolis方法求最大LE.JPG2) Wolf方法Wolf方法求最大LE.JPGWolf方法的Matlab程序如下：functionlambda_1=lyapunov_wolf(data,N,m,tau,P)%該函數(shù)用來計(jì)算時(shí)間序列的最大Lyapunov指數(shù)-Wolf方法%m:嵌入維數(shù)！一般選大于等于10%tau:時(shí)間延遲！一般選與周期相當(dāng)，如我選2000%data:時(shí)間序列！可以選1000;%N:時(shí)間序列長(zhǎng)度滿足公式：M=N-(m-1)*tau=24000-(10-1)*1000=5000%P:時(shí)間序列的平均周期，選擇演化相點(diǎn)距當(dāng)前點(diǎn)的位置差，即若當(dāng)前相

13、點(diǎn)為I,則演化相點(diǎn)只能在|IJ卜P的相點(diǎn)中搜尋！P=周期=600%lambda_1:返回最大lyapunov指數(shù)值min_point=1;%&&要求最少搜索到的點(diǎn)數(shù)MAX_CISHU=5;%&&最大增加搜索范圍次數(shù)%FLYINGHAWK%求最大、最小和平均相點(diǎn)距離%最大相點(diǎn)距離%最小相點(diǎn)距離max_d=0;min_d=1.0e+100;avg_dd=0;Y=reconstitution(data,N,m,tau);%相空間重構(gòu)可將此段程序加到整個(gè)程序中，在時(shí)間循環(huán)內(nèi)，可以保存時(shí)間序列的地方。見完整程序。M=N-(m-1)*tau;%重構(gòu)相空間中相點(diǎn)的個(gè)數(shù)fori

14、=1:(M-1)forj=i+1:Md=0;fork=1:md=d+(Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);endd=sqrt(d);ifmax_d<dmax_d=d;endifmin_d>dmin_d=d;endavg_dd=avg_dd+d;endend%Sf在 min_epsmax_eps中找不至U演化相%演化相點(diǎn)與當(dāng)前相點(diǎn)距離的最小限%&&演化相點(diǎn)與當(dāng)前相點(diǎn)距離的最大限avg_d=2*avg_dd/(M*(M-1);%平均相點(diǎn)距離dlt_eps=(avg_d-min_d)*0.02;點(diǎn)時(shí)，對(duì)max_eps的放寬幅度min_eps=min_

15、d+dlt_eps/2;max_eps=min_d+2*dlt_eps;%從p+iM-1個(gè)相點(diǎn)中找與第一個(gè)相點(diǎn)最近的相點(diǎn)位置(Loc_DK)及其最短距離DKDK=1.0e+100;%第i個(gè)相點(diǎn)到其最近距離點(diǎn)的距離Loc_DK=2;%第i個(gè)相點(diǎn)對(duì)應(yīng)的最近距離點(diǎn)的下標(biāo)fori=(P+1):(M-1)%限制短暫分離，從點(diǎn)P+1開始搜索d=0;fork=1:md=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%以下計(jì)算各相點(diǎn)對(duì)應(yīng)的李氏數(shù)保存到lmd()數(shù)組中%

16、i為相點(diǎn)序號(hào)，從1至iJ(M-1),也是i-1點(diǎn)的演化點(diǎn)；Loc_DK為相點(diǎn)i-1對(duì)應(yīng)最短距離的相點(diǎn)位置，DK為其對(duì)應(yīng)的最短距離%Loc_DK+1為L(zhǎng)oc_DK的演化點(diǎn)，DK1為i點(diǎn)到Loc_DK+1點(diǎn)的距離，稱為演化距離%前i個(gè)log2(DK1/DK)的累計(jì)和用于求i點(diǎn)的lambda值sum_lmd=0;%存放前i個(gè)log2(DK1/DK)的累計(jì)和fori=2:(M-1)%計(jì)算演化距離DK1=0;fork=1:mDK1=DK1+(Y(k,i)-Y(k,Loc_DK+1)*(Y(k,i)-Y(k,Loc_DK+1);endDK1=sqrt(DK1);old_Loc_DK=Loc_DK;%保存原

17、最近位置相點(diǎn)old_DK=DK;%計(jì)算前i個(gè)log2(DK1/DK)的累計(jì)和以及保存i點(diǎn)的李氏指數(shù)if(DK1=0)&(DK=0)sum_lmd=sum_lmd+log(DK1/DK)/log(2);endlmd(i-1)=sum_lmd/(i-1);此處可以保存不同相點(diǎn)i對(duì)應(yīng)的李氏指數(shù)，見完整程序。%以下尋找i點(diǎn)的最短距離:要求距離在指定距離范圍內(nèi)盡量短,與DK1的角度最小point_num=0;%&&在指定距離范圍內(nèi)找到的候選相點(diǎn)的個(gè)數(shù)cos_sita=0;%&&夾角余弦的比較初值要求一定是銳角zjfwcs=0;%&&增加范圍次數(shù)wh

18、ile(point_num=0)%*搜索相點(diǎn)forj=1:(M-1)ifabs(j-i)<=(P-1)%&&候選點(diǎn)距當(dāng)前點(diǎn)太近，跳過！continue;end%*計(jì)算候選點(diǎn)與當(dāng)前點(diǎn)的距離dnew=0;fork=1:mdnew=dnew+(Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);enddnew=sqrt(dnew);if(dnew<min_eps)|(dnew>max_eps)%&&不在距離范圍，跳過！continue;end%*計(jì)算夾角余弦及比較DOT=0;fork=1:mDOT=DOT+(Y(k,i)-Y(k,j)*(Y(k

19、,i)-Y(k,old_Loc_DK+1);endCTH=DOT/(dnew*DK1);ifacos(CTH)>(3.14151926/4)%&&不是小于45度的角，跳過！continue;endifCTH>cos_sita%&&新夾角小于過去已找到的相點(diǎn)的夾角，保留cos_sita=CTH;Loc_DK=j;DK=dnew;endpoint_num=point_num+1;endifpoint_num<=min_pointmax_eps=max_eps+dlt_eps;zjfwcs=zjfwcs+1;ifzjfwcs>MAX_CISHU

20、%&&超過最大放寬次數(shù)，改找最近的點(diǎn)DK=1.0e+100;forii=1:(M-1)ifabs(i-ii)<=(P-1)%&&候選點(diǎn)距當(dāng)前點(diǎn)太近，跳過！continue;endd=0;fork=1:md=d+(Y(k,i)-Y(k,ii)*(Y(k,i)-Y(k,ii);endd=sqrt(d);if(d<DK)&(d>min_eps)DK=d;Loc_DK=ii;endendbreak;endpoint_num=0;%&&擴(kuò)大距離范圍后重新搜索cos_sita=0;endendend%取平均得到最大李雅普諾夫指數(shù)(此

21、處只有一個(gè)值，若為正說明體系是混沌的，若為負(fù)則說明體系是周期性的確定性運(yùn)動(dòng))lambda_1=sum(lmd)/length(lmd);程序中用到的reconstitution函數(shù)如下：此段程序可直接放在時(shí)間循環(huán)內(nèi)部，即可以保存時(shí)間序列的地方。見完整程序范例。functionX=reconstitution(data,N,m,tau)%該函數(shù)用來重構(gòu)相空間%m為嵌入空間維數(shù)%tau為時(shí)間延遲%data為輸入時(shí)間序列%N為時(shí)間序列長(zhǎng)度%X為輸出,是m*n維矩陣M=N-(m-1)*tau;%相空間中點(diǎn)的個(gè)數(shù)forj=1:M%相空間重構(gòu)fori=1:mX(i,j)=data(i-1)*tau+j);

22、endend這里聲明一下，這些程序并非我自己編寫的，均是轉(zhuǎn)載，其使用我已經(jīng)驗(yàn)證過，絕對(duì)可以運(yùn)行！(3)小數(shù)據(jù)量方法說小數(shù)據(jù)量方法是目前最實(shí)用、應(yīng)用最廣泛的方法應(yīng)該不為過吧，呵呵！小數(shù)據(jù)量方法求最大Lyapunov指數(shù).JPG上面兩種方法，即Wolf方法和小數(shù)據(jù)量方法，在計(jì)算LE之前，都要求對(duì)時(shí)間序列進(jìn)行重構(gòu)相空間，重構(gòu)相空間的優(yōu)良對(duì)于最大LE的計(jì)算精度影響非常大！因此重構(gòu)相空間的幾個(gè)參數(shù)的確定就非常重要。(1)時(shí)間延遲主要推薦兩種方法自相關(guān)函數(shù)法、CC方法自相關(guān)函數(shù)法對(duì)一個(gè)混沌時(shí)間序列，可以先寫出其自相關(guān)函數(shù)，然后作出自相關(guān)函數(shù)關(guān)于時(shí)間t的函數(shù)圖像。根據(jù)數(shù)值試驗(yàn)結(jié)果，當(dāng)自相關(guān)函數(shù)下降到初始值的

23、11/e時(shí)，所得的時(shí)間t即為重構(gòu)相空間的時(shí)間延遲。CC方法可以同時(shí)計(jì)算出時(shí)間延遲和時(shí)間窗口，個(gè)人推薦使用這種方法！(2)平均周期平均周期的計(jì)算可以采用FFT方法。在matlab幫助中有一個(gè)太陽(yáng)黑子的例子，現(xiàn)摘錄如下：loadsunspot.dat%裝載數(shù)據(jù)文件year=sunspot(:,1);%讀取年份信息%讀取黑子活動(dòng)數(shù)據(jù)%繪制原始數(shù)據(jù)圖% 快速 FFT 變換wolfer=sunspot(:,2);plot(year,wolfer)title('SunspotData')Y=fft(wolfer);N=length(Y);%FFT變換后數(shù)據(jù)長(zhǎng)度Y(1)=;%去掉Y的第一個(gè)數(shù)

24、據(jù)，它是所有數(shù)據(jù)的和power=abs(Y(1:N/2).A2;%求功率譜nyquist=1/2;freq=(1:N/2)/(N/2)*nyquist;%求頻率plot(freq,power),gridon%繪制功率譜圖xlabel('cycles/year')title('Periodogram')period=1./freq;%年份(周期)plot(period,power),axis(04002e7),gridon繪制年份功率譜曲線ylabel('Power')xlabel('Period(Years/Cycle)')mp,

25、index=max(power);%求最高譜線所對(duì)應(yīng)的年份下標(biāo)period(index)%由下標(biāo)求出平均周期(3)嵌入維數(shù)目前嵌入維數(shù)的主要計(jì)算方法是采用Grassberge和Procaccia提出的GP算法計(jì)算出序列的關(guān)聯(lián)維數(shù)d,然后利用嵌入維數(shù)m>=2d+1,選取合適的嵌入維數(shù)。GP算法程序如下：functionln_r,ln_C=G_P(data,N,tau,min_m,max_m,ss)%thefunctionisusedtocalculatecorrelationdimentionwithG-Palgorithm%計(jì)算關(guān)聯(lián)維數(shù)的GP算法%data:thetimeseries時(shí)間

26、序列時(shí)間序列長(zhǎng)度 時(shí)間延遲 最小的嵌入維數(shù)%N:thelengthofthetimeseries%tau:thetimedelay%min_m:theleastembeddeddimentionm%max_m:thelargestembeddeddimentionm最大的嵌入維數(shù)%ss:thestepsizeofr的步長(zhǎng)%skyhawkform=min_m:max_mY=reconstitution(data,N,m,tau);%reconstitutestatespaceM=N-(m-1)*tau;%thenumberofpointsinstatespacefori=1:M-1forj=i+

27、1:Md(i,j)=max(abs(Y(:,i)-Y(:,j);%calculatethedistanceofeachtwoend%pointsinstatespacd算狀態(tài)空間中每?jī)牲c(diǎn)之間的距離endmax_d=max(max(d);%themaxdistanceofallpoints得到所有點(diǎn)之間的最大距離d(1,1)=max_d;min_d=min(min(d);%themindistanceofallpoints得到所有點(diǎn)間的最短距離delt=(max_d-min_d)/ss;%thestepsizeofr得到r的步長(zhǎng)fork=1:ssr=min_d+k*delt;C(k)=corre

28、lation_integral(Y,M,r);%calculatethecorrelationintegralln_C(m,k)=log(C(k);%lnC(r)ln_r(m,k)=log(r);%lnrfprintf('%d/%d/%d/%dn',k,ss,m,max_m);endplot(ln_r(m,:),ln_C(m,:);holdon;endfid=fopen('lnr.txt','w');fprintf(fid,'%6.2f%6.2fn',ln_r);fclose(fid);fid=fopen('lnC.txt

29、','w');fprintf(fid,'%6.2f%6.2fn',ln_C);fclose(fid);程序中的correlation_integral函數(shù)如下：functionC_I=correlation_integral(X,M,r)%thefunctionisusedtocalculatecorrelationintegral%C_I:thevalueofthecorrelationintegral%X:thereconstitutedstatespace,Misam*Mmatrix%m:theembeddingdemention%M:Misthenumberofembeddedpointsinm-dimensionalsapce%r:theradiusoftheHeavisidefunction,sigma/2<r<2sigma%calculatethesumofallthevaluesofHeaviside%skyhawksum_H=0;fori=1:M%fprintf('%d/