《matlab编写的Lyapunov指数计算程序.doc》由会员分享，可在线阅读，更多相关《matlab编写的Lyapunov指数计算程序.doc（22页珍藏版）》请在金锄头文库上搜索。
1、matlab编写的Lyapunov指数计算程序matlab编写的Lyapunov指数计算程序 已有 2406 次阅读 2009-12-29 08:37 |个人分类:其它|系统分类:科普集锦|关键词:李雅普诺夫指数 一、计算连续方程Lyapunov指数的程序其中给出了两个例子：计算Rossler吸引子的Lyapunov指数 计算超混沌Rossler吸引子的Lyapunov指数http:/ Texp,Lexp=lyapunov(n,tstart,stept,tend,ystart,ioutp);global DS;global P;global calculation_progress first
2、_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,% Determining Lyapunov Exponents from a T
3、ime 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/mexmat/kvm/matds/% Input parameters:% n - number of equation%
4、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 valu
5、es 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 ioutp0 then each ioutp-th point will be prin
6、t.% 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
7、 = 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
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
9、,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
10、 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 st
11、epsnit = 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,NumberTitle,off);set(Fig_Lyap,CloseRequestFcn,);hold on;
12、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 PlotLyapi=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); 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);
