matlab编写的Lyapunov指数计算程序

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;
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;
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');
% 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)
% Lexp - Lyapunov exponents to each time value.
%
% Users have to write their own ODE functions for their specified
% systems and use handle of this fmeter.
while t<tend
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);
% 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');
% 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)];
matlab编写的lyapunov指数计算程序文档精品是我多年收藏整理的好文档都是文档中的经典精品极品确实值得下载收藏
matlab编写的Lyapunov指数计算程序
已有2406次阅读2009-12-29 08:37 |个人分类:其它|系统分类:科普集锦|关键词:李雅普诺夫指数
一、计算连续方程Lyapunov指数的程序
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;
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;
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;
else
% This function is a part of MATDS program - toolbox for dynamical system investigation
% See: http://www.math.rsu.ru/mexmat/kvm/matds/
%
% Input parameters:
znorm(j)=sqrt(znorm(j));
for k=1:n1 y(n1*j+k)=y(n1*j+k)/znorm(j); end;
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
% "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.
%
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;
%
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;
% 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
%
% [T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);
%
% See files: lorenz_ext, run_lyap.
%
% --------------------------------------------------------------------
%
% 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,
% 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
% 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
set(Fig_Lyap,'CloseRequestFcn','');
hold on;
box on;
timeplot = tstart+(tend-tstart)/10;
axis([tstart timeplot -1 1]);
title('Dynamics of Lyapunov exponents');
% 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:
% 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/
% 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:
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;
%
% 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 |
其中给出了两个例子：
计算Rossler吸引子的Lyapunov指数
计算超混沌Rossler吸引子的Lyapunov指数
/downloads39/sourcecode/math/detail132231.html
二、recnstitution重构相空间，在非线性系统分析中有重要的作用