Lyapunov指数计算的ode系统-MATLAB

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
% 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 = rx - y - xz = f2
% dz/dt = xy - bz = f3
%
% The Jacobian of system:
% | -sigma sigma 0 |
% J = | r-z -1 -x |
% | y x -b |
%
% Then, the variational equation has a form:
%
% F = JY
% 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)=-xz+Rx-y; f(3)=xy-BETAz;
%
% Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
%
% f(4:12)=JacY;
%
% 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 © 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);%步数

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(n1i+j)=y(n1j+i); end;% 4 7 10 5 8 11 6 9 12
end; % 4 5 6 7 8 9 10 11 12

%
% 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(n1j+1)=y0(n1j+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(n1l+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(n1j+i)=y0(n1i+j);
end;
end;

end;