- function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
- %
- % 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 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/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 所调用求微分方程的MATLAB函数名
- % tstart - start values of independent value (time t)时间的初始值
- % stept - step on t-variable for Gram-Schmidt renormalization
- % procedure.进行Gram_Schmidt正交化的迭代次数
- % 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.每隔ioutp个点输出相应的时间和李雅普诺夫指数数值
- %
- % 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.
- %用户如果要计算其他方程，需要写出自己的常微分方程组来替换lorenz_ext文件内容得到相应的函数表达式。
- % 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] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);
-
- 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=[Lexp; lp];
- Texp=[Texp; t];
- end;
- if (mod(ITERLYAP,ioutp)==0)
- fprintf('t=%6.4f',t);
- for k=1:n1 fprintf(' %10.6f',lp(k)); end;
- fprintf('\n');
- end;
- for i=1:n1
- for j=1:n1
- y(n1*j+i)=y0(n1*i+j);
- end;
- end;
- end;
|
发表于 2014-7-1 11:28
显身卡
楼主