matlab编写的Lyapunov指数计算程序

matlab编写的Lyapunov指数计算程序

已有 2406 次阅读 2009-12-29 08:37 |个人分类:其它|系统分类:科普集锦|关键词:李雅普诺夫指数

一、计算连续方程Lyapunov指数的程序 其中给出了两个例子：

计算Rossler吸引子的Lyapunov指数

计算超混沌Rossler吸引子的Lyapunov指数

http://www.pudn.com/downloads39/sourcecode/math/detail132231.html

二、recnstitution重构相空间，在非线性系统分析中有重要的作用 function [Texp,Lexp]=lyapunov(n,tstart,stept,tend,ystart,ioutp); 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, and J. A. Vastano,

% \% 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/mexmat/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 integrator 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 print, % if ioutp>0 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 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 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)=Jac*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/

% lyapunov.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 %

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).n_lyapunov; n1=n; n2=n1*(n_exp+1); neq=n2;

% Number of steps

nit = round((tend-tstart)/stept)+1;

% Memory allocation

y=zeros(n2,1); cum=zeros(n2,1); y0=y;

gsc=cum; znorm=cum;

% Initial values

y(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','NumberTitle','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_exp

PlotLyap{i}=plot(tstart,0); end;

uu=findobj(Fig_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]); end

ITERLYAP = 0;

% Main loop

calculation_progress = 1;

while t

tt = t + stept;

ITERLYAP = ITERLYAP+1; if tt>tend, tt = tend; end;

% Solutuion of extended ODE system

% [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y); while calculation_progress == 1

[T,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))

y(1,1:n)=TRJ_bufer(bufer_i,1:n); t = Time_bufer(bufer_i); calculation_progress = 1; else

calculation_progress = 0; end; end;

if (calculation_progress == 99) break; else

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

for j=2:n_exp for 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:n1

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

%

% update 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(PlotLyap{k},'Xdata'); Yd=get(PlotLyap{k},'Ydata');

if timeplot

u1=min(u1,min(Yd)); u2=max(u2,max(Yd)); end;

Xd=[Xd t]; Yd=[Yd lp(k)];

set(PlotLyap{k},'Xdata',Xd,'Ydata',Yd); end;

if timeplot

timeplot = timeplot+(tend-tstart)/20; figure(Fig_Lyap);

axis([tstart timeplot u1 u2]); end;

drawnow;

% Output modification

if ITERLYAP==1 Lexp=lp; Texp=t;

else

Lexp=[Lexp; lp]; Texp=[Texp; t]; end;

if (mod(ITERLYAP,ioutp)==0) for k=1:n_exp

txtstring{k}=['\\lambda_' int2str(k) '=' num2str(lp(k))]; end

legend(Fig_Lyap_Axes,txtstring,3); end; end;

ss=warndlg('Attention! Plot of lyapunov exponents will be closed!','Press OK to continue!'); uiwait(ss);

delete(Fig_Lyap);

fprintf('\\n \\n Results of Lyapunov exponents calculation: \\n t=%6.4f',t); for k=1:n_exp fprintf(' L%d=%f; ',k,lp(k)); end; fprintf('\\n');

if ~isempty(driver_window) if ishandle(driver_window) delete(driver_window); driver_window = []; end; end;

calculation_progress = 0;

update_ds; 三、wolf 方法计算李雅普诺夫指数 http://www.pudn.com/downloads52/sourcecode/math/detail178303.html

四、给出了分形计算的源代码的matlab程序，可以迅速帮助大家进行分形的分析与计算，参数容易设置

function [Texp,Lexp]=new1lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp,d); %

% 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, and J. A. Vastano,

% \% 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/mexmat/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 integrator 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 print, % if ioutp>0 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 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 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)=Jac*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/

% lyapunov.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 steps

nit = round((tend-tstart)/stept);

% Memory allocation

y=zeros(n2,1); cum=zeros(n1,1); y0=y; gsc=cum; znorm=cum;

% Initial values

y(1:n)=ystart(:);

for i=1:n1 y((n1+1)*i)=1.0; end;

t=tstart;

% Main loop

for 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 by 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:n1

for 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:n1

for l=1:(j-1)

y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l); end; 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 modification

if ITERLYAP==1 Lexp=lp; Texp=t; else

Lexp=lp; Texp=t; end;

for i=1:n1 for j=1:n1

y(n1*j+i)=y0(n1*i+j); end;

end; end;

五、小数据量法计算 Lyapunov 指数的 Matlab 程序 - (mex 函数，超快) http://www.pudn.com/downloads63/sourcecode/math/detail221870.html

六、计算lyapunov指数的程序 program lylorenz

parameter(n=3,m=12,st=100) integer::i,j,k

real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3 y(1)=10. y(2)=1. y(3)=0. a=10. b=8./3. r=28. t=0. h=0.01

!!!!!initial conditions do i=n+1,m

y(i)=0. end do do i=1,n

y((n+1)*i)=1. s(i)=0 end do

open(1,file='lorenz1.dat') open(2,file='ly lorenz.dat')

do 100 k=1,st !!!!!!!!st iterations call rgkt(m,h,t,y,f,yc,y1)

!!!!normarize vector 123 z=0. do i=1,n

do j=1,n

z(i)=z(i)+y(n*j+i)**2 enddo

if(z(i)>0.)z(i)=sqrt(z(i)) do j=1,n

y(n*j+i)=y(n*j+i)/z(i)

enddo end do

!!!!generate gsr coefficient k1=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 3 do i=1,n

y(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 them do i=2,n z(i)=0. do j=2,n

z(i)=z(i)+y(n*j+i)**2 enddo

if(z(i)>0.)z(i)=sqrt(z(i)) do j=2,n

y(n*j+i)=y(n*j+i)/z(i) end do end do

!!!!!!!update lyapunov exponent do i=1,n

if(z(i)>0)s(i)=s(i)+log(z(i)) enddo

100 continue do i=1,n

s(i)=s(i)/(1.*st*h*1000.) write(2,*)s(i) enddo end

!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine rgkt(m,h,t,y,f,yc,y1)

real y(m),f(m),y1(m),yc(m),a,b,r integer::i,j do j=1,1000

call df(m,t,y,f) t=t+h/2.0 do i=1,m

yc(i)=y(i)+h*f(i)/2.0

y1(i)=y(i)+h*f(i)/6.0 end do call df(m,t,yc,f) do i=1,m

yc(i)=y(i)+h*f(i)/2.0

y1(i)=y1(i)+h*f(i)/3.0 end do call df(m,t,yc,f) t=t+h/2.0 do i=1,m

yc(i)=y(i)+h*f(i)

y1(i)=y1(i)+h*f(i)/3.0 end do

call df(m,t,yc,f) do i=1,m

y(i)=y1(i)+h*f(i)/6.0 end do

if(j>500)write(1,*)t,y(1),y(2),y(3) end do return end

!!!!!!!!!!!!!!!!!!!!!!!! subroutine df(m,t,y,f) real y(m),a,b,r,f(m) common a,b,r a=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,2

f(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) enddo

return end

七、计算各种混沌系统李雅普洛夫指数的MATLAB 源程序 http://www.pudn.com/downloads50/sourcecode/math/detail172856.html