Lorenz 系统的相图的模拟程序Matlab实现

%这是Lorenz 系统的相图的模拟程序振动资讯HO
Fl^h

%   Adapted for Demo by Ned Gulley, 6-21-93; jae Roh, 10-15-96振动资讯:f Fz:L$t Id

%   Copyright 1984-2000 The MathWorks, Inc.

} @q/d6to S0

%   $Revision: 5.11 $  $Date: 2000/06/01 03:46:51 $振动资讯*f,|)N} B8KlV*E

% Author: Thomas Lee

4Ui,@_\`i0

% E-mail: lixf1979@126.com振动资讯 ~&ZM.Za
e fge1]3A4yW

% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China振动资讯\bg2s/O1WZ/[

 

p?*Ck(A,i.iQWw0

function ydot = lorenzeq(t,y)

&|
i;FR3C S)]U0

%LORENZEQ Equation of the Lorenz chaotic attractor.振动资讯n`a.yvF

%   ydot = lorenzeq(t,y).

#~IO@M6a&@n3Y0

%   The differential equation is written in almost linear form.

G'IAFSo!V{#@C9z0

global SIGMA RHO BETA振动资讯|_'v:?,Bt/J

A = [ -BETA    0     y(2) 

C^U!y6J7a8sWb0

         0  -SIGMA   SIGMA 振动资讯n ^8Kz:I4g

      -y(2)   RHO    -1  ];振动资讯Xl@H.w

ydot = A*y;振动资讯$j5|*q3w;C-K

 振动资讯6Jfju3qI

 振动资讯7YO:uzl1x#Ow

function lorenz(action)振动资讯{a3B\yY(i} O'E;^

%LORENZ Plot the orbit around the Lorenz chaotic attractor.

']'{_'E7}1fCIQj[$J0

%   This demo animates the integration of  the   

eL:il!Q%E4I;T0

%   three coupled nonlinear differential equations振动资讯!B3ByRNG:o
k

%   that define the "Lorenz Attractor", a chaotic 振动资讯#Nf:DC*|c6B&s

%   system first described by Edward Lorenz of    振动资讯_1K*FH?
L6N?

%   the Massachusetts Institute of Technology.    振动资讯*{n;Gb u_!tG R}

%   As the integration proceeds you will see a    振动资讯;r,zB'w:m2I5VN

%   point moving around in a curious orbit in    

~kU3nT5?/c0

%   3-D space known as a strange attractor. The   振动资讯llP-?VO{/f*nb

%   orbit is bounded, but not periodic and not   

P:Vcj*Fg#B0

%   convergent (hence the word "strange").       

E-ior2X_[0

%   Use the "Start" and "Stop" buttons to control 振动资讯}K+M)b1VT

%   the animation.                                振动资讯!A+dm2v1c6R

振动资讯h;kKM9I5OJOw

% Possible actions: 振动资讯%wN`^%z_a/r3X(Mn

;yon#?1\0

% initialize

vl"ZTo0

% close振动资讯-mBTn(x(GA

% Information regarding the play status will be held in振动资讯/x JDT]R%P W+Z

% the axis user data according to the following table:

;X1q*mR%V8G[~ w!vC0

play= 1;振动资讯
EM RIN

stop=-1;

5G$|Y;_I iw$a0

if nargin<1,振动资讯D([k+R@"A#@&z

    action='initialize';

7jX&[t'R K,Q0

end;

:V;bo4h5[d)rFn0

if strcmp(action,'initialize'),

:c0k^u&oi-t/ikXYAL0

    oldFigNumber=watchon;振动资讯:A2l gc+~5DMq

    figNumber=figure( ...振动资讯
~*~4Oz!?({Oi

        'Name','Lorenz Attractor', ...振动资讯ld^H
Uj;r`

        'NumberTitle','off', ...振动资讯2O*yi'h)mJ1RN

        'Visible','off');振动资讯 L;a,vxz

    colordef(figNumber,'black')

E/Vh+Y\0

    axes( ...

yWf-u.r}^I0

        'Units','normalized', ...振动资讯%E1x ~qyf a;^K

        'Position',[0.05 0.10 0.75 0.95], ...振动资讯 [m0FT%I_!m.r2x

        'Visible','off');振动资讯
C'H \'O;Q;C0O
i

   text(0,0,'Press the "Start" button to see the Lorenz demo', ...振动资讯$C)cA7l$u

        'HorizontalAlignment','center');

0R${D6\7\g0

    axis([-1 1 -1 1]);振动资讯Az^:\E0d+Nx

   %===================================

(~!|}Q*j,k$@!T;b0

    % Information for all buttons

(A/p(_8C0i0

    labelColor=[0.8 0.8 0.8];

9H/V%dq/s}0

    yInitPos=0.90;振动资讯o{IS0s} I

    xPos=0.85;

X:tc0e? F#QG;BH_UN0

    btnLen=0.10;振动资讯"fD"g)\\Y,e&d*Zb

    btnWid=0.10;

)h0TOr3H1[,_0

    % Spacing between the button and the next command's label

w-N8V$G)V;Wj9b(n0

    spacing=0.05;振动资讯m}J8U @,BW

   %====================================

$g;O4y)c1@.y Y2RO0

    % The CONSOLE frame

FC+E F1G A-[0

    frmBorder=0.02;

\T#q"^4gp0

    yPos=0.05-frmBorder;振动资讯jvx:J;Q/V$u0}

    frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];振动资讯?'z6wrh.m?5^

    h=uicontrol( ...

\ {'j9QtffN0

        'Style','frame', ...

`F!a,PE"]d/D.W1t h6Hk(L0

        'Units','normalized', ...

H&w)Pj+qK9v Cf0

        'Position',frmPos, ...

'E$}uqkR6DWz0

        'BackgroundColor',[0.50 0.50 0.50]);振动资讯_{0No;c9[$d

      %====================================

^%ep!U)f:F$AJC:C*V0

    % The START button振动资讯GS4pC@l

    btnNumber=1;

d3C ?L-DPj0

    yPos=0.90-(btnNumber-1)*(btnWid+spacing);

+T/C2V6w-V.W0

    labelStr='Start';

U2f)aE2[P0

    cmdStr='start';

+dfJ]I7r:W0

    callbackStr='lorenz(''start'');';振动资讯y)vfCE)f

    % Generic button information振动资讯8vq7@5R^'I O

    btnPos=[xPos yPos-spacing btnLen btnWid];

LKZW)l8f0

    startHndl=uicontrol( ...振动资讯3p#F7^k
Vnw

        'Style','pushbutton', ...

q!~9L;ve#|*\^7Q*T6q-`0

        'Units','normalized', ...振动资讯WS*e'i"[E

        'Position',btnPos, ...振动资讯Uo.W x3YO/\%x.x

        'String',labelStr, ...

TU F,V3XcQE0

        'Interruptible','on', ...振动资讯g[w9[7|c

        'Callback',callbackStr);

^0S-i3m/u:hN0

    %====================================

6_bl r*YB8I[||0

    % The STOP button

5O8q_U8_+w0

    btnNumber=2;振动资讯6a(g m}M2|

    yPos=0.90-(btnNumber-1)*(btnWid+spacing);

C.Kr c`)jNvy0

    labelStr='Stop';振动资讯{Hn2] w)s&\

    % Setting userdata to -1 (=stop) will stop the demo.

a0d:g{YG9xt0

    callbackStr='set(gca,''Userdata'',-1)';振动资讯;D4b"m$g1au(B~"c

    % Generic  button information振动资讯7dHir3S(A#Y

    btnPos=[xPos yPos-spacing btnLen btnWid];

n$j,x&~_xi0

    stopHndl=uicontrol( ...

2Mh0o9R-yj0

        'Style','pushbutton', ...

\:}(L,mr/DSss*t&d)d I.j0

        'Units','normalized', ...振动资讯 h'R}}8uX7zW"w

        'Position',btnPos, ...振动资讯8T4n$Wo$AOe%F;d

        'Enable','off', ...

2Q%e1o7i"? T.`0

        'String',labelStr, ...振动资讯7DzNQ(Fw@#q/tt^8a

        'Callback',callbackStr);振动资讯:]$I9TiP%C1G

    %====================================

z;nj/w4I
m0

    % The INFO button振动资讯A8lM8L"VB AK

    labelStr='Info';振动资讯0SN+zi/Gb!D]

    callbackStr='lorenz(''info'')';振动资讯]"h"p*hJ

    infoHndl=uicontrol( ...

d8~{m5wS\!P-Y_-}%p0

        'Style','push', ...

#]mi~+qW0

        'Units','normalized', ...振动资讯 yG?oWj Iy

        'position',[xPos 0.20 btnLen 0.10], ...振动资讯#|1z.TK*A)qn

        'string',labelStr, ...

l*y!]9v8d|_0

        'call',callbackStr);振动资讯%K?"]!B#wt

    %====================================振动资讯c^}&^@0O)L

    % The CLOSE button振动资讯qvk^9W8v C

    labelStr='Close';振动资讯Uy:Bf [X'E2JZ

    callbackStr='close(gcf)';

(Dq+T~} lrc0

    closeHndl=uicontrol( ...

vY`L1^ B%~t0

        'Style','push', ...振动资讯WiON*n+j

        'Units','normalized', ...振动资讯~'Z%i!B+l9XQ

        'position',[xPos 0.05 btnLen 0.10], ...

(q+qHR7olx;jd`^0

        'string',labelStr, ...

(_lKp!^#y0

        'call',callbackStr);振动资讯3`,q^~ n
a"d,Lh.F

      % Uncover the figure振动资讯}3sWK|]A q$^r

    hndlList=[startHndl stopHndl infoHndl closeHndl];

9^0o6n%nJ.Q0

    set(figNumber,'Visible','on', ...振动资讯 s_Oz$l#_

        'UserData',hndlList);

Xp[b9P5F&b,~?$x0

    watchoff(oldFigNumber);振动资讯)\(J5};}4k

    figure(figNumber);振动资讯)tp0g4K,T!f
e
Y `(xc

elseif strcmp(action,'start'),

VV D!`Q+Q0

    axHndl=gca;振动资讯-`Y$AaAjm&J [ Q K

    figNumber=gcf;

M
I.T*I
as/J1K0

    hndlList=get(figNumber,'UserData');振动资讯Q!fa[Z? B|"rg

    startHndl=hndlList(1);

Q3o l*g|wE)h#dH0

    stopHndl=hndlList(2);

OW!L"G}t+q)G0

    infoHndl=hndlList(3);振动资讯:](sB#p1I p"I2u~.X

    closeHndl=hndlList(4);振动资讯+vu/s@U
p5b

    set([startHndl closeHndl infoHndl],'Enable','off');

bGm$E!I*A*?0

    set(stopHndl,'Enable','on');振动资讯5E'd(}.We8M

    % ====== Start of Demo

E9Ksd8Ke0

    set(figNumber,'Backingstore','off');振动资讯$Hw+qE)|*gN+{

    % The graphics axis limits are set to values known 振动资讯I?N V3V+R5Y

    % to contain the solution.

\ b"f:o t0D q
g0

    set(axHndl, ...

i!OB q5A6R'[}0

        'XLim',[0 40],'YLim',[-35 10],'ZLim',[-10 40], ...

)~wB5rg.HS0

        'Userdata',play, ...振动资讯{/L!J8X1^kbJ

        'XTick',[],'YTick',[],'ZTick',[], ...

${)QO6a'@:_0

        'Drawmode','fast', ...

G-?0G"]6Dz0

        'Visible','on', ...

xs3_d0b)WyAi0

        'NextPlot','add', ...

9N7|9]'? h0XD'h:]0

        'Userdata',play, ...

g&X"s
\#D1m$o&nn0

        'View',[-37.5,30]);

Jr2@8]:q2E!k;a:@y0

    xlabel('X');振动资讯-g~b/|8w'jUt{5W

    ylabel('Y');振动资讯O!eA8H.UhSrT?

    zlabel('Z');

Z8NYd^ m}C0

    % The values of the global parameters are振动资讯{
@jei!~`Q

    global SIGMA RHO BETA振动资讯jAIdUCP(fK k(c

    SIGMA = 10.;

,M[Z5KT!c/{5J3|bD0

    RHO = 28.;振动资讯\J0_X.n1B

    BETA = 8./3.;振动资讯aafNO!p4} xJZ

    % The orbit ranges chaotically back and forth around two different points,

Sv+cssy0

    % or attractors.  It is bounded, but not periodic and not convergent.振动资讯*l`d8L-F h;V1E u

    % The numerical integration, and the display of the evolving solution,振动资讯XN!|J,m"{-e"^ C;A,{

    % are handled by the function ODE23P.振动资讯|8}Ic(]?

     FunFcn='lorenzeq';

lOiE7?~^0P0

    % The initial conditions below will produce good results

cg t"CH;_u}m.SZ0

    %y0 = [20 5 -5];

K,]{.wKp|0

    % Random initial conditions

2X%_ |:\5M0

    y0(1)=rand*30+5;

#L[%~:kH0

    y0(2)=rand*35-30;

A/X7G1P/gj3K-M
X0

    y0(3)=rand*40-5;振动资讯.AeR;juL

    t0=0;振动资讯-rWF&R7f'i;hg
@

    tfinal=100;振动资讯e4\T5M1w@

    pow = 1/3;

%S@Ca5B5c(IX0

    tol = 0.001;振动资讯4HY1q!eWD3|

    t = t0;振动资讯Z)DUJC)Wx!XWLu

    hmax = (tfinal - t)/5;振动资讯3i;b$WT-K }^

    hmin = (tfinal - t)/200000;

&m/b'e&L8d"Jm Bn(Q0

    h = (tfinal - t)/100;振动资讯y:H{;P9\

    y = y0(:);振动资讯+eh7e#Y/~ U5J

    tau = tol * max(norm(y,'inf'),1);振动资讯B.U3m
E[y

    % Save L steps and plot like a comet tail.振动资讯*]$_l H,l c8oo3m&V

    L = 50;

"rk!t*nz8T3k:fV0

    Y = y*ones(1,L);

WT!u-qu+Ya0

     cla;振动资讯x1ueS0~:S

    head = line( ...振动资讯Q!? l)M'Fwm7n

        'color','r', ...

2`5C mLw6O9v0

        'Marker','.', ...振动资讯 j,sY-GqT

        'markersize',25, ...

k#M'V-B)Z0

        'erase','xor', ...

N5M|]~+|dh0

        'xdata',y(1),'ydata',y(2),'zdata',y(3));振动资讯/NFA9fx5D$K&d7J$_

    body = line( ...

.Bp+d2ZiGBz0

        'color','y', ...振动资讯7Z \l \3h

        'LineStyle','-', ...振动资讯/YN8xb"i

        'erase','none', ...

C6v:ZL}e0

        'xdata',[],'ydata',[],'zdata',[]);

d~$I:z8ry9Cun0

    tail=line( ...

t dct#S*[1C0

        'color','b', ...

P V#n#c.xkn*n
Zyq_0

        'LineStyle','-', ...振动资讯 [J&{'{9WXl

        'erase','none', ...

JS*k.LF*RJxR/Bx6k0

        'xdata',[],'ydata',[],'zdata',[]);振动资讯 VGA AGV)|s|

     % The main loop振动资讯^cJ6I#g)@P

    while (get(axHndl,'Userdata')==play) & (h >= hmin)

GfG;S}pZ;]0

        if t + h > tfinal, h = tfinal - t; end

$eJRD&Eol0

        % Compute the slopes

/q#? \}d n#Nv0

        s1 = feval(FunFcn, t, y);

U]
\]{Q!Jf0

        s2 = feval(FunFcn, t+h, y+h*s1);振动资讯{fFXya#sJ

        s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4);

$]KS"Mq{UnB"E0

         % Estimate the error and the acceptable error振动资讯8\M:~3_O5td

        delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');振动资讯 G!U OM1XD

        tau = tol*max(norm(y,'inf'),1.0);振动资讯&bdu/}}

        % Update the solution only if the error is acceptable

^"S
\7K!yCy VZx0

        ts = t;

*D&J)`N)N*~5fN |0

        ys = y;

jxHU6gF"^0

        if delta <= tau振动资讯*A_C2B*a j

            t = t + h;

ylIT.r6@c6jn[0

            y = y + h*(s1 + 4*s3 + s2)/6;

#I%^A3wi+T0

            % Update the plot

8FzFL+J-g0

            Y = [y Y(:,1:L-1)];

,}y$S+wC0

            set(head,'xdata',Y(1,1),'ydata',Y(2,1),'zdata',Y(3,1))振动资讯+oan;acB7V

            set(body,'xdata',Y(1,1:2),'ydata',Y(2,1:2),'zdata',Y(3,1:2))振动资讯H"y)Qo"m?`;iC\

            set(tail,'xdata',Y(1,L-1:L),'ydata',Y(2,L-1:L),'zdata',Y(3,L-1:L))

p h-n K4GSPuy*_2{c0

            drawnow;

:se$i;d^ swdE"Ap0