function edlorenz(action,kind), % This demo animates the integration of the three coupled nonlinear % differential equations that define the "Lorenz Attractor", a chaotic % system first described by Edward Lorenz of the MIT. % As the integration proceeds you will see: % Picture on the left: a point moving around in a curious orbit in 3-D % space known as a strange attractor. The orbit is bounded, but not % periodic and not convergent. C represents horizontally averaged % temperature mode. % Picture on the right: point is in time-parameter axes associated with the % stream function mode A. It can be switched to display the tent map. % Picture in the center: the temperature difference between actual % temperature and temperature in absence of convection. % % References: % Lorenz, E.N., Deterministic non-periodic flow, % J. Atmos. Sci., 20, 130-141, 1963 % Ott, Edward, Chaos in dynamical systems, % Cambridge University Press, 1993. % File edlorenz.m modified from lorenz.m Matlab-file % by Svetlana Panasyuk, for educational purposes of MIT class % on GeoSystems % last modified: August 2000 % Related questions can be addressed to: % Svetlana Panasyuk % e-mail: panasyuk@alum.mit.edu if nargin<1,action = 'initialize';kind = 1;end; %====================================================================== if strcmp(action,'initialize'),% set up the input graphical window %====================================================================== oldFigNumber = watchon; pos = [50 50 560 420;50 100 700 700]; figNumber = figure('Position',pos(kind,:),'Name','Lorenz Attractor', ... 'NumberTitle','off','Sharecolors','off','Renderer','painters', ... 'Visible','off'); MaaxHndl = axes('Units','normalized','Position',[.05 .10 .75 .95], ... 'Visible','off'); %------------------------------------------------------------------------ % Information for general descriptive and immediate action buttons: %------------------------------------------------------------------------ %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The CONSOLE frames, their titles and the desciptive texts: % if kind == 1, text(0,0,'Set initial conditions and Press the "Start" button ', ... 'HorizontalAlignment','center'); text(0,-0.1,'to begin the integration', ... 'HorizontalAlignment','center'); axis([-1 1 -1 1]); frmHndl = uicontrol('Style','frame','Units','normalized', ... 'Position',[.78 .08 .21 .4],'BackgroundColor',[1 1 1]*.5); t = uicontrol('Style','text','Units','normalized', ... 'Position',[.8 .425 .13 .04],'String','Input','BackgroundColor',[1 1 1]*.5); %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The INFO/Start/Repeat/Close buttons % strings = str2mat('Info','Start','Repeat','Close'); btnpos = [.05 .02 .18 .05;.795 .02 .18 .05;.55 .02 .18 .05;.3 .02 .18 .05]; clbacks = str2mat(['edlorenz(''info'',1)'],['edlorenz(''start'',1)'],... ['edlorenz(''repeat'',1)'],'close(gcf)'); for n = 1:length(strings(:,1)), var = uicontrol('Style','push','Units','normalized', ... 'Position',btnpos(n,:),'String',strings(n,:),'CallBack',clbacks(n,:)); eval([strings(n,1:4) 'Hndl = var;']); end; %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The InputInfo (II) bottons and the Input bottons % strings = str2mat('A','B','C','t','r'); A = 5;B = -5;C = 20;t = 10;r = 28; for n = 1:5, InInHndl(n) = uicontrol('Style','push','Units','normalized', ... 'Position',[.79 .45-.07*n .09 .05],'String',strings(n,:), ... 'CallBack',['edlorenz(''InputInfo' num2str(n) ''')']); InpuHndl(n) = uicontrol('Style','Edit','Units','normalized', ... 'Position',[.88 .45-.07*n .1 .05],'String',num2str(eval(strings(n,:)))); end; else InfoHndl=-1;InInHndl=-ones(1,5);InpuHndl=-ones(1,5); StarHndl=-1;RepeHndl=-1;ClosHndl=-1;axis([-1 1 -1 1]); end; %------------------------------------------------------------------------ % The axes for the plots: BaaxHndl displays parameter space trajectory % TmaxHndl displays the temperature variations vs time, % AdaxHndl displays the advection temperature variations %------------------------------------------------------------------------ ABC = uicontrol('Style','Text','Units','normalized',... 'Position',[.4 .93 .4 .05],'HorizontalAlignment','left'); strings = str2mat('Baax','Tmax','Adax','Teax'); btnpos = [.42 .6 .3 .3;.07 .6 .3 .3;.02 .07 0.8 0.4;.78 .6 .21 .3]; labels = str2mat('C','A','time','temperature variations (B)',' ',' ',... 'max(i)','max(i+1)'); for n = 1:length(strings(:,1)), var = axes('Position',btnpos(n,:),'XTick',[],'YTick',[], ... 'DrawMode','normal','NextPlot','add','Visible','off'); xlabel(labels(2*n-1,:));ylabel(labels(2*n,:)); eval([strings(n,1:4) 'Hndl = var;']); end; set(BaaxHndl,'XTick',0,'YTick',0,'ZTick',0,'View',[-37.5,30]); set(AdaxHndl,'Box','off'); %---------------------------------------------------------------------- % Uncover the figure: 1-14 BtnHndl, 15-19 AxesHndl , 20 ABC output %---------------------------------------------------------------------- hndlList = [InInHndl InpuHndl InfoHndl StarHndl RepeHndl ClosHndl ... MaaxHndl BaaxHndl TmaxHndl AdaxHndl TeaxHndl ABC]; set(figNumber,'Visible','on','UserData',hndlList); watchoff(oldFigNumber); figure(figNumber); if kind == 2,edlorenz('start',2);end; %=========================================================================== elseif strcmp(action,'start'), % To start the demo %=========================================================================== global SIGMA BETA RHO;if kind == 2,global tfinal y0 ANIM;end; %---------------------------------------------------------------------- % gather all the input information and prepare for the run %---------------------------------------------------------------------- strings = str2mat('InInHnd','InpuHndl','InfoHndl','StarHndl','RepeHndl', ... 'ClosHndl','MaaxHndl','BaaxHndl','TmaxHndl','AdaxHndl','TeaxHndl','ABC'); hndlList = get(gcf,'UserData'); for n = 3:length(strings(:,1)), eval([strings(n,:) ' = hndlList(8+n);']); end; InInHnd = hndlList(1:5);InpuHndl = hndlList(6:10); set(ABC,'String',' '); axes(MaaxHndl);cla; text(-.5,-.25,'Advection Temperature Field'); set([BaaxHndl TmaxHndl TeaxHndl],'Visible','off','Visible','on'); %--------------------------------------------------------------------- % Start of Demo % The orbit ranges chaotically back and forth around two different points, % or attractors. It is bounded, but not periodic and not convergent. %---------------------------------------------------------------------- %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The initial conditions & run parameters % %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ SIGMA = 10;BETA = 8/3; if kind == 1,% nice results are for y0 = [20 5 -5] y0(1,1) = str2num(get(InpuHndl(3),'String')); y0(2,1) = str2num(get(InpuHndl(1),'String')); y0(3,1) = str2num(get(InpuHndl(2),'String')); RHO = str2num(get(InpuHndl(5),'String')); tfinal = str2num(get(InpuHndl(4),'String')); elseif kind == 2,% given: RHO tfinal y0 ANIM end; pow = 1/3; t0 = 0;t = t0; hmax = (tfinal - t)/5;hmin = (tfinal - t)/200000;h = (tfinal - t)/100; y = y0; tol = 0.001;tau = tol*max(norm(y,'inf'),1); CABT = [y0;t0]; %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The main loop, do while until hmin % while h >= hmin, if t + h > tfinal, h = tfinal - t; end % Compute the slopes s1 = [-BETA 0 y(2);0 -SIGMA SIGMA;-y(2) RHO -1]*y;ys = y+h*s1; s2 = [-BETA 0 ys(2);0 -SIGMA SIGMA;-ys(2) RHO -1]*ys;ys= y+h*(s1+s2)/4; s3 = [-BETA 0 ys(2);0 -SIGMA SIGMA;-ys(2) RHO -1]*ys; % Estimate the error and the acceptable error delta = norm(h*(s1 - 2*s3 + s2)/3,'inf'); tau = tol*max(norm(y,'inf'),1.0); ts = t;ys = y; if delta <= tau,% Update the solution only if the error is acceptable t = t + h;y = y + h*(s1 + 4*s3 + s2)/6; CABT = [CABT [y;t]]; end if delta ~= 0.0,% Update the step size h = min(hmax, 0.9*h*(tau/delta)^pow); end end; %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % End of Demo % set(MaaxHndl,'UserData',CABT); ys = round(CABT(1:3,length(CABT))*10)/10; set(ABC,'String',['Final: A=' num2str(ys(2)) ' B=' num2str(ys(3)) ... ' C=' num2str(ys(1)) ]); edlorenz('repeat',kind); %=========================================================================== elseif strcmp(action,'repeat'); %=========================================================================== global SIGMA BETA RHO;if kind == 2,global tfinal y0 ANIM;end; strings = str2mat('InInHnd','InpuHndl','InfoHndl','StarHndl','RepeHndl', ... 'ClosHndl','MaaxHndl','BaaxHndl','TmaxHndl','AdaxHndl','TeaxHndl','ABC'); hndlList = get(gcf,'UserData'); for n = 3:length(strings(:,1)), eval([strings(n,:) ' = hndlList(8+n);']); end; InInHnd = hndlList(1:5);InpuHndl = hndlList(6:10); CABT = get(MaaxHndl,'UserData');len = length(CABT); L = (len<11)*round(len/2) + (len>=11)*10;% Save L steps and plot like a comet tail. mima = [min(CABT.');max(CABT.')];nums = find(diff(mima)==0); for ii = 1:length(nums),mima(:,nums(ii)) = [-1;1]*1e-14 + mima(1,nums(ii));end; %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % Define the comet on both upper axes and calculate&plot the tent map % axes(BaaxHndl);cla; set(BaaxHndl,'Drawmode','fast','Visible','off','Visible','on', ... 'XLim',mima([1 2],1),'YLim',mima([1 2],2),'ZLim',mima([1 2],3), ... 'XGrid','on','YGrid','on','ZGrid','on','nextplot','add'); xlabel('C');ylabel('A');zlabel('B'); head = line('color','r','Marker','.','MarkerSize',25,'EraseMode','xor',... 'xdata',CABT(1,1),'ydata',CABT(2,1),'zdata',CABT(3,1)); body = line('color','r','linestyle','-','EraseMode','none',... 'LineWidth',1,'xdata',[],'ydata',[],'zdata',[]); tail = line('color','b','linestyle','-','EraseMode','none',... 'LineWidth',1,'xdata',[],'ydata',[],'zdata',[]); axes(TmaxHndl);cla; set(TmaxHndl,'Drawmode','fast','YLim',mima([1 2],3),'XLim',[0 mima(2,4)],... 'YTick',[0],'YGrid','on','nextplot','add'); xlabel('time');ylabel('temperature variations (B)'); thead = line('color','r','Marker','.','MarkerSize',20,'EraseMode',... 'xor','xdata',CABT(4,1),'ydata',CABT(3,1)); tbody = line('color','r','linestyle','-','EraseMode','none',... 'xdata',[],'ydata',[]); ttail = line('color','b','linestyle','-','EraseMode','none',... 'xdata',[],'ydata',[]); axes(TeaxHndl);cla; set(TeaxHndl,'Drawmode','fast','YTick',[]); xlabel('max(i)');ylabel('max(i+1)');title('1D tent map'); y = abs(CABT(3,len:-1:1).');y_ = y.*(y>mean(y)); yy = y_(1:length(y)-1)-y_(2:length(y)); num1 = (yy>0);num2 = (yy<0); ln = length(num1)-1;NN = find(num2(1:ln).*num1(2:ln+1));ln = length(NN); if ln>2, plot(y(NN(2:ln)),y(NN(1:ln-1)),'ro'); mimay=[min(y(NN(2:ln))) min(y(NN(1:ln-1)));max(y(NN(2:ln))) max(y(NN(1:ln-1)))]; set(TeaxHndl,'AspectRatio',[NaN 1],'XLim',[mimay(1,1) mimay(2,1)], ... 'YLim',[mimay(1,2) mimay(2,2)]); if RHO > 40, text(mimay(1,1)+1,round(mimay(2,1)/1.2),'tent-map algorithm might'); text(mimay(1,1)+1,round(mimay(2,1)/3),'be unstable for this Ra#'); end; axis([min(y(NN)) max(y(NN)) min(y(NN)) max(y(NN))]); else axis([0 CABT(4,len) min(CABT(3,:)) max(CABT(3,:))]); text(1,0,'not enough data'); end; %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The parameters to display advection field and define colormap % l = 2/sqrt(4/BETA-1);coef1 = (4+l*l)/sqrt(2)/l;coef2 = (pi*coef1).^3/sqrt(2)/l; xx = [0:7]/7/8;yy = ([0:15]/15).'; lenx = length(xx);leny = length(yy);lenxy = lenx*leny; s = coef2*sin(2*pi*yy)*ones(size(xx)); cs = sqrt(2)*coef2*sin(pi*yy)*cos(2*pi*l*xx); axes(AdaxHndl);cla; junk = zeros(leny,lenx);Theta = junk;Ksi = junk;junk(:) = [1:lenxy].'; pc=pcolor([junk junk(:,8:-1:1) junk junk(:,8:-1:1)]); set(pc,'EraseMode','none','EdgeColor','none','cdata',[]); Theta = CABT(1,1)*s - CABT(3,1)*cs; mimaT = [min(min(Theta)) max(max(Theta))]; if diff(mimaT) == 0,Theta = ones(size(Theta)); else Theta = Theta-mimaT(1);Theta = floor(Theta/max(max(Theta))*(lenxy-1))+1; end; colormap('default');colormap(hot); MyColorMap = (interp1([1:64]/64,colormap,[3:130]/130)); set(gcf,'MinColormap',length(MyColorMap),'colormap',MyColorMap); %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % The drawing loop % junk = [CABT(:,1)*ones(1,L) CABT]; if kind == 1,ANIM = 1;elseif kind == 2, ANIM = (len-1)*(ANIM~=1)+1;end for ii = L + [ANIM:len], if abs(junk(1,ii))>0.23 & abs(junk(3,ii))>0.23,% to ensure stable colormap Theta = junk(1,ii)*s - junk(3,ii)*cs; Theta = Theta-min(min(Theta)); Theta = floor(Theta/max(max(Theta))*(lenxy-1))+1; set(pc,'cdata',[Theta Theta(:,8:-1:1) Theta Theta(:,8:-1:1)]); end; vec1 = ii-[1 0];vec2 = ii-[L L-1]; set(thead,'xdata',junk(4,ii),'ydata',junk(3,ii)); set(tbody,'xdata',junk(4,vec1),'ydata',junk(3,vec1)) set(ttail,'xdata',junk(4,vec2),'ydata',junk(3,vec2)) set(head,'xdata',junk(1,ii),'ydata',junk(2,ii),'zdata',junk(3,ii)) set(body,'xdata',junk(1,vec1),'ydata',junk(2,vec1),'zdata',junk(3,vec1)) set(tail,'xdata',junk(1,vec2),'ydata',junk(2,vec2),'zdata',junk(3,vec2)) drawnow; end; vec1 = len + [-L:-1];vec2 = [1:len-L];% plot all over again, for print option set(tail,'xdata',CABT(1,vec2),'ydata',CABT(2,vec2),'zdata',CABT(3,vec2)); set(body,'xdata',CABT(1,vec1),'ydata',CABT(2,vec1),'zdata',CABT(3,vec1)) set(head,'xdata',CABT(1,len),'ydata',CABT(2,len),'zdata',CABT(3,len)); set(ttail,'xdata',CABT(4,vec2),'ydata',CABT(3,vec2)); set(tbody,'xdata',CABT(4,vec1),'ydata',CABT(3,vec1)); set(thead,'xdata',CABT(4,len),'ydata',CABT(3,len)); axes(MaaxHndl); texts = str2mat('t_0=','t=','A=','B=','C=','A=','B=','C=','Ra='); strings = num2str([CABT(4,[1 len]).';CABT([2 3 1],1);CABT([2 3 1],len);RHO],3); pos = [.95 -.4;.95 -.7;1.2 -.4;1.2 -.5;1.2 -.6;1.2 -.7;1.2 -.8;1.2 -.9;1.2 -1]; for n=1:length(pos(:,1)),text(pos(n,1),pos(n,2),[texts(n,:) strings(n,:)]);end; %=========================================================================== elseif strcmp(action,'info'); %=========================================================================== ttlStr = 'Lorenz Attractor Info'; hlpStr = ... [' This demo animates the integration of the three ' ' coupled nonlinear differential equations that ' ' define the "Lorenz Attractor", a chaotic system ' ' first described by Edward Lorenz of the MIT. ' ' As the integration proceeds you will see: ' ' Picture on the left: a point moving around in a ' ' curious orbit in 3-D space known as a strange ' ' attractor.The orbit is bounded, but not periodic ' ' and not convergent. C represents horizontally ' ' averaged temperature mode. ' ' Picture on the right: point is in time-parameter ' ' associated with the stream function mode A. ' ' Picture in the center: the temperature difference' ' between actual temperature and temperature in ' ' absence of convection. ' ' ' 'File edlorenz.m modified from lorenz.m Matlab-file' ' by SvetLana Panasyuk while at MIT, for' 'educational purposes of MIT class on GeoSystems ' ' panasyuk@alum.mit.edu ']; helpfun(ttlStr,hlpStr); %=========================================================================== elseif strcmp(action,'InputInfo1'); %=========================================================================== RttlStr = 'A-parameter Info'; RhlpStr = ... [' A = 5 is a good example ' ' ']; helpfun(RttlStr,RhlpStr); %=========================================================================== elseif strcmp(action,'InputInfo2'); %=========================================================================== RttlStr = 'B-parameter Info'; RhlpStr = ... [' B = -5 is a good example ' ' ']; helpfun(RttlStr,RhlpStr); %=========================================================================== elseif strcmp(action,'InputInfo3'); %=========================================================================== RttlStr = 'C-parameter Info'; RhlpStr = ... [' C = 20 is a good example ' ' ']; helpfun(RttlStr,RhlpStr); %=========================================================================== elseif strcmp(action,'InputInfo4'); %=========================================================================== RttlStr = 'Maximal time Info'; RhlpStr = ... [' t = 10 is a good example ' ' ']; helpfun(RttlStr,RhlpStr); %=========================================================================== elseif strcmp(action,'InputInfo5'); %=========================================================================== RttlStr = 'Rayleigh Number Info'; RhlpStr = ... [' r = Ra/Ra_c is Rayleigh number relative to its ' ' critical value (which equals to 657.5 for the ' ' onset of thermal convection in a fluid layer ' ' heated from below. ' ' ' ' important r values are: ' ' 0chaotic attractor' ' 24.06 24.74 Lorenz Chaotic Attractor stable ']; helpfun(RttlStr,RhlpStr); end;