www.pudn.com > NumericalComputingwithMatlabCode.zip > lorenzgui.m, change:2003-10-27,size:7117b


function lorenzgui 
%LORENZGUI   Plot the orbit around the Lorenz chaotic attractor. 
%   This function animates the integration of the three coupled 
%   nonlinear differential equations that define the Lorenz Attractor, 
%   a chaotic system first described by Edward Lorenz of MIT. 
%   As the integration proceeds you will see a point moving in 
%   an orbit in 3-D space known as a strange attractor. 
%   The orbit ranges around two different critical points, or attractors. 
%   The orbit is bounded, but may not be periodic and or convergent. 
% 
%   The mouse and arrow keys change the 3-D viewpoint.  Uicontrols 
%   provide "pause", "resume", "stop", "restart", "clear", and "close". 
% 
%   A listbox provides a choice among five values of the parameter rho. 
%   The first value, 28, is the most common and produces the chaotic 
%   behavior.  The other four values values produce periodic behaviors 
%   of different complexities.  A change in rho becomes effective only 
%   after a "stop" and "restart". 
% 
%   Reference: Colin Sparrow, "The Lorenz Equations: Bifurcations, 
%   Chaos, and Strange Attractors", Springer-Verlag, 1982. 
 
if ~isequal(get(gcf,'name'),'Lorenz Gui') 
    
   % This is first entry, just initialize the figure window. 
 
   rhos = [28 99.65 100.5 160 350]; 
   shg 
   clf reset 
   p = get(gcf,'pos'); 
   set(gcf,'color','black','doublebuff','on','name','Lorenz Gui', ... 
      'menu','none','numbertitle','off', ... 
      'pos',[p(1) p(2)-(p(3)-p(4))/2 p(3) p(3)]) 
 
   % Callback to erase comet by jiggling figure position 
 
   klear = ['set(gcf,''pos'',get(gcf,''pos'')+[0 0 0 1]), drawnow,' ... 
            'set(gcf,''pos'',get(gcf,''pos'')-[0 0 0 1]), drawnow']; 
 
   % Uicontrols 
 
   paws = uicontrol('style','toggle','string','start', ... 
      'units','norm','pos',[.02 .02 .10 .04],'value',0, ... 
      'callback','lorenzgui'); 
   stop = uicontrol('style','toggle','string','close', ... 
      'units','norm','pos',[.14 .02 .10 .04],'value',0, ... 
      'callback','cameratoolbar(''close''), close(gcf)'); 
   clear = uicontrol('style','push','string','clear', ... 
      'units','norm','pos',[.26 .02 .10 .04], ... 
      'callback',klear); 
   rhostr = sprintf('%6.2f|',rhos); 
   rhopick = uicontrol('style','listbox','tag','rhopick', ... 
      'units','norm','pos',[.82 .02 .14 .14], ... 
      'string',rhostr(1:end-1),'userdata',rhos,'value',1); 
 
else 
 
   % The differential equation is ydot = A(y)*y 
   % With this value of eta, A is singular. 
   % The eta's in A will be replaced by y(2) during the integration. 
 
   rhopick = findobj('tag','rhopick'); 
   rhos = get(rhopick,'userdata'); 
   rho = rhos(get(rhopick,'value')); 
   sigma = 10; 
   beta = 8/3; 
   eta = sqrt(beta*(rho-1)); 
   A = [ -beta    0     eta 
            0  -sigma   sigma  
         -eta   rho    -1  ]; 
    
   % The critical points are the null vectors of A. 
   % The initial value of y(t) is near one of the critical points. 
    
   yc = [rho-1; eta; eta]; 
   y0 = yc + [0; 0; 3]; 
    
   % Integrate forever, or until the stop button is toggled. 
    
   tspan = [0 Inf]; 
   opts = odeset('reltol',1.e-6,'outputfcn',@lorenzplot,'refine',4); 
   ode45(@lorenzeqn, tspan, y0, opts, A); 
 
end 
 
 
% ------------------------------ 
 
function ydot = lorenzeqn(t,y,A) 
%LORENZEQN  Equation of the Lorenz chaotic attractor. 
%   ydot = lorenzeqn(t,y,A). 
%   The differential equation is written in almost linear form. 
%      ydot = A*y 
%   where 
%      A = [ -beta    0     y(2) 
%               0  -sigma   sigma  
%            -y(2)   rho    -1  ]; 
 
A(1,3) = y(2); 
A(3,1) = -y(2); 
ydot = A*y; 
 
 
% ------------------------------ 
 
function fin = lorenzplot(t,y,job,A) 
%LORENZPLOT   Plot the orbit of the Lorenz chaotic attractor. 
 
persistent Y 
 
if isequal(job,'init') 
 
   % Initialize axis and comet, R = axis settings, L = length of comet. 
 
   rho = A(3,2); 
   switch rho 
      case 28,    R = [  5  45  -20  20  -25  25];  L = 100; 
      case 99.65, R = [ 50 150  -35  35  -60  60];  L = 240; 
      case 100.5, R = [ 50 150  -35  35  -60  60];  L = 120; 
      case 160,   R = [100 220  -40  40  -75  75];  L = 165; 
      case 350,   R = [285 435  -55  55 -105 105];  L =  80; 
      otherwise,  R = [100 250  -50  50 -100 100];  L = 150; 
   end 
   set(gcf,'pos',get(gcf,'pos')+[0 0 0 1]) 
   drawnow 
   set(gcf,'pos',get(gcf,'pos')-[0 0 0 1]) 
   drawnow 
   if get(gca,'userdata') ~= rho, delete(gca), end 
   set(gca,'color','black','pos',[.03 .05 .93 .95],'userdata',rho) 
   axis(R); 
   axis off 
 
   comet(1) = line(y(1),y(2),y(3),'linestyle','none','marker','.', ... 
      'erasemode','xor','markersize',25); 
   comet(2) = line(NaN,NaN,NaN,'color','y','erasemode','none'); 
   comet(3) = line(NaN,NaN,NaN,'color','y','erasemode','none'); 
   Y = y(:,ones(L,1)); 
 
   uics = flipud(get(gcf,'children')); 
   paws = uics(1); 
   stop = uics(2); 
   set(paws,'string','pause','callback','','value',0); 
   set(stop,'string','stop','callback','','value',0); 
 
   beta = -A(1,1); 
   eta = sqrt(beta*(rho-1)); 
   yc = [rho-1; eta; eta]; 
   line(yc(1),yc(2),yc(3),'linestyle','none','marker','o','color','g') 
   line(yc(1),-yc(2),-yc(3),'linestyle','none','marker','o','color','g') 
 
   ax = [R(2) R(1) R(1) R(1) R(1)]; 
   ay = [R(3) R(3) R(4) R(3) R(3)]; 
   az = [R(5) R(5) R(5) R(5) R(6)]; 
   p = .9; 
   q = 1-p; 
   grey = [.4 .4 .4]; 
   line(ax,ay,az,'color',grey); 
   text(p*R(1)+q*R(2),R(3),p*R(5),sprintf('%3.0f',R(1)),'color',grey) 
   text(q*R(1)+p*R(2),R(3),p*R(5),sprintf('%3.0f',R(2)),'color',grey) 
   text(R(1),p*R(3)+q*R(4),p*R(5),sprintf('%3.0f',R(3)),'color',grey) 
   text(R(1),q*R(3)+p*R(4),p*R(5),sprintf('%3.0f',R(4)),'color',grey) 
   text(R(1),R(3),p*R(5)+q*R(6),sprintf('%3.0f',R(5)),'color',grey) 
   text(R(1),R(3),q*R(5)+p*R(6),sprintf('%3.0f',R(6)),'color',grey) 
   fin = 0; 
 
   cameratoolbar('setmode','orbit') 
   uicontrol('style','text','units','norm','pos',[.38 .02 .34 .04], ... 
      'foreground','white','background','black','fontangle','italic', ... 
      'string','Click on axis to rotate view') 
 
elseif isequal(job,'done') 
 
   fin = 1; 
 
else 
 
   % Update comet 
 
   L = size(y,2); 
   Y(:,end+1:end+L) = y; 
   comet = flipud(get(gca,'children')); 
   set(comet(1),'xdata',Y(1,end),'ydata',Y(2,end),'zdata',Y(3,end)); 
   set(comet(2),'xdata',Y(1,2:end),'ydata',Y(2,2:end),'zdata',Y(3,2:end)) 
   set(comet(3),'xdata',Y(1,1:2),'ydata',Y(2,1:2),'zdata',Y(3,1:2)) 
   Y(:,1:L) = []; 
   drawnow; 
 
   % Pause and restart 
 
   uics = flipud(get(gcf,'children')); 
   paws = uics(1); 
   stop = uics(2); 
   rhopick = uics(4); 
   rho = A(3,2); 
   while get(paws,'value')==1 & get(stop,'value')==0 
      set(paws,'string','resume'); 
      drawnow; 
   end 
   set(paws,'string','pause') 
   fin = get(stop,'value') | get(rhopick,'value')==rho; 
   if fin 
      set(paws,'value',0,'string','restart','callback','lorenzgui') 
      set(stop,'value',0,'string','close', ... 
         'callback','cameratoolbar(''close''), close(gcf)') 
   end 
end