findlyap.m

计算各种混沌系统李雅普洛夫指数的MATLAB 源程序。

function Line=findlyap(MainHandles)
%FINDLYAP  determines the Lyapunov exponents and dimension
%
%    The alogrithm employed in this toolbox for determining Lyapunov
%    exponents is according to the algorithms proposed in
%
%    [1] 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.
%
%    [2] J. P. Eckmann and D. Ruelle, "Ergodic Theory of Chaos and Strange
%        Attractors," Rev. Mod. Phys., Vol. 57, pp. 617-656, 1985.
%
%    The algorithm given in [1] is used for first-order systems while
%    the QR-based algorithm proposed in [2] is applied for higher order
%    systems.

%   by Steve W. K. SIU, July 5, 1998.


%Print the results to the file every iteration
PrintStep=10;
%Display the results on the screen every iteration
DisplayStep=1;

%Clear the current axis objects
cla;
v=axis;
xc=(v(2)+v(1))/2;	%Center of the axis box [xc,y]
y=(v(4)+v(3))/2;
aW=v(2)-v(1);		%Axis width
x=xc-1/4*aW;

%Display the text "Loading... Please wait!"
th=text(x,y,'Loading... Please wait!','Color','r','FontSize',15,...
   	'FontAngle','italic');
drawnow;
%Clear the bottom text if any
set(MainHandles(1),'String','');

%Get the data stored in 'UserData' of the setting button
DATA=get(MainHandles(4),'UserData');
%Restore the data input by the user
output=DATA(1);		%Output checkbox: "on"=1; "off"=0
LEout=DATA(2);			%Output Lyapunov exponents to file: "yes"=1; "no"=0
LEprec=DATA(3);      %Precision of the output Lyapunov exponents
LDout=DATA(4);			%Output Lyapunov dimension to file: "yes"=1; "no"=0
LDprec=DATA(5);      %Precision of the Lyapunov dimension
IntMethod=DATA(6);	%Integration method: 1=Discrete map, 2=ODE45, 3=ODE23
							% 4=ODE113, 5=ODE23S, 6=ODE15S
InitialTime=DATA(7);	%Initial time
FinalTime=DATA(8);	%Final time
TimeStep=DATA(9);		%Time step
RelTol=DATA(10);		%Relative tolerance
AbsTol=DATA(11);		%Absolute tolerance
plot1=DATA(12);		%Plot imediately: 1="checked", 0="unchecked"
plot2=DATA(13);		%Plot according to specified iterations
ItrNum=DATA(14);		%No. of iteration for updating the plot
Color=DATA(15);		%Line Color option
							%Line color: 1="blue", 2="balck", 3="green", 4="red",
                     %5="yellow", 6="magenta", 7="cyan"
DiscardItr=DATA(16);	%Iterations to be discarded
UpdateStepNum=DATA(17);	%Lyapunov exponents updating steps
linODEnum=DATA(18);	%No. of linearized ODEs
ic=DATA(19:length(DATA));	%Initial conditions

%Get the output file and ODE function names stored in
%'UserData' of the "Start" button.
%Handles(2) is the handle of "Start" button
NAMES=get(MainHandles(2),'UserData');
OutputFile=rmspace(NAMES(1,:));
odefun=rmspace(NAMES(2,:));

%Construct a look-up table for the line colors
COLORS='bkgrymc';
%Map the "Line color" pop-up menu position to the look-up table
LineColor=COLORS(Color);	

%Construct a look-up table for integration methods
Methods=char('Discrete map', 'ode45','ode23','ode113','ode23s','ode15s');
ODEsolver=strcat(Methods(IntMethod,:));

%Dimension of the linearized system (total: d x d ODEs)
d=sqrt(linODEnum);
%Initial conditions for the linearized ODEs
Q0=eye(d);
IC=[ic(:);Q0(:)];
ICnum=length(IC);		%Total no. of initial coniditions
%One iteration: Duration for updating the LEs
Iteration=UpdateStepNum*TimeStep;	
DiscardTime=DiscardItr*Iteration+InitialTime;

%MATLAB's ODE functions will give the intermediate solutions if 
%the duration between the initial time and the final time is only
%one time step, this will slow down the whole iteration process. 
%To avoid this, reduce the time step by half.
if (UpdateStepNum==1 & IntMethod~=3)
   TimeStep=TimeStep/2;
end

T1=InitialTime;
T2=T1+Iteration;
TSpan=[T1:TimeStep:T2];
%Absolute tolerance of each components is set to the same value
options=odeset('RelTol',RelTol,'AbsTol',ones(1,ICnum)*AbsTol);

%Initialize variables
n=0;			%Iteration counter
k=0;			%Effective iteration counter
				% (discarded iterations are not counted)
h=1;			%No. of line handles sets (1 set = d line handles)
delLine=0;  %Indicator for deleting the drawn lines            
Sum=zeros(1,d);
xData=[];
yData=[];
Line=[];
bufferSize=10000; %Max. no. of data can be stored in the buffer before creating a new line.
if ( output & (LEout | LDout) )
   %If the output file cannot be opened, warn the user
   msg=['Unable to open "' sprintf(OutputFile) '".'];
   Warn='errordlg(msg,''ERROR'',''replace''); problem=1;';
   eval('fid=fopen(OutputFile,''wt''); problem=0;',Warn)
   %Construct a look-up table for precision format
   Prec=char('%.4f','%.6f','%.8f','%.10f','%.12f');
   %Map the pop-up menu position to its corresponding format
   LEprecision=strcat(Prec(LEprec,:));
   LDprecision=strcat(Prec(LDprec,:));
   if ~problem
      fprintf(fid,'Time');
      if LEout
         for i=1:d
            Str1=['\tLE%d'];
            fprintf(fid,Str1,i);
         end
      end
      if LDout
         Str2=['\tLD'];
         fprintf(fid,Str2); 
      end
      fprintf(fid,'\n');
   end
else
   problem=0;
end

%Start the stop watch
tic;
%Get the state of the "Stop" button
stop=get(MainHandles(3),'UserData');
if DiscardTime>0
   %Display the text "Discarding transient steps..."
   set(MainHandles(1),'String','Discarding transient steps...');
end

A=[];

%String that contains the integration command
IntegrationStr=['[t,X]=',ODEsolver,'(odefun,TSpan,IC,options);'];
%Main loop
while (~stop & ~problem)
   n=n+1;
   %Integration
   if IntMethod>1
      eval(IntegrationStr);
   else	%If it is a discrete map
      for i=1:UpdateStepNum
         X(i,:)=(feval(odefun,IC))';
      end
   end
   [rX,cX]=size(X);
   L=cX-linODEnum;      %No. of initial conditions for 
                        %the original system
   for i=1:d
      m1=L+1+(i-1)*d;
      m2=m1+d-1;
      A(:,i)=(X(rX,m1:m2))';
   end
   %QR decomposition
   if d>1
      %The algorithm for 1st-order system doesn't require
      %QR decomposition
      [Q,R]=qr(A);
      if T2>DiscardTime
         Q0=Q;
      else
         Q0=eye(d);
      end
   else
      R=A;
   end
      
   
   %Delete the text "Loading...Please wait!" 
   %before the first iteration
   if n==1
      delete(th);
      drawnow;
      %Display the final time
      set(MainHandles(11),'String',FinalTime);
   end
   
  %Any zero diagonal element will cause overflow
  %in the following calculation, so discard this step.
   permission=1;
   for i=1:d
      if R(i,i)==0
         permission=0;
         break;
      end
   end
  %To determine the Lyapunov exponents
   if (T2>DiscardTime & permission)
      k=k+1;
      T=k*Iteration;
      TT=n*Iteration+InitialTime;
      %There are d Lyapunov exponents
      Sum=Sum+log(abs(diag(R))');
      lambda=Sum/T;
      
      %Sort the Lyapunov exponents in descenting order
      Lambda=fliplr(sort(lambda));
      %To calculate the Lyapunov dimension (or Kaplan-Yorke dimension)
      LESum=Lambda(1);			
      LD=0;
      if (d>1 & Lambda(1)>0)
         for N=1:d-1
            if Lambda(N+1)~=0
               LD=N+LESum/abs(Lambda(N+1));
               LESum=LESum+Lambda(N+1);
               if LESum<0
                  break;
               end
            end
         end
      end
      %Store the [x,y] data for plotting
      [rxD,cxD]=size(xData);
      [ryD,cyD]=size(yData);
      if rxD<=bufferSize
         xData=[xData;TT];
         yData=[yData;lambda];
      else
         %When the buffers are full, refresh them
         %Max. size of buffers = 10000 data
         xData=[xData(rxD);TT];
         yData=[yData(ryD,:);lambda];
         h=h+1;		%add one set of line handles
         delLine=0;	%After refreshing the buffers, the
                     % previous drawn lines must not be deleted
      end
      
      if ( output & ~problem & (LEout | LDout) & rem(k,PrintStep)==0)
         fprintf(fid,'%.2f',TT);
         if LEout
            for i=1:d
               Str1=['\t',LEprecision];
               fprintf(fid,Str1,Lambda(i));
            end
         end
         if LDout
            Str2=['\t',LDprecision];
            fprintf(fid,Str2,LD);
         end
         fprintf(fid,'\n');
      end

      %Draw lines immediately if "update the plot immediately" was chosen.
      if (plot1==1 | plot2==1 & rem(k,ItrNum)==0)
         if delLine
            %Clear the previous drawn line if any
            delete(Line(h,:));
         end
      	%Draw d lines      
         for i=1:d
            %Set "Erase Mode" to "none" for increasing speed (less refresh)
            Line(h,i)=line('EraseMode','none','Color',LineColor,...
                      'xData',xData,'yData',yData(:,i));
            %Force MATLAB to draw immediately
            drawnow;
         end
         delLine=1;	%Set a flag to indicate that the lines now can be deleted
      end
      %Display the calculated Lyapunov exponents
      if rem(k,DisplayStep)==0
         set(MainHandles(1),'String',[num2str(Lambda),blanks(3),'( ',num2str(LD),' )']);
      end
   end
   

   %To see whether "Stop" is pressed
   stop=get(MainHandles(3),'UserData');
   
   %Display current time, and time used
   set(MainHandles(10),'String',num2str(round(T2)));		%Current time
   set(MainHandles(12),'String',num2str(round(toc)));		%Used time
   drawnow;
   
   %If calculation is finished or "stop" button is pressed, exit the loop.
   if (stop | T2>=FinalTime)
      %Reset the "Erase mode" to normal
      set(Line,'EraseMode','normal');
      %Show the final results (for making sure the final results being shown if DisplayStep>1)
      if (T2>DiscardTime & permission)
         set(MainHandles(1),'String',[num2str(Lambda),blanks(3),'( ',num2str(LD),' )']);
      end
      break;
   end
   %Update the initial conditions and time span for the next iteration
   if IntMethod>1
      ic=X(rX,1:L);
      T1=T1+Iteration;
      T2=T2+Iteration;
      TSpan=[T1:TimeStep:T2];
   else %For discrete map
      ic=X(UpdateStepNum,1:L);
      T2=T2+Iteration;
   end
   IC=[ic(:);Q0(:)];
end		%End of main loop

if T2>=FinalTime
   set(MainHandles([2,4,13,15]),'Enable','On');
   set(MainHandles(3),'Enable','Off');
end
if (output & ~problem)
   fclose(fid);
end

%----------------Subroutine-----------------------------------
function outStr=rmspace(inStr)
%RMSPACE		Function for removing the beginning and ending
%				spaces of a string

%Remove spaces at the end of the string
outStr=strcat(inStr);
%Delete spaces at the beginning of the string
if ~isempty(outStr)
   while isspace(outStr(1))
   	outStr=outStr(2:length(outStr));
   end
end