130 JOURNAL OF THE ATMOSPHERIC SCIENCES VOLUME 20 Deterministic Nonperiodic Flow1 EDWARD N.LORENZ Massachusetts Islilute of Technology (Manuscript received 18 November 1962， in revised form 7 January 1963) ABSTRACT Finite systems of deterministic ordinary nonlinear differential equations may be designed to represent forced dissipative hydrodynamic flow. Solutions of these equations can be identified with trajectories in phase space. For those systems with bounded solutions, it is found that nonperiodic solutions are ordinarily unstable with respect to small modifications, so that slightly differing initial states can evolve into consider- ably different states. Systems with bounded solutions are shown to possess bounded numerical solutions. A simple system representing cellular convection is solved numerically. All of the solutions are found to be unstable, and almost all of them are nonperiodic. The feasibility of very-long-range weather prediction is examined in the light of these results. 1.Introduction Thus there are occasions when more than the statistics Certain hydrodynamical systems exhibit steady-state of irregular flow are of very real concern. In this study we shall work with systems of deter- flow patterns, while others oscillate in a regular periodic ministic equations which are idealizations of hydro- fashion. Still others vary in an irregular, seemingly haphazard manner, and, even when observed for long dynamical systems. We shall be interested principally in periods of time, do not appear to repeat their previous nonperiodic solutions, i.e., solutions which never repeat history. their past history exactly, and where all approximate repetitions are of finite duration. Thus we shall be in- These modes of behavior may all be observed in the volved with the ultimate behavior of the solutions, as familiar rotating-basin experiments, described by Fultz, opposed to the transient behavior associated with et al. (1959)and Hide (1958). In these experiments, a arbitrary initial conditions. cylindrical vessel containing water is rotated about its A closed hydrodynamical system of finite mass may axis, and is heated near its rim and cooled near its center ostensibly be treated mathematically as a finite collec- in a steady symmetrical fashion. Under certain condi- tion of molecules-usually a very large finite collection tions the resulting flow is as symmetric and steady as the -in which case the governing laws are expressible as a heating which gives rise to it. Under different conditions finite set of ordinary differential equations. These equa- a system of regularly spaced waves develops, and pro- tions are generally highly intractable, and the set of gresses at a uniform speed without changing its shape. molecules is usually approximated by a continuous dis- Under stil different conditions an irregular flow pattern tribution of mass. The governing laws are then expressed forms, and moves and changes its shape in an irregular as a set of partial differential equations, containing such nonperiodic manner. quantities as velocity, density, and pressure as de- Lack of periodicity is very common in natural sys- pendent variables. tems, and is one of the distinguishing features of turbu- It is sometimes possible to obtain particular solutions lent flow. Because instantaneous turbulent flow patterns of these equations analytically, especially when the are so irregular, attention is often confined to the sta- solutions are periodic or invariant with time, and, in- tistics of turbulence, which, in contrast to the details of deed, much work has been devoted to obtaining such turbulence, often behave in a regular well-organized solutions by one scheme or another. Ordinarily, how- manner. The short-range weather forecaster, however, ever, nonperiodic solutions cannot readily be deter- is forced willy-nilly to predict the details of the large- mined except by numerical procedures. Such procedures scale turbulent eddies-the cyclones and anticyclones- involve replacing the continuous variables by a new which continually arrange themselves into new patterns. finite set of functions of time, which may perhaps be the values of the continuous variables at a chosen grid of IThe research reported in this work has been sponsored by the points, or the coefficients in the expansions of these Geophysics Research Directorate of the Air Force Cambridge variables in series of orthogonal functions.The govern- Research Center, under Contract No. AF 19(604)-4969. ing laws then become a finite set of ordinary differential
131EDWARD N.LORENZMARCH1963an M-dimensional Euclidean spaceTwhose coordinatesequations again, although a far simpler set than the oneareXi,.,Xm.Eachpointinphasespacerepresentsawhich governs individual molecular motions.In any real hydrodynamical system, viscous dissipa-possible instantaneous state of the system. A statewhich is varying in accordancewith (1)is representedtion is always occurring, unless the system is moving asby a moving particle in phase space, traveling along aa solid, and thermal dissipation is always occurring,trajectory in phase space. For completeness, the positionunless the system is at constant temperature. For cer-of a stationary particle, representing a steady state, istain purposes many systems may be treated as conserva-tive systems, in which the total energy, or some otherincluded as a trajectory.Phase space has been a useful concept in treatingquantity, does not vary with time. In seeking the ulti-mate behavior of a system, the use of conservativefinite systems, and has been used by such mathema-equations is unsatisfactory, since the ultimate value ofticians as Gibbs (1902) in his development of statisticalmechanics,Poincare (1881)in his treatment of the solu-any conservative quantitywould then havetoequal thetions of differential equations, and Birkhoff (1927) inarbitrarily chosen initial value. This difficulty may beobviated by including the dissipative processes, therebyhis treatise on dynamical systems.From the theory of differential equations (e,g., Fordmaking the equations nonconservative, and also includ-1933, ch. 6), it follows, since the partial derivativesing external mechanical or thermal forcing, thus pre-venting the system from ultimately reaching a state ofaF/ax, are continuous, that if tois any time, and ifrest. If the system is to be deterministic, the forcingXio, ...Xmo is any point in F, equations (1) possess afunctions, if not constant with time, must themselvesuniquesolutionvary according to some deterministic rule.(2)X, f(X1o,*-,Xaro,0), i1, **, M,In this work, then, we shall deal specifically withfinite systems of deterministic ordinary differentialvalid throughout some time interval containing to, andequations, designed to represent forced dissipativesatisfying the conditionhydrodynamical systems. We shall study the propertiesof nonperiodic solutions of these equations.(3)f(X10,**,Xmo,fo)-Xio, i=1, -, M.It is not obvious that such solutions can exist at allThe functions f. are continuous in Xio, .., X mo and t.Indeed, in dissipative systems governed by finite sets ofHence there is a unique trajectory through each point oflinear equations, a constant forcing leads ultimately tor. Two or more trajectories may,however, approach thea constant response, while a periodic forcing leads to asamepoint or the same curve asymptoticallyas t-→periodic response. Hence, nonperiodic flow has some-or as t-→--oo. Moreover, since the functions f: aretimes been regarded as the result of nonperiodic orcontinuous, the passage of time defines a continuousrandomforcing.deformation of any region of r into another region.The reasoning leading to these concludions is notIn the familiar case of a conservative svstem, whereapplicable when the governing equations are nonlinear.somepositivedefinitequantityQ,whichmayrepresentIf the equations contain terms representing advectionsome form of energy, is invariant with time, each trathe transport of some property of a fluid by the motionjectory is confined to one or another of the surfaces ofof the fluid itself-a constant forcing can lead to aconstant Q.These surfaces may take the form of closedvariable response.In the rotating-basin experimentsconcentric shells.already mentioned, both periodic and nonperiodic flowIf, on the other hand, there is dissipation and forcing,result from thermal forcing which, within the limits ofand if, whenever O equals or exceeds some fixed valueexperimental control, is constant.Exact periodic solu-Qi, the dissipation acts to diminishQmore rapidly thentions of simplified systems of equations, representingthe forcing can increase Q, then (-dQ/dt) has a positivedissipative flow with constant thermal forcing, havelower bound whereQ≥Qr,and each trajectory mustbeen obtained analytically by the writer (1962a).Theultimately become trapped in the region where Q<Qrwriter (1962b)has also found nonperiodic solutions ofTrajectories representing forced dissipative flow maysimilar systems of equations bynumericalmeans.therefore differ considerably from those representingconservativeflow.2.PhasespaceForced dissipative systems of this sort are typifiedConsider a system whose state may be described by Mby the systemvariables Xr, ., X. Let the system be governed by(4)dX/dt-EauxX,Xx-EbuX,+c,the set of equationsi.k7dX:/dt-F(X,..Xm),i=1,.,M,(1)whereaX,X,Xvanishes identically,buX,X,ispositive definite, and ci, .", c are constants.Ifwhere time t is the single independent variable, and theQ--E"(5)functions F,possess continuous first partial derivatives.Suchasystemmaybestudied bymeansofphasespace
132JOURNAL OF THE ATMOSPHERIC SCIENCESVOLUME20and if er, *.", em are the roots of the equationsabsolute-value signs denote distance in phase space.Because F is continuously deformed as t varies, every(6)E (bu+bj)ej=ci,point on the trajectory through Po is also a limit point2of P(), and the set of limit points of P(t) forms a tra-jectory, or a set of trajectories, called the limiting traitfollows from (4)thatjectories of P(t).A limiting trajectory is obviously contained within R in its entirety.dQ/dt- bueej-Zbu(X;-ei)(X,-e).(7If a trajectory is contained among its own limiting6i.itrajectories, it will be called central; otherwise it will becalled noncentral.A central trajectory passes arbitrarilyThe right side of (7) vanishes only on the surface ofclosely arbitrarily often to any point through which itan ellipsoid E, and is positive only in the interior of Ehas previously passed, and, in this sense at least, sepa-The surfaces of constant Q are concentric spheres. If Srate sufficiently long segments of a central trajectorydenotes a particular one of these spheres whose interiorare statistically similar.Anoncentral trajectory remainsR contains the ellipsoid E, it is evident that each tra-a certain distance away from any point through whichjectory eventually becomes trapped within R.it has previously passed. It must approach its entire setof limit points asymptotically, although it need not3.Theinstability of nonperiodic flowapproach any particular limiting trajectory asymptoti-cally.Its instantaneous distancefrom its closestlimitIn this section we shall establish one of the mostpoint is therefore a transient quantity, which becomesimportant properties of deterministic nonperiodic flow,arbitrarily small as t-→ o.namely,its instability with respect to modifications ofA trajectory P(t) will be called stable at a point P(t)small amplitude. We shall find it convenient to do thisif any other trajectory passing sufficiently close to P(ti)by identifying the solutions of the governing equationsat time ti remains close to P(o) as t-→ co; i.e., P(t) iswith trajectories in phase space. We shall use suchstable at P(t) if for any e>0 there exists a (e,t)>0symbols as P() (variable argument) to denote trajec-such that if [Pi(t)-P(t)/ti, [Pi(t2)tories, and such symbols as P or P(to) (no argument or-P(t2)/0 there exists a s(e)>0 and a timediffers in that Birkhoff was concerned mainly with con-to(e) such that if ti>to and |P(t)-P(t)lti.Pi(ta)-P(t)i0arbitrarily closely arbitrarily often.More precisely,Poisand for any time interval To, there exists a t(e,ro)>roa limit point of P(t)if forany e>0and any timet,thereexists a time t(e,ti)>t such that |P(t2)-Po/tu,[P(t+)-P(t2)
133EDWARDN.LORENZMARCH1963if it is stable at all, its very stability is one of its tran-<e, Periodic trajectories are special cases of quasi-sient properties, which tends to die out as time properiodic trajectories.A trajectory which is not quasi-periodic will be calledgresses.Inviewoftheimpossibilityofmeasuringinitialconditions precisely,and thereby distinguishing betweennonperiodic.If P(0) is nonperiodic,P(t+-) may bea central trajectory and a nearby noncentral trajectory,arbitrarily closetoP(t)forsometimetandsomeall nonperiodic trajectories are effectively unstable fromarbitrarilylargetimeintervalr,but,ifthisisso,P(t+)the point of view of practical prediction.cannot remain arbitrarily close to P(t) as t→ oo, Non-periodic trajectories are of course representations of4.Numerical integration of nonconservative sys-deterministic nonperiodic flow,and form the principaltemssubject of this paper.Periodic trajectories are obviously central. Quasi-The theorems of the last section can be of importanceperiodic central trajectories include multiple periodiconly if nonperiodic solutions of equations of the typetraiectorieswithincommensurableperiods,whilequasi-considered actually exist. Since statistically stationaryperiodic noncentral trajectories include those whichnonperiodic functions of time are not easily describedapproach periodic trajectories asymptotically.Non-analytically, particular nonperiodic solutions can prob-periodic trajectories may be central or noncentral.ably befound most readily by numerical procedures.InWe can now establish the theorem that a trajectorythis section we shall examine a numerical-integrationwith a stable limiting trajectory is quasi-periodic.Forprocedure which is especially applicable to systems ofif Po()isa limiting trajectory of P(t),two distinctequations of the form (4).In a later section we shall usepoints P(t) and P(tr++), with r arbitrarily large, maythis procedure to determine a nonperiodic solution of abe found arbitrary close to any point Po(to). Since Po(0)simple set of equations.is stable, P() and P(t+-) must remain arbitrarilyTo solve (1) numerically we may choose an initialclose to Po(t+t-ti), and hence to each other, as t-→ co,time to and a time increment Zi, and letand P(t) is quasi-periodic.(8)Xi,n=X,(to+nAt).Itfollows immediatelythat a stable central trajectoryis quasi-periodic, or, equivalently, that a nonperiodicWe then introduce the auxiliary approximationscentral trajectory is unstable.(9)Xi(n+1)=Xi,n+F(Pn)A,The result has far-reaching consequences when thesystem being considered is an observable nonperiodic(10)Xi(n+2)= X (n+1)+F(P(n+1)At,system whose future state we may desire to predict. Itwhere Pn and P(n+1) are the points whose coordinates areimplies that two states differing by imperceptibleamountsmayeventuallyevolveinto two considerably(Xi,n,,XMn)and(Xi(n+1),**,XM(n+1)different states. If, then, there is any error whatever inThe simplest numerical procedure for obtainingobserving the present state-and in any real systemapproximate solutions of (1) is the forward-differencesuch errors seem inevitablean acceptable predictionprocedure,of an instantaneous state in the distant future mayXi,n+1=Xi(n+1).(11)well be impossible.As for noncentral trajectories, it follows that a uni-In many instances better approximations to the solu-formly stable noncentral trajectory is quasi-periodic, or,tions of (1) may be obtained by a centered-differenceeguivalently,a nonperiodic noncentral trajectory is notprocedureuniformly stable.The possibility of a nonperiodic non-Xi,n+1=Xi,n-1+2F(P.)At.(12)central trajectory which is stable but not uniformlystable still exists. To the writer, at least, such trajec-This procedure is unsuitable, however, when the detertories, although possible on paper, do not seem charac-ministic nature of (1) is a matter of concern, since theteristic of real hydrodynamical phenomena. Any claimvalues of Xy.n, ..., X M,n do not uniquely determine thethat atmospheric flow,for example, is represented by avalues ofXi,n+1,..",X.n+1.trajectory of this sort would lead to the improbableA procedure which largely overcomes the disadvant-conclusion that we ought to master long-range fore-ages of both the forward-difference and centered-differ-casting as soon as possible, because, the longer we wait,ence procedures is thedouble-approximation procedure,the more difficult our task will become.defined by the relationIn summary,wehave shown that, subject to theX,n+1=Xi,n+[F(P)+F(P(a+)JAl.(13)conditions of uniqueness, continuity, and boundednessprescribed at the beginning of this section, a centralHerethe coefficientof Atis anapproximation tothetimetrajectory, which in a certain sense is free of transientderivative of X,at time le+(n+)t.From (9) andproperties, is unstable if it is nonperiodic. A noncentral(10), it follows that (13) may be rewrittentrajectory, which is characterized by transient prop-erties, is not uniformly stable if it is nonperiodic, and,Xi.n+1=(Xtn+Xi(n+2).(14)
134JOURNAL OF THE ATMOSPHERIC SCIENCESVOLUME20A convenient scheme for automatic computation is thetotically approaching a steady state.A similar resultsuccessive evaluation of X(n+1), Xi(n+2), and X.n+1holds when the double-approximation procedure (14) isaccording to (0), (10) and (14). We have used thisapplied to a conservative system.procedure in all the computations described in thisstudy.5.The convection equations of SaltzmanIn phase space a numerical solution of (1) must beIn this section we 'shall introduce a system of threerepresented by a jumping particle rather than a con-ordinary differential equations whose solutions affordtinuously moving particle. Moreover, if a digital com-the simplest example of deterministic nonperiodic flowputer is instructed to represent each number in itsof which the writer is aware.The system is a simplifica-memory by a preassigned fixed number of bits, onlytion of one derived by Saltzman (1962)to study finite-certain discrete points in phase space will ever be oc-amplitude convection. Although our present interest iscupied. If the numerical solution is bounded, repetitionsin the nonperiodic nature of its solutions, rather thanmust eventually occur, so that, strictly speaking, everyinits contributionstotheconvection problem,weshallnumerical solution is periodic. In practice this considera-describe its physical background briefly.tion may be disregarded, if the number of differentRayleigh (1916) studied the flow occurring in a layerpossible states is far greater than the number of itera-of fuid of uniform depth H,when the temperaturetions ever likely to be performed. The necessity fordifference between the upper and lower surfaces isrepetition could be avoided altogether by the somewhatmaintained at a constant value AT. Such a systemuneconomical procedure of letting the precision ofpossesses a steady-state solution in which there is nocomputation increaseas n increases.motion, and the temperature varies linearly with depth,Consider now numerical solutions of equations (4),If this solution is unstable, convection should develop.obtained by the forward-difference procedure (11).ForIn the case where all motions are parallel to thesuch solutions,x-z-plane, and no variations in the direction of they-axis occur, the governing equations may be writtenQn+1=Qn+(dQ/dt)nAI+ZF2(Pn)AR.(15)(see Saltzman,1962)aa0Let S' be any surface of constant Q whose interior Ra(,V)-(17)++gacontains the ellipsoid E where dQ/dt vanishes, and letata(x,2)aS be any surface of constant Q whose interior R con-tains S'.aa(+,0)ATaySince F? and dQ/dt both possess upper bounds in(18)0-+V20.atHaxa(x,2)R', we may choose At so small that Pn+i lies in R ifPn lies in R'. Likewise, since F possesses an upperHere is a stream function for the two-dimensionalbound and dO/dt possesses a negative upper bound inmotion, is the departure of temperature from thatR-R',wemaychooset so small that Qn+i<Q. if Puoccurring in the state of no convection, and the con-lies in R--R'.Hence Atmaybe chosen so small that anystants g,a,y,and x denote,respectively,theaccelerationjumping particle which has entered R remains trappedof gravity, the coefficient of thermal expansion, thewithin R, and the numerical solution does not blow up.kinematic viscosity,and the thermal conductivity.TheA blow-up may still occur, however, if initially theproblem is most tractable when both the upper andparticle is exterior to R.lower boundaries are taken to be free, in which caseConsider now the double-approximation procedure and vy vanish at both boundaries.(14). The previous arguments imply not only thatRayleigh found that fields of motion of the formP(n+1) lies within Rif P,lies within R,but also thatP(n+2) lies within R if P(n+1) lies within R. Since the(19)=o sin (raH-1a)sin (H-12),region Ris convex, it follows that Pn+1, as given by (14),(20)=,cos (πaH-1a)sin (-12),lies within Rif Pn lies within R.Hence if A is chosen sosmall that the forward-difference procedure does notwould develop if the quantityblow up, the double-approximation procedure also doesRagaHT-1k-1,(21)not blow up.We note in passing that if we apply the forward-now called the Rayleigh number, exceeded a critical valuedifference procedure to a conservative system whereR。=*a-2(1+a2)3(22)dQ/dt-0 everywhere,TheminimumvalueofRe,namely27/4,occurs(16)Qn+1=Qn+F2(Pn)Ae.when a?=2.Saltzman (1962) derived a set of ordinary differentialIn this case, for any fixed choice of At the numericalequations by expanding and in double Fourier seriessolution ultimately goes to infinity, unless it is asymp-in and s, with functions of t alone for coeffcients, and
135MARCH1963EDWARDN.LORENZsubstituting these series into (17) and (18). He arrangedwill convert them to this form. One of the simplest ofthe right-hand sides of the resulting equations in double-these is thetransformationFourier-series form, by replacing products of trigono-(28)X'-X,Y'-Y,Z-Z-r-0.metric functions of (or ) by sums of trigonometricfunctions,and then equated coefficients of similar func-Solutions of (25)-(27) therefore remain bounded withintions of x and z. He then reduced the resulting infinitearegionRasr→co,and thegeneral results of Sectionssystem to a finite system by omitting reference to all2,3and 4applyto theseequations.but a specified finite set of functions of t, in the mannerThe stability of a solution X(+), Y(+), Z(+) may beproposed by the writer (1960).formally investigated by considering the behavior ofHe then obtained time-dependent solutions by nu-small superposed perturbations o(+), yo(+), zo(+). Suchmerical integration. In certain cases all except three ofperturbations are temporarily governed by the linear-the dependent variables eventually tended to zero, andized equationsthese three variables underwent irregular, apparentlynonperiodic fluctuations.0GTxoxo-0These same solutions would have been obtained if the-1 -Xo -| (r-2)yo(29)series had at the start been truncated to include a totalXY-6501Z0of three terms. Accordingly, in this study we shall letSincethecoeficientsin(29)varywithtime,unlessa(1+a2)-k-=XV2 sin (raH-a) sin (μH-1g),(23)the basic state X, Y, Z is a steady-state solution ofR-1RAT-19YV2cos(aH-1)sin(H-1z)(25)-(27), a general solution of (29) is not feasible.However, the variation of the volume V。 of a small-Zsin (2H-1z),(24)region in phase space, as each point in the region iswhere X, Y, and Z are functions of time alone. Whendisplaced in accordance with (25)-(27),is determinedexpressions (23)and (24)are substituted into (17)andby the diagonal sum of the matrix of coefficients;(18), and trigonometric terms other than those occurringspecificallyin (23) and (24) are omitted, we obtain the equationsVo=--(a+b+1)Vo.(30)X--ax+oY,(25)This is perhaps most readily seen by visualizing theY=-Xz+rX-Y(26)motion in phase space as the flow of a fuid, whosedivergence isZ=XY-bz.(27)ax.ayaz.Here a dot denotes a derivative with respect to the+(31)=-(g+6+1).dimensionlesstime=H-2(1+a2)xt,whileg-1yisaxayazthe Prandil number,r=-R-1Ra, and b=4(1+a*)-1,Hence each small volume shrinks to zero as T o, atExcept for multiplicative constants, our variables X,Y,a rate independent of X, Y, and Z.This does not implyand Z are the same as Saltzman's variables A, D, andthat each small volume shrinks to a point; it may simplyG.Equations (25),(26),and (27)are the convectionbecome flattened into a surface, It follows that theequationswhose solutions weshall study.volume of the region initially enclosed by the surfaceSIn these equations X is proportional to the intensityshrinks to zero at this same rate, so that all trajectoriesof the convective motion, while Y is proportional to theultimately become confined to a specific subspacetemperature difference between the ascending and de-having zero volume. This subspace contains all thosescending currents, similar signs of X and Y denotingtrajectories which lie entirely within R, and so containsthat warm fluid is rising and cold fluid is descending.all central trajectories.The variable Z is proportional to the distortion of theEquations (25)-(27) possess the steady-state solutionvertical temperature profile from linearity,a positiveX-Y-Z0, representing the state of no convection.value indicating that the strongest gradients occur nearWith this basic solution, the characteristic equation ofthe boundaries.the matrix in (29) isEquations (25)-(27)maygive realistic results whenthe Rayleigh number is slightly supercritical, but their[+2+(+1)入+(1))=0.(32)solutions cannot be expected to resemble those of (17)and (18) when strong convection occurs, in view of theThis equation has three real roots when r>O; all areextremetruncation.negative when r1.Thecriterion for the onset of convection is therefore r=1,6.Applicationsoflineartheoryor Ra-Re, in agreement with Rayleigh's result.Whenr>1,equations (25)-(27)possess twoadditionalAlthough equations (25)-(27), as they stand, do nothavetheform of (4),a numberof lineartransformationssteady-state solutions XY-±Vb(r-1), Z--1
136JOURNAL OFTHEATMOSPHERICSCIENCESVOLUME 20For either ofthese solutions, the characteristic equation7. Numerical integration of the convection equa-of the matrix in (29) istions入3+(g+6+1)入2+(r+)入+2b(r—1)=0.To obtain numerical solutions of the convection equa-.(33)tions, we must choose numerical values for the con-This equation possesses one real negative root and twostants.Following Saltzman (1962),we shall let =10complex conjugate roots when r>i; the complex con-and a,so thatb=8/3.The critical Rayleigh numberjugate roots are pure imaginary if the product of thefor instability of steady convection then occurs whencoefficients of 入?and 入 equals the constant term,orr=470/19=24.74.We shall choose the slightly supercritical value r-28.r=a(a+b+3)(α-b-1)-1.(34)The states of steady convection are then represented bythe points (6v2, 6V2, 27) and (-6V2, -6V2, 27) in phaseThis is the critical value of rfor the instability of steadyspace, while the state of no convection corresponds toconvection. Thus if gb+1, steady convection is unstable for suffi-We have used the double-approximation procedureciently high Rayleigh numbers. This result of coursefor numerical integration, defined by (9), (10),and (14),applies only to idealized convection governed by (25)-Thevalue Ar=0.01 has been chosen for the dimension(27), and not to the solutions of the partial differentiallesstimeincrement. The computations have beenequations (17)and (18).performed onaRoyal McBeeLGP-30electronic com-The presence of complex roots of (34) shows that ifunstable steady convection is disturbed, the motion willTABILE 2. Numerical solution of the convection equationsoscillate in intensity. What happens when the disturb-Values ofX, Y,Z are given at every iteration Nfor which zances become large is not revealed by linear theory. Topossesses a relativemaximum,for thefirst 6000 iterations.investigate finite-amplitude convection, and to studythe subspace to which trajectories are ultimately con-NXYNxNYZfined, we turn to numerical integration.00450174005504833029011700750352010700910083028730980123007603650168-009202880084317101340082038302300092-008402893268015500690435TABLE 1.Numerical solution of the convection equations.02900292-009200833333-011400790342Values of X, Y, Z are given at every fifth iteration N, for the0354-0093-008302923400011700770350first 160 iterations.04160093008302933468-01250083036104780094-008202953541-01290073037805400296_00940082NxYz3625-0146-007404130602-0095-00820298369501270079037000000000001000000664-009600830300377201360072039400050004001200000726-0097-008303023853-0144-0077040700100009002000000789009700810304392601290072038000150016003600020851009900830307401401480068042100200030006600070914~0100-008103094082-0120-00740359002500540115002409770100-0080031241530129-00780375003000930192007410400080-010203154233040401440082003501500268020111030104-00810319430701350081038500400195023403971167-010500790323441701620069045000450055048312310174-0107-0079032844800106008103240050009700670415129500820333-011145440109008203290055002500930340136101110077033946090110008003340060-0020008902981427-01160079467503470112007603410275149500650046-0084-01200077035747410118008103490070_0083026200720371-00611566-012548100120007403600075-008602561643-0070-0139007703964881013000810376008000770091025517220140007504014963006801410406008500950258179800840135-00720391503501330381-00810090-0089-0098026618820146007404135124015100760422009500930098027519520127-0078037051920119007503580100-0094009302832029013500700393526200830372012901050297-0092-0086211001460083534004080140007903970110-0088-007902862183-0128007003795419013700670399028101150083-0073226801440066041554950140008103940120-0078-0070027323370126007903685576-0072040501410125007102642412-00750137008103895649008201350384025725010130-0074-0075-0153-008004235752016000740443013500800252256901190357-007600765816-0110-00810332-0087025126390140~007901290082037158810082033901130145271700830093025401360070039559480114007503460150009802622796-0088-0143-007904020271287101550092009901340076038800960281296201600094015200720426
137MARCH1963EDWARDN.LORENZputing machine. Approximately one second per itera-tion, aside from output time, is required.For initial conditions we have chosen a slight de-parture from the state of no convection, namely (O,1,O).Table 1 has been prepared by the computer. It gives thevalues of V (the number of iterations),X, Y, and Z atevervfifthiterationforthefirst160iterations.Intheprinted output (but not in the computations) the valuesof X, Y, and Z are multiplied by ten, and thenonlythose figures to the left of the decimal point are printed.Thus the states of steady convection would appear as0084,0084,0270and-0084,-0084,0270,whilethestate of no convection would appear as 0000, 0000, 0000.The initial instability of the state of rest is evident. All20three variables grow rapidly, as the sinking cold fluidis replaced by even colder fluid from above, and therisingwarmfluid bywarmerfluidfrombelow,so thatbystep 35 the strength of the convection far exceeds thatFic.1.Numerical solution of the convection equations.Graphof steady convection. Then Y diminishes as the warmof Y as a function of time for the first 1000 iterations (uppercurve), second 1000 iterations (middle curve), and third 1000fluid is carried over the top of the convective cells, soiterations (lowercurve).that by step 50, when X and y have opposite signs,warm fiuid is descending and cold fluid is ascending.Themotion thereupon ceases and reverses its direction, asindicated by the negative values of Xfollowing step 60By step 85 the system has reached a state not far fromthat of steady convection.Between steps 85 and 150 itexecutes a complete oscillation in its intensity, theslight amplification being almost indetectable.The subsequent behavior of the system is illustratedin Fig.1, which shows the behavior of Y for the first3000 iterations.After reaching its early peak near step35 and then approaching equilibrium near step 85, itundergoes systematic amplified oscillations until nearstep 1650. At this point a critical state is reached, andthereafter Y changes sign at seemingly irregular inter-vals,reaching sometimesone,sometimestwo,and some-times three or more extremes of one sign before changingsign again.Fig. 2 shows the projections on the X-y- and Y-Z-planes in phase space of the portion of the trajectorycorresponding to iterations 1400-1900.The states ofsteady convection are denoted by C and C'. The firstportion of the trajectory spirals outward from thevicinity of C,as the oscillations about the state ofsteady convection, which have been occurring since step85,continue togrow.Eventually,near step1650,itcrosses the X-Z-plane, and is then deflected toward theneighborhood of C. It temporarily spirals about C, butcrosses the X-Z-plane after one circuit, and returns tothe neighborhood of C, where it soon joins the spiralover which it has previously traveled. Thereafter itcrosses from one spiral to the other at irregular intervals.Fig.3,in which the coordinates areY and Z,isbasedFrG.2.Numerical solution of the convection _equations.upon the printed values of X, Y, and Z at every fifthProjections on the X-y-plane and the Y-Z-plane in phase spaceiteration for the first 6000 iterations.These values deter-ofthesegmentof thetrajectoryextending from iteration1400toiteration 1900.Numerals"14,."15,"etc.,denotepositionsatmine X as a smooth single-valued function of Y and Ziterations 1400,1500, etc. States of steady convection are denotedover much of the range of Y and Z; they determine Xby C and C
138JOURNALOFTHEATMOSPHERICSCIENCESVOLUME20N400200-200Fic. 3. Isopleths of X as a function of Y and Z (thin solid curves),and isopleths of the lower of twovalues of X, where two values occur (dashed curves), for approximate surfaces formed by all points onlimiting trajectories.Heavy solid curve, and extensions as dotted curves, indicate natural boundaries ofsurfaces.as one of two smooth single-valued functions over theholes near C and C'also represent regions which cannotremainder of the range.In Fig.3 the thin solid lines arebe occupied after they have once been abandoned.Returning to Fig. 2, we find that the trajectory apisopleths of X,and wheretwo values of X exist,thedashed lines are isopleths of the lower value. Thus,parently leaves one spiral only after exceeding somewithin the limits of accuracy of the printed values, thecritical distance from the center. Moreover, the extenttrajectory is confined to a pair of surfaces which appearto which this distance is exceeded appears to determinetomerge in the lower portion of Fig.3.The spiral aboutthe point at which the next spiral is entered; this in turnC lies in the upper surface, while the spiral about C' liesseems to determine the number of circuits tobeexecutedin the lower surface.Thus it is possible for the trajectorybefore changing spirals again.to pass back and forth from one spiral to the otherIt therefore seems that some single feature of a givenwithoutintersecting itself.circuit should predict the same feature of the followingAdditional numerical solutions indicate that othercircuit. A suitable feature of this sort is the maximumtrajectories, originating at points well removed fromvalue of Z, which occurs when a circuit is nearly com-these surfaces, soon meet these surfaces. The surfacespleted. Table 2 has again been prepared by the com-puter, and shows the values of X, Y, and Z at only thosetherefore appear to be composed of all points lying onlimitingtrajectories.iterations N for which Z has a relative maximum. TheBecause the origin represents a steady state, no tra-succession of circuits about C and C'is indicated by thesuccession of positive and negative values of X and Y.jectory can pass through it. However, two trajectoriesEvidently X and Y change signs following a maximumemanate from it, i.e., approach it asymptotically as co.The heavy solid curve in Fig. 3, and its ex-which exceeds some critical value printed as about 385.Fig. 4 has been prepared from Table 2. The abscissatensions as dotted curves, are formed by these two tra-isM,thevalueofthenthmaximumof Z,whilethejectories.Trajectories passing close to the origin willordinate is Mn+t, the value of the following maximum.tend to follow the heavy curve, but will not cross it, sothat the heavy curve forms a natural boundary to theEach point represents a pair of successive values of Ztaken from Table 2. Within the limits of the round-offregion which a trajectory can ultimately occupy. The
139MARCH1963EDWARDN.LORENZThe two three-phase sequences differ qualitatively inin tabulating Z, there is a precise two-to-one relationthat the former possesses two numbers, and the latterbetween Mn and Mn+1.The initial maximum Mi=483only one number,exceeding2.Thus thetrajectory corre-is shown as if it had followed a maximum Mo=385,sponding to the former makes two circuits about C,sincemaxima near385arefollowed bycloseapproachesfollowed by one about C (or vice versa).The trajectoryto the origin, and then by exceptionally large maxima.correspondingtothelattermakesthreecircuits aboutC,It follows that an investigator,unaware of the naturefollowed by threeabout C',so that actuallyonlyZof the governing equations, could formulate an empiricalvaries in three phases, while X and Y vary in sixprediction scheme from the "data""pictured in Figs. 2Now consider a sequence where M。is not a rationaland 4. From the value of the most recent maximum offraction. In this case (36) shows that Mn+ cannot equalZ, values at future maxima may be obtained by repeatedapplications of Fig.4.Values of X,Y, and Z betweenmaxima of Z may be found from Fig. 2, by interpolatingMn+1between neighboring curves. Of course, the accuracy ofpredictions made by this method is limited by theexactness of Figs. 2 and 4, and, as we shall see, by theaccuracy with which the initial values of X, Y, and zare observed.Some of the implications of Fig.4 are revealed by:considering an idealized two-to-one correspondence be-tween successive members of sequences Mo, Mi,..consisting of numbers between zero and one. Thesesequences satisfy therelationsifMn1.The correspondence defined by (35) is shown in Fig. 5,which is an idealization of Fig. 4. It follows from re-peated applications of (35) that in any particularsequence,FrG.4.Corresponding values of relative maximumof Z(abscissa) and subsequent relative maximum of Z (ordinate)M,=mn±2"Mo,(36)occurring during the first 6000 iterations.where mn is an even integer.Consider first a sequence where Mo=w/2p, where u isodd.Inthis caseM,1=,and thesequence terminates.These sequences form a denumerable set,and correspond to the trajectories which score direct hits upon thestateofnoconvection.Next consider a sequence where Mow/2p, where uand are relatively prime odd numbers. Then if k>0,Mp+I+=ws/u, where wx and are relatively prime andwis even.Sincefor any the number of proper fractionsu/ is finite, repetitions must occur, and the sequence0.5His periodic. These sequences also form a denumerableset, and correspond to periodic trajectories.The periodic sequences having a given number ofdistinct values, or phases, are readily tabulated. Inparticular there are a single one-phase, a single two-phase, and two three-phase sequences, namely,2/3, *,2/5, 4/5, ...,o0.52/7, 4/7, 6/7, .,FIG.5.The function Ma+1-2M.if M., serving as an idealization of the locus of points in Fig. 4