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,
- % "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 = 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<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);
- 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))<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
- 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<t
- 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<t
- 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;
-
|
[复制链接]
发表于 2006-9-2 06:07
提升卡
变色卡
显身卡
楼主
发表于 2006-10-7 18:22
