matlab編寫的Lyapunov指數(shù)計算程序匯總

1、matlab編寫的lyapunov指數(shù)計算程序匯總 當(dāng)前離線 威望45 點(diǎn)c3p幣59779 元貢獻(xiàn)值2318 點(diǎn)推廣邀請能量5815 焦耳最后登錄2012-3-10注冊時間2000-7-9帖子42852精華39積分59779閱讀權(quán)限200uid21ip卡 狗仔卡 管理員修身 齊家 賺錢 玩天下威望45 點(diǎn)c3p幣59779 元貢獻(xiàn)值2318 點(diǎn)推廣邀請能量5815 焦耳最后登錄2012-3-10注冊時間2000-7-9帖子42852精華39積分59779閱讀權(quán)限200uid21 串個門 加好友 打招呼 發(fā)消息主題帖 發(fā)表于 2009-2-22 20:16:50 |只看該作者 以下整理的內(nèi)容，

2、版權(quán)歸原作者，轉(zhuǎn)自。申明：以下各程序?yàn)閭€人在網(wǎng)絡(luò)上收集的lyapunov指數(shù)計算程序，未經(jīng)過驗(yàn)證，不保證程序的正確性和計算結(jié)果的正確性，請大家見諒，也歡迎大家探討！計算連續(xù)方程lyapunov指數(shù)的程序，比較好用的其中給出了兩個例子：計算rossler吸引子的lyapunov指數(shù)計算超混沌rossler吸引子的lyapunov指數(shù)/downloads39/sourcecode/math/detail132231.html=基于volterra濾波器混沌時間序列多步預(yù)測-參考文獻(xiàn)：1、張家樹.混沌時間序列的volterra自適應(yīng)預(yù)測

3、.物理學(xué)報.2000.032、scott c.douglas, teresa h.-y. meng, normalized data nonlinearities for lms adaptation. ieee trans.sign.proc. vol.42 1994 -文件說明：1、original_multisteppred_main.m 程序主文件,直接運(yùn)行此文件即可2、original_train.m 訓(xùn)練函數(shù)3、original_test.m 測試函數(shù)4、lorenzdata.dll 產(chǎn)生lorenz離散序列5、normalize_1.m 歸一化6、phasparecon.m 相空

4、間重構(gòu)7、phaspa2voltcoef.m 構(gòu)造 volterra 自適應(yīng) fir 濾波器的輸入信號矢量 un8、traintestsample_2.m 將特征矩陣前 train_num 個為訓(xùn)練樣本，其余為測試樣本9、fir_nlms.dll nlms自適應(yīng)算法/downloads45/sourcecode/math/detail150408.html=recnstitution重構(gòu)相空間，在非線性系統(tǒng)分析中有重要的作用復(fù)制內(nèi)容到剪貼板代碼:function texp,lexp=lyapunov(n,tstart,stept,tend,ystart,io

5、utp);global ds;global p;global calculation_progress first_call;global driver_window;global trj_bufer time_bufer bufer_i;% lyapunov exponent calcullation for ode-system.% the alogrithm employed in this m-file for determining lyapunov% exponents was proposed in% a. wolf, j. b. swift, h. l. swinney, an

6、d j. a. vastano,% determining lyapunov exponents from a time series, physica d,% vol. 16, pp. 285-317, 1985.% for integrating ode system can be used any matlab ode-suite methods. % this function is a part of matds program - toolbox for dynamical system investigation% see: http:/www.math.rsu.ru/mexma

7、t/kvm/matds/% input parameters:% n - number of equation% rhs_ext_fcn - handle of function with right hand side of extended ode-system.% this function must include rhs of ode-system coupled with % variational equation (n items of linearized systems, see example). % fcn_integrator - handle of ode inte

8、grator function, for example: ode45 % tstart - start values of independent value (time t)% stept - step on t-variable for gram-schmidt renormalization procedure.% tend - finish value of time% ystart - start point of trajectory of ode system.% ioutp - step of print to matlab main window. ioutp=0 - no

9、 print, % if ioutp0 then each ioutp-th point will be print.% output parameters:% texp - time values% lexp - lyapunov exponents to each time value.% users have to write their own ode functions for their specified% systems and use handle of this function as rhs_ext_fcn - parameter. % example. lorenz s

10、ystem:% dx/dt = sigma*(y - x) = f1% dy/dt = r*x - y - x*z = f2% dz/dt = x*y - b*z = f3% the jacobian of system: % | -sigma sigma 0 |% j = | r-z -1 -x |% | y x -b |% then, the variational equation has a form:% % f = j*y% where y is a square matrix with the same dimension as j.% corresponding m-file:%

11、 function f=lorenz_ext(t,x)% sigma = 10; r = 28; beta = 8/3;% 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);% f=zeros(9,1);% f(1)=sigma*(y-x); f(2)=-x*z+r*x-y; f(3)=x*y-beta*z;% jac=-sigma,sigma,0; r-z,-1,-x; y, x,-beta;% % f(4:12)=jac*y;% run lyapunov exponen

12、t calculation:% % t,res=lyapunov(3,lorenz_ext,ode45,0,0.5,200,0 1 0,10); % % see files: lorenz_ext, run_lyap. % % -% copyright (c) 2004, govorukhin v.n.% this file is intended for use with matlab and was produced for matds-program% http:/www.math.rsu.ru/mexmat/kvm/matds/% lyapunov.m is free software

13、. lyapunov.m is distributed in the hope that it % will be useful, but without any warranty. % n=number of nonlinear odes% n2=n*(n+1)=total number of odes%options = odeset(reltol,ds(1).rel_error,abstol,ds(1).abs_error,maxstep,ds(1).max_step,.outputfcn,odeoutp,refine,0,initialstep,0.001);n_exp = ds(1)

14、.n_lyapunov;n1=n; n2=n1*(n_exp+1);neq=n2;% number of stepsnit = round(tend-tstart)/stept)+1;% memory allocation y=zeros(n2,1); cum=zeros(n2,1); y0=y;gsc=cum; znorm=cum;% initial valuesy(1:n)=ystart(:);for i=1:n_exp y(n1+1)*i)=1.0; end;t=tstart;fig_lyap = figure;set(fig_lyap,name,lyapunov exponents,n

15、umbertitle,off);set(fig_lyap,closerequestfcn,);hold on;box on;timeplot = tstart+(tend-tstart)/10;axis(tstart timeplot -1 1);title(dynamics of lyapunov exponents);xlabel(t);ylabel(lyapunov exponents);fig_lyap_axes = findobj(fig_lyap,type,axes);for i=1:n_expplotlyapi=plot(tstart,0); end;uu=findobj(fig

16、_lyap,type,line); for i=1:size(uu,1)set(uu(i),erasemode,none) ;set(uu(i),xdata,ydata,);set(uu(i),color,0 0 rand);enditerlyap = 0;% main loopcalculation_progress = 1;while ttend, tt = tend; end;% solutuion of extended ode system % t,y = feval(fcn_integrator,rhs_ext_fcn,t t+stept,y); while calculation

17、_progress = 1t,y = integrator(ds(1).method_int,ode_lin,t tt,y,options,p,n,neq,n_exp);first_call = 0;if calculation_progress = 99, break; end;if ( t(size(t,1)tt ) & (calculation_progress=0)y=y(size(y,1),:);y(1,1:n)=trj_bufer(bufer_i,1:n);t = time_bufer(bufer_i);calculation_progress = 1;elsecalcul

18、ation_progress = 0;end;end;if (calculation_progress = 99) break; elsecalculation_progress = 1;end;t=tt;y=y(size(y,1),:);first_call = 0;% construct new orthonormal basis by gram-schmidt%znorm(1)=0.0;for j=1:n1 znorm(1)=znorm(1)+y(n1+j)2; end;znorm(1)=sqrt(znorm(1);for j=1:n1 y(n1+j)=y(n1+j)/znorm(1);

19、 end;for j=2:n_expfor k=1:(j-1)gsc(k)=0.0;for l=1:n1 gsc(k)=gsc(k)+y(n1*j+l)*y(n1*k+l); end;end;for k=1:n1for l=1:(j-1)y(n1*j+k)=y(n1*j+k)-gsc(l)*y(n1*l+k);end;end;znorm(j)=0.0;for k=1:n1 znorm(j)=znorm(j)+y(n1*j+k)2; end;znorm(j)=sqrt(znorm(j);for k=1:n1 y(n1*j+k)=y(n1*j+k)/znorm(j); end;end;% upda

20、te running vector magnitudes%for k=1:n_exp cum(k)=cum(k)+log(znorm(k); end;% normalize exponent%rescale = 0;u1 =1.e10;u2 =-1.e10;for k=1:n_exp lp(k)=cum(k)/(t-tstart); % plot xd=get(plotlyapk,xdata);yd=get(plotlyapk,ydata);if timeplottu1=min(u1,min(yd);u2=max(u2,max(yd);end;xd=xd t; yd=yd lp(k);set(

21、plotlyapk,xdata,xd,ydata,yd);end;if timeplot0 then each ioutp-th point will be print.% output parameters:% texp - time values% lexp - lyapunov exponents to each time value.% users have to write their own ode functions for their specified% systems and use handle of this function as rhs_ext_fcn - para

22、meter. % example. lorenz system:% dx/dt = sigma*(y - x) = f1% dy/dt = r*x - y - x*z = f2% dz/dt = x*y - b*z = f3% the jacobian of system: % | -sigma sigma 0 |% j = | r-z -1 -x |% | y x -b |% then, the variational equation has a form:% % f = j*y% where y is a square matrix with the same dimension as

23、j.% corresponding m-file:% function f=lorenz_ext(t,x)% sigma = 10; r = 28; beta = 8/3;% 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);% f=zeros(9,1);% f(1)=sigma*(y-x); f(2)=-x*z+r*x-y; f(3)=x*y-beta*z;% jac=-sigma,sigma,0; r-z,-1,-x; y, x,-beta;% % f(4:12)=ja

24、c*y;% run lyapunov exponent calculation:% % t,res=lyapunov(3,lorenz_ext,ode45,0,0.5,200,0 1 0,10); % % see files: lorenz_ext, run_lyap. % % -% copyright (c) 2004, govorukhin v.n.% this file is intended for use with matlab and was produced for matds-program% http:/www.math.rsu.ru/mexmat/kvm/matds/% l

25、yapunov.m is free software. lyapunov.m is distributed in the hope that it % will be useful, but without any warranty. % n=number of nonlinear odes% n2=n*(n+1)=total number of odes%n1=n; n2=n1*(n1+1);% number of stepsnit = round(tend-tstart)/stept);% memory allocation y=zeros(n2,1); cum=zeros(n1,1);

26、y0=y;gsc=cum; znorm=cum;% initial valuesy(1:n)=ystart(:);for i=1:n1 y(n1+1)*i)=1.0; end;t=tstart;% main loopfor iterlyap=1:nit% solutuion of extended ode system t,y = unit(t,stept,y,d); t=t+stept;y=y(size(y,1),:);for i=1:n1 for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;end;% construct new orthonormal basis b

27、y gram-schmidt%znorm(1)=0.0;for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)2; end;znorm(1)=sqrt(znorm(1);for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;for j=2:n1for k=1:(j-1)gsc(k)=0.0;for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;end;for k=1:n1for l=1:(j-1)y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);end

28、;end;znorm(j)=0.0;for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)2; end;znorm(j)=sqrt(znorm(j);for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;end;% update running vector magnitudes%for k=1:n1 cum(k)=cum(k)+log(znorm(k); end;% normalize exponent%for k=1:n1 lp(k)=cum(k)/(t-tstart); end;% output modificationif

29、 iterlyap=1lexp=lp;texp=t;elselexp=lp;texp=t;end;for i=1:n1 for j=1:n1y(n1*j+i)=y0(n1*i+j);end;end;end;=小數(shù)據(jù)量法計算 lyapunov 指數(shù)的 matlab 程序 - (mex 函數(shù)，超快)-參考文獻(xiàn)：michael t.rosenstein,a practical method for calculating largest lyapunov exponents from small sets,physica d,1993,65: 117-134文件說明：-1、largestlyapun

30、ov_example1.m 程序主文件1，直接運(yùn)行此文件即可2、largestlyapunov_example2.m 程序主文件2，直接運(yùn)行此文件即可3、lorenzdata.dll 產(chǎn)生 lorenz 離散數(shù)據(jù)4、phasparecon.m 相空間重構(gòu)5、lyapunov_luzhenbo.dll lyapunov 計算主函數(shù)6、lyapunov_buffer.dll lyapunov 計算緩存/downloads63/sourcecode/math/detail221870.html=計算lyapunov指數(shù)的程序復(fù)制內(nèi)容到剪貼板代碼:program l

31、ylorenzparameter(n=3,m=12,st=100)integer:i,j,kreal y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3y(1)=10.y(2)=1.y(3)=0.a=10.b=8./3.r=28.t=0.h=0.01!initial conditionsdo i=n+1,my(i)=0.end dodo i=1,ny(n+1)*i)=1.s(i)=0end doopen(1,file=lorenz1.dat)open(2,file=ly lorenz.dat)do 100 k=1,st !st iterations

32、call rgkt(m,h,t,y,f,yc,y1)!normarize vector 123z=0. do i=1,ndo j=1,nz(i)=z(i)+y(n*j+i)*2enddoif(z(i)0.)z(i)=sqrt(z(i)do j=1,ny(n*j+i)=y(n*j+i)/z(i)enddoend do!generate gsr coefficientk1=0.k2=0.k3=0.do i=1,n k1=k1+y(3*i+1)*y(3*i+2)k2=k2+y(3*i+3)*y(3*i+2)k3=k3+y(3*i+1)*y(3*i+3)end do!conduct new vecto

33、r2 and 3do i=1,ny(3*i+2)=y(3*i+2)-k1*y(3*i+1) y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)end do!generate new vector2 and 3,normarize themdo i=2,nz(i)=0.do j=2,nz(i)=z(i)+y(n*j+i)*2enddoif(z(i)0.)z(i)=sqrt(z(i)do j=2,ny(n*j+i)=y(n*j+i)/z(i)end doend do!update lyapunov exponentdo i=1,nif(z(i)0)s(i)=s(i)

34、+log(z(i)enddo100 continuedo i=1,ns(i)=s(i)/(1.*st*h*1000.)write(2,*)s(i)enddoend!subroutine rgkt(m,h,t,y,f,yc,y1)real y(m),f(m),y1(m),yc(m),a,b,rinteger:i,jdo j=1,1000call df(m,t,y,f)t=t+h/2.0do i=1,myc(i)=y(i)+h*f(i)/2.0y1(i)=y(i)+h*f(i)/6.0 end do call df(m,t,yc,f)do i=1,myc(i)=y(i)+h*f(i)/2.0y1(

35、i)=y1(i)+h*f(i)/3.0 end do call df(m,t,yc,f)t=t+h/2.0do i=1,myc(i)=y(i)+h*f(i)y1(i)=y1(i)+h*f(i)/3.0end docall df(m,t,yc,f)do i=1,my(i)=y1(i)+h*f(i)/6.0end doif(j500)write(1,*)t,y(1),y(2),y(3)end doreturnend!subroutine df(m,t,y,f)real y(m),a,b,r,f(m)common a,b,ra=10.b=8./3.r=28.f(1)=a*(y(2)-y(1)f(2)

36、=y(1)*(r-y(3)-y(2)f(3)=y(1)*y(2)-b*y(3)do i=0,2f(4+i)=a*y(7+i)-y(4+i)f(7+i)=y(4+i)*(r-y(3)-y(7+i)-y(1)*y(10+i)f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)*y(7+i)enddoreturnend =c-c方法計算時間延遲和嵌入維數(shù)計算lyapunov指數(shù)計算關(guān)聯(lián)維數(shù)混沌時間序列預(yù)測c-c方法計算時間延遲和嵌入維數(shù)主程序：c_cmethod.m, c_cmethod_independent.m子函數(shù)：correlation_integral.m(計算關(guān)聯(lián)積分)d

37、isjoint.m(將時間序列拆分成t個不相關(guān)的子序列)heaviside.m（計算時間序列的海維賽函數(shù)值）參考文獻(xiàn)nonlinear dynamics, delay times, and embedding windows。計算lyapunov指數(shù)：largest_lyapunov_exponent.m（用呂金虎的方法計算最大lyapunov指數(shù)）參考文獻(xiàn)：基于lyapunov指數(shù)改進(jìn)算法的邊坡位移預(yù)測。lyapunov_wolf.m(用wolf方法計算最大lyapunov指數(shù))計算關(guān)聯(lián)維數(shù)：g_p.m(g-p算法)混沌時間序列預(yù)測主函數(shù)mainpre_by_jiaquanyijie_1.m

38、(該程序用加權(quán)一階局域法對數(shù)據(jù)進(jìn)行進(jìn)行一步預(yù)測) mainpre_by_jiaquanyijie_n.m(該程序用加權(quán)一階局域法對數(shù)據(jù)進(jìn)行進(jìn)行n步預(yù)測)mainpre_by_lya_1.m(基于最大lyapunov指數(shù)的一步預(yù)測)mainpre_by_lya_n.m(基于最大lyapunov指數(shù)的n步預(yù)測)nearest_point.m(計算最后一個相點(diǎn)的最近相點(diǎn)的位置及最短距離)子函數(shù)jiaquanyijie.m(該函數(shù)用加權(quán)一階局域法(xx)、零級近似(yy)和基于零級近似的加權(quán)一階局域法(zz)對時間數(shù)據(jù)進(jìn)行一步預(yù)測)pre_by_lya.m(基于最大lyapunov指數(shù)的預(yù)測方法)pr

39、e_by_lya_new.m(改進(jìn)的基于最大lyapunov指數(shù)的預(yù)測方法)如果需要該工具包，大家到/downloads25/ . rs/detail82624.html下載=求取lyapunov指數(shù)的小數(shù)據(jù)量方法，采用混合編程lylorenz.f90的程序如下：program lylorenzparameter(n=3,m=12,st=100)integer:i,j,kreal y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3y(1)=10.y(2)=1.y(3)=0.a=10.b=8./3.r=28.t=0

40、.h=0.01!initial conditionsdo i=n+1,my(i)=0.end dodo i=1,ny(n+1)*i)=1.s(i)=0end doopen(1,file=lorenz1.dat)open(2,file=ly lorenz.dat)do 100 k=1,st !st iterationscall rgkt(m,h,t,y,f,yc,y1)!normarize vector 123z=0. do i=1,ndo j=1,nz(i)=z(i)+y(n*j+i)*2enddoif(z(i)0.)z(i)=sqrt(z(i)do j=1,ny(n*j+i)=y(n*j+i

41、)/z(i)enddoend do!generate gsr coefficientk1=0.k2=0.k3=0.do i=1,n k1=k1+y(3*i+1)*y(3*i+2)k2=k2+y(3*i+3)*y(3*i+2)k3=k3+y(3*i+1)*y(3*i+3)end do!conduct new vector2 and 3do i=1,ny(3*i+2)=y(3*i+2)-k1*y(3*i+1) y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)end do!generate new vector2 and 3,normarize themdo i=2

42、,nz(i)=0.do j=2,nz(i)=z(i)+y(n*j+i)*2enddoif(z(i)0.)z(i)=sqrt(z(i)do j=2,ny(n*j+i)=y(n*j+i)/z(i)end doend do!update lyapunov exponentdo i=1,nif(z(i)0)s(i)=s(i)+log(z(i)enddo100 continuedo i=1,ns(i)=s(i)/(1.*st*h*1000.)write(2,*)s(i)enddoend!subroutine rgkt(m,h,t,y,f,yc,y1)real y(m),f(m),y1(m),yc(m),

43、a,b,rinteger:i,jdo j=1,1000call df(m,t,y,f)t=t+h/2.0do i=1,myc(i)=y(i)+h*f(i)/2.0y1(i)=y(i)+h*f(i)/6.0 end do call df(m,t,yc,f)do i=1,myc(i)=y(i)+h*f(i)/2.0y1(i)=y1(i)+h*f(i)/3.0 end do call df(m,t,yc,f)t=t+h/2.0do i=1,myc(i)=y(i)+h*f(i)y1(i)=y1(i)+h*f(i)/3.0end docall df(m,t,yc,f)do i=1,my(i)=y1(i)

44、+h*f(i)/6.0end doif(j500)write(1,*)t,y(1),y(2),y(3)end doreturnend!subroutine df(m,t,y,f)real y(m),a,b,r,f(m)common a,b,ra=10.b=8./3.r=28.f(1)=a*(y(2)-y(1)f(2)=y(1)*(r-y(3)-y(2)f(3)=y(1)*y(2)-b*y(3)do i=0,2f(4+i)=a*y(7+i)-y(4+i)f(7+i)=y(4+i)*(r-y(3)-y(7+i)-y(1)*y(10+i)f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)