Abstract

The goal of this paper is to provide a complete picture of the long-term behavior of Lorenz’s maximum simplification equations along with the corresponding meteorological interpretation for all initial conditions and all values of the parameter.

1. Introduction

In this paper our goal is to analyze the properties of equilibrium points and long-term behavior of the solution of the Lorenz’s (1960) maximum simplification equations (MSE). Lorenz (1960) introduced the maximum simplification equations and demonstrated that they could be used to understand certain atmospheric phenomena, such as the index cycle. The dynamical equations were simplified to a point where reasonable predictions of atmospheric conditions could not be expected, however, certain physical processes and features of interest were retained. This MSE consists of a set of three nonlinear and coupled system of three ordinary differential equations (ODE) that govern the evolution of the amplitudes of three modes arising from the double Fourier series expansion of the solution of the barotropic vorticity equation (BVE). The importance of MSE stems from the fact that while it is the simplest possible (minimal set) approximation to BVE, using it, Lorenz (1960) was able to capture many of the qualitative features of the flow including the exchange of kinetic energy between the zonal and meridional components. Motivated by this success, the basic spirit of this ground-breaking work of Lorenz (1960) namely approximating the solution of complex realistic problems using a low-order model (LOM) has become a standard arsenal in problem solving in geophysical domain. Refer to Saltzman (1962), Lorenz (1963), De Swart (1989), Wiin-Nielson (1992), Gluhovsky et al. (2002), and the references therein for details.

Thompson (1957) used a set of equations similar to MSE to explore the limits of deterministic predictability that stem from the uncertainty in the initial conditions. Epstein (1969) used MSE to study the stochastic dynamic prediction. Recently, Lakshmivarahan et al. (2003) used MSE to compare the performance of first-order and second-order methods for data assimilation. To our knowledge, long-term behavior of Lorenz’s MSE has not been explored. In this paper we examine the growth of small errors superimposed on equilibrium points of MSE leading to a set of interesting heteroclinic orbits. The necessary conditions for obtaining a solution that approaches an equilibrium point are established. In addition, we provide a complete catalog of the behavior of MSE starting from any arbitrary initial point. The impacts of changing various model parameters on the solutions with respect to the index cycle are also examined.

In section 2, MSE and their key properties are derived. The set of all equilibria for MSE are identified and a set of necessary conditions for the orbits to approach an equilibrium point are derived in section 3. A comprehensive analysis of stability of equilibria is covered in section 4. A complete picture of the evolution of small errors superimposed on equilibria is given in section 5. A complete classification of the behavior of the solution starting from an arbitrary point is covered in section 6. Meteorological implications of the above analysis are given in section 7. The concluding section 8 provides a summary.

2. The maximum simplification equations—The model

Our model equations are the Lorenz’s so-called maximum simplification equations given by

 
formula

where A = A(t), F = F(t), and G = G(t) are the time varying amplitudes of the three spectral components obtained from the standard BVE and , , and Ġ refer to their time derivatives. See Lorenz (1960) for details of the derivation. The streamfunction ψ = ψ(x, y, t) and the vorticity ∇2ψ = ∇2ψ(x, y, t) corresponding to (2.1) are given by

 
formula

and

 
formula

where x and y are the usual space coordinates, ∇2 is the standard Laplacian operator, 2π/l = L is the distance between zonal wind maximum, and 2π/k = K is the wavelength of the disturbance. That is, K and L denote the spatial scales of the flow and the disturbance. As an example, L = 10 000 km and K = 5000 km, that is, l = 6.285 × 10−4 and k = 1.257 × 10−3. The expression for the easterly wind u = u(x, y, t) = −∂ψ/∂y and northerly wind υ = υ(x, y, t) = ∂ψ/∂x are given by

 
formula

and

 
formula

a. Dependence on the parameter r

It can be verified that the model (2.1) depends only on the ratio r = k/l = L/K and can be written as

 
formula

where

 
formula
 
formula

and

 
formula

A plot of the values of a = a(r), b = b(r), and c = c(r) as a function of r > 0 is given in Fig. 1. It can be verified that b = c when r = 41/3 = 0.7598 and a = c when r = 43 = 1.3161. Further, b → 0+ (from above) as r → 0 and a → 0 (from below) as r increases.

Fig. 1.

Variation of the coefficients a, b, c of the model in (2.6) with respect to r.

Fig. 1.

Variation of the coefficients a, b, c of the model in (2.6) with respect to r.

Let η = (A, F, G)T ∈ ℝ3 be a column (real) vector of three spectral components, where T denotes the transpose. That is, η1 = A, η2 = F, and η3 = G. Then (2.6) can be succinctly written as

 
formula

where g: ℝ3 → ℝ3 is such that g(η) = [g1(η), g2(η), g3(η)]T with g1(η) = aFG, g2(η) = bAG, and g3(η) = cAF.

It should be interesting to note that the Eq. (2.6) is equivalent to the Euler’s equation of motion for a rigid body and is known as the Euler gyrostat (Obukhov 1973). Derivation and analysis of Euler gyrostat can be found in textbooks on classical mechanics (Landau and Lifshitz 1960; Wittenburg 1977; Arnol’d 1989). Recently Gluhovsky and his collaborators (Gluhovsky and Agee 1997; Gluhovsky and Tong 1999; Gluhovsky et al. 2002) have shown that a Volterra gyrostat (of which the Euler gyrostat is a special case) can be used as a building block for generating a wide variety of low-order models of interest in several application domains.

Implicit in the above description are the two sets of coordinate systems. First is the three-dimensional physical space consisting of the two space variables, x and y, and the time variable t, in which the physical quantities such as the streamfunction ψ, vorticity ∇2ψ, and the wind components u and υ are described. The second is the three-dimensional spectral space in which the spectral amplitudes A, F, and G evolve as functions of time according to (2.6) or (2.10). The goal is to characterize the long-term behavior of the physical quantities from those of the spectral components.

b. Special case of r = 1

It can be verified that a = −½, b = ½ and c = 0 when r = 1. In this case, Ġ = 0 and G(t) = G(0) for all t. Hence, if G(0) ≠ 0, (2.6) reduces to a linear harmonic oscillator

 
formula

whose solution is periodic (Hirsch and Smale 1974). When G(0) = 0 and r = 1, then η(t) ≡ 0. In the following it is tacitly assumed that r ≠ 1.

c. Conservation of kinetic energy and vorticity

Let 𝗗 = diag(d1, d2, d3) be a diagonal matrix where di > 0, i = 1, 2, and 3. Consider the functional f: ℝ3 → ℝ defined by

 
formula

It can be verified that the time derivative of this functional evaluated along the trajectory of (2.6) is given by

 
formula

Thus, given any parameter vector p = (a, b, c)T of the model, for any d = (d1, d2, d3)T that is orthogonal to p, (η) ≡ 0 and for this choice of the vector d or matrix 𝗗, f (η) is conserved by the model. By direct computation it can be verified that the mean kinetic energy

 
formula

and the mean square vorticity

 
formula

are conserved by the model (2.1). That is, E(t) = E(0) and V(t) = V(0) for all t ≥ 0. Refer to Lorenz (1960) and Krishnamurti et al. (1998) for more detail.

d. Solution space and degree of freedom

Let η(0) = [A(0), F(0), G(0)]T be the given initial condition for (2.10). Then (2.15) implies that the solution η(t) of (2.10) is confined to the surface of the V ellipsoid, which is centered at the origin and is given by

 
formula

where

 
formula

Clearly, fV(η) depends only on η(0) and not on r. Similarly, from (2.14) it follows that the solution η(t) of (2.10) is confined to the surface of the E ellipsoid, which is also centered at the origin, and is given by

 
formula

where

 
formula

Substituting (2.19) into (2.18) it can be verified that fE(η) depends only on η(0) and r.

Combining these, it follows that the solution η(t) of (2.10) must lie in the set

 
formula

which is the intersection of the two concentric ellipsoids. That is, this set 𝒮 is the solution space for the model Eq. (2.10). Refer to Fig. 2 for an example of the two intersecting ellipsoids.

Fig. 2.

A view of intersecting ellipsoids and the solution space.

Fig. 2.

A view of intersecting ellipsoids and the solution space.

Taking derivatives of both sides of (2.16) and (2.18) with regard to time, we get

 
formula

where fV(η) = (1/V(0))(A, F, 2G)T is the gradient of fV(η) and

 
formula

where fE(η) = [1/2E(0)][(A/l2), (F/k2), (2G/l2 + k2)]T is the gradient of fE(η). It can be verified that the gradient vectors fV and fE are linearly independent unless r = 1 and G = 0. Thus, the flow of the dynamics in (2.10) is always orthogonal to the plane defined by these two gradient vectors. Since fV and fE define the outward normal to the surface of the ellipsoids, the flow always lies along the tangent to the points in 𝒮 in (2.20).

Lorenz (1960) has shown that (2.6) can be solved analytically where the solution A(t), F(t), and G(t) are expressed as elliptic functions of time. In principle, we can use the tables of elliptic functions (Abramowitz and Stegun 1965) to characterize the long-term behavior of the solution of (2.6). However, following Lorenz (1960), we use the standard numerical methods to solve (2.6). While Lorenz (1960) used the central difference scheme that involves three time steps, we use the fourth-order adaptive Runge–Kutta (RK) method.

e. Boundedness of the trajectories

An immediate consequence of the property described in section 2b is that the trajectories η(t) of (2.10) remains bounded. Let η(0) = [A(0), F(0), G(0)]T be the initial state. Then from (2.14) and (2.15) we get

 
formula

For example, when A(0) = 10−3 = G(0) and F(0) = 5, K = 5000 and L = 10 000, it can be verified that |A(t)| ≤ 2.5, |F(t)| ≤ 5.0, and |G(t)| ≤ 3.5355, for all t ≥ 0.

This property in turn implies that the model (2.1) is (externally) stable in the sense of Lyapunov (Hirsch and Smale 1974; Sparrow 1982). However, in view of the nonlinear interaction between modes, there is transfer of energy between modes leading to (internal) instability.

An example illustrating this transfer of energy between modes is given in Fig. 3 where the system (2.6) is solved using the fourth-order adaptive Runge–Kutta method. This figure provides a complete picture of the model trajectories, namely the variation of A(t), F(t), and G(t) with time; two dimensional plots of A(t) − F(t), A(t) − G(t) and F(t) − G(t), and the three-dimensional view of the solution.

Fig. 3.

Solution of (2.6) with initial condition η(0) = (0.001, 1.0, 0.001)T, K = 5000, L = 10 000 (hence r = 2).

Fig. 3.

Solution of (2.6) with initial condition η(0) = (0.001, 1.0, 0.001)T, K = 5000, L = 10 000 (hence r = 2).

3. Equilibrium points of the model

Recall that if = g(η) represents a dynamical system, then the set of points at which g(η) = 0 are called the equilibrium points of the system (Hirsch and Smale 1974; Devaney 1989; Sparrow 1982). Accordingly, we define three sets of equilibria for the system in (2.6). Let

 
formula

It can be verified that ɛA, ɛF and ɛG constitute the set of all equilibrium points of (2.6). It can be seen that the origin is the only point common to all these three sets.

Given that the solution of the model (2.10) always lies in the set 𝒮 in (2.20), a question of interest at this point is: Which of the above equilibrium points lie arbitrarily close to the solution space 𝒮? Stated in other words, when does the model solution starting from an arbitrary initial point ever approach an equilibrium point?

In answering this question, first consider η* = (A*, 0, 0)TɛA. It can be verified that η* ∈ 𝒮𝒜 when

 
formula

that is, when

 
formula

in which case A* = ±2V = ±2lE. Thus, (3.2) is a necessary condition for an orbit to include η* ∈ ɛA.

Considering η* = (0, F*, 0)TɛF and repeating the above argument, it can be verified that

 
formula

is a necessary condition for an orbit to include η* ∈ ɛF where F* = ±2V = ±2kE. Similarly, considering η* = (0, 0, G*)TɛG, it readily follows that

 
formula

is necessary for the orbit to include η* ∈ ɛG, where G* = ±2V = ±2(l2 + k2)E.

Stated in other words, the key to understanding the behavior of orbits of the model in (2.10) is the value of the ratio V/E. Lorenz (1960) recognized the importance of V/E (= 2k2 as the threshold for sub and supercritical disturbances, for example), but did not explain why it was significant. We see here that the necessary conditions for obtaining any solution including an equilibrium point as part of the orbit are set by specific values of the ratio V/E.

4. Stability of equilibria

Recall that the local stability properties of an equilibrium point, say η*, is governed by the distribution of the eigenvalues of the Jacobian, 𝗗g(η) of the mapping g(η) in (2.10) evaluated at η* (Hirsch and Smale 1974; Devaney 1989; Lewis et al. 2006, chapter 32). The Jacobian of g(η) is given by

 
formula

While the model (2.6) is symmetric with respect to the variables A, F, and G, since the coefficients a, b, and c are different, we consider three cases (refer to Fig. 1).

a. Stability of equilibrium points in ɛA

Let η* ∈ ɛA. Then, the Jacobian (4.1) becomes

 
formula

Let r < 1. Then, a < 0, b > 0, and c > 0. It can be verified that the eigenvalues of 𝗗g(η*) are given by

 
formula

and the corresponding eigenvectors are given by

 
formula

For latter reference, define the diagonal matrix Λ(1) of eigenvalues and 𝗩(1), the nonsingular matrix of the corresponding eigenvectors as

 
formula

and

 
formula

where (it follows from the definition)

 
formula

Notice that while the eigenvalues λ(1)2 and λ(1)3 in (4.3) depend on A and the ratio r, the eigenvectors v(1)2 and v(1)3 are independent of A. That is, every point in ɛA shares the same set of eigenvectors that depends only on the value of the ratio r. For purposes of comparison, the eigenvalues and vectors of 𝗗g(η*) when η* belongs to each of the three equilibria are summarized in Table 1.

Table 1.

Comparison of the eigenvalues and eigenvectors of Jacobian at the equilibria

Comparison of the eigenvalues and eigenvectors of Jacobian at the equilibria
Comparison of the eigenvalues and eigenvectors of Jacobian at the equilibria

Let r > 1. Then a < 0, b > 0, and c < 0. It can be verified that the complex eigenvalues of 𝗗g(η*) are given by

 
formula

and the corresponding complex eigenvectors are v(1)1, w(1)2 = u1 + iu2 and w(1)3 = u1iu2, respectively, where

 
formula

We leave it to the reader to verify that this eigen decomposition in the complex domain can be equivalently expressed in real domain using the standard canonical form (Hirsch and Smale 1974) as

 
formula

where

 
formula

and

 
formula

That is, the eigenvectors that are the columns of 𝗩(1) are given in (4.8). In this paper we use this latter form of eigen decomposition when the eigenvalues are complex. Refer to Table 1 for a summary.

The product bc and the ratio c/b in (4.3) and (4.4) are given by

 
formula

and

 
formula

Variation of bc and c/b as a function of r are given in Figs. 4a,b, respectively. It follows from (4.3)(4.4) and (4.8)(4.9) that the eigenvalues and eigenvectors depend on two factors: the value of the first coordinate A of the equilibrium point η* = (A, 0, 0)T and the value of the ratio r. For completeness, we consider several subcases.

Fig. 4.

Variation of bc and c/b with respect to r.

Fig. 4.

Variation of bc and c/b with respect to r.

1) Case 1: Let A > 0 and r < 1

In this case, it can be verified that λ(1)2 > 0 and λ(1)3 < 0, and the equilibrium point, η* = (A, 0, 0) is called a saddle point, with the eigenvector v(1)2 corresponding to the unstable (repelling) direction and the eigenvector v(1)3 to the stable (attracting) direction (Hirsch and Smale 1974). Refer to Table 2 for a summary. In verifying these claims, let η = η* + ξ where ξ is a small perturbation. Then, using the first-order Taylor series expansion, we obtain

 
formula

or

 
formula

as the local dynamics of evolution of the perturbation. Now change the variables from the standard (A, F, G) orthogonal coordinate system to oblique coordinate system defined by the eigenvectors v(1)1, v(1)2 and v(1)3 in (4.6). Define a new vector z = (z1, z2, z3)T where

 
formula

Substituting (4.13) into (4.12) and using (4.7) we get

 
formula

That is,

 
formula

and z2(t) = z2(0)eλ(1)t2 and z3(t) = z3(0)eλ(1)t3. Since λ(1)2 > 0, indeed z2 grows exponentially and λ(1)3 < 0 implies that z3 decreases exponentially as well. Refer to Fig. 5 for a geometric exposition of these ideas. Thus, any perturbation with a nonzero component along v(1)2 will grow. Hence, in this case η* = (A, 0, 0) is an unstable equilibrium point. Referring to Fig. 5, the angle θ that v(1)2 makes with the positive F axis (same as the angle v(1)3 makes with the negative F axis) is given by

 
formula

Thus, as r increases from 0 to 1, θ decreases from 90° to 0°.

Table 2.

Properties of eigenvalues of 𝗗g(η*), η* ∈ ɛA

Properties of eigenvalues of 𝗗g(η*), η* ∈ ɛA
Properties of eigenvalues of 𝗗g(η*), η* ∈ ɛA
Fig. 5.

A geometric view of the eigenvectors v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/a)]T at an equilibrium point η* = (A, 0, 0) with (a) A = 1 and r = 0.5 [hence b = 0.1, c = 0.75, (c/b) = 2.7386] and (b) A = −1, r = 0.5; v(1)1 is perpendicular to the plane containing v(1)2 and v(1)3. The equilibrium point is a saddle point.

Fig. 5.

A geometric view of the eigenvectors v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/a)]T at an equilibrium point η* = (A, 0, 0) with (a) A = 1 and r = 0.5 [hence b = 0.1, c = 0.75, (c/b) = 2.7386] and (b) A = −1, r = 0.5; v(1)1 is perpendicular to the plane containing v(1)2 and v(1)3. The equilibrium point is a saddle point.

2) Case 2: Let A < 0 and r < 1

In this case since λ(1)2 < 0 and λ(1)3 > 0, the role of the eigenvalue–eigenvector pairs (λ(1)2, v(1)2) and (λ(1)3, v(1)3) are interchanged and the above analysis and conclusions hold. Thus, the equilibrium point η* = (A, 0, 0) is an unstable saddle point. Refer to Fig. 5.

(i) Effect of changing A

Thus, keeping r(<1) fixed, the force of attraction along v(1)3 and the force of repulsion along v(1)2 as measured by the magnitudes of the eigenvalues |λ(1)3| = Abc and |λ(1)2| = Abc, respectively, diminish to zero as A decreases from 1 to 0. As A increases in value in the negative direction, while the magnitudes of these forces increase, their role is reversed, with v(1)2 being the new direction of attraction and v(1)3 being the new direction of repulsion. This reversal of the role between v(1)2 and v(1)3 occurs when A = 0.

(ii) Effect of changing r

Let A > 0 be fixed. From (4.15) it follows that as r → 0+, c/b → ∞, and θ → 90°. Combining this with (4.4) we get a degenerate limit, namely

 
formula

As r increases form 0+, the angle θ decreases from 90°. When r → 1, (c/b) → 0+ and θ → 0+. Again we get a degenerate limit

 
formula

When r = 1, c = 0 and the character of the equilibrium η* = (A, 0, 0)T changes from being a saddle point to a center. Refer to the next subsection for more details.

3) Case 3: Let A > 0 and r > 1 or A < 0 and r > 1

From (4.8) it follows that the nonzero eigenvalues are complex conjugate with zero real part. In this case the equilibrium point is called a center for the corresponding linearized problem (Hirsch and Smale 1974). A geometric view of the two new eigenvectors in (4.9) is shown in Fig. 6. It turns out that while the new eigenvector v(1)3 = (0, 1, 1)T is independent of r and makes a fixed 45° angle with the positive F axis, the new eigenvector v(1)2 = (0, −|(b/c)|, |(c/b)|)T makes an angle θ with the negative F axis where

 
formula

As r increases from one to infinity, |c/b| increases form 0 to ½ and θ from 0 to 26.57°.

Fig. 6.

A geometric view of the eigenvectors v(1)3 = (0, 1, 1)T and v(1)2 = [0, −|(b/c)|, |(c/b)|]T at the equilibrium point η* = (A, 0, 0)T with A = 1 and r = 2.0 [hence b = 1.6, c = −0.75 |(c/b)| = 0.469 and |(b/c)| = 1.46]; v(1)1 is perpendicular to the plane containing v(1)2 and v(1)3. The equilibrium point is a center.

Fig. 6.

A geometric view of the eigenvectors v(1)3 = (0, 1, 1)T and v(1)2 = [0, −|(b/c)|, |(c/b)|]T at the equilibrium point η* = (A, 0, 0)T with A = 1 and r = 2.0 [hence b = 1.6, c = −0.75 |(c/b)| = 0.469 and |(b/c)| = 1.46]; v(1)1 is perpendicular to the plane containing v(1)2 and v(1)3. The equilibrium point is a center.

(i) Change of character of the equilibrium η* = (A, 0, 0)T

As r changes from 1 to 1+, it follows from (4.17) that the old eigenvectors v(1)2 and v(1)3 in (4.4) become degenerate and cease to exist. Simultaneously, a new set of eigenvectors v(1)2 and v(1)3 in (4.9) are created. Of these, v(1)3 = (0, 1, 1)T is fixed. When r → 1+, |c/b| → 0 and hence new v(1)2 → (0, −1, 0)T. Comparing this with (4.17), it looks as though the old v(1)3 is reincarnated as the new v(1)2 as r changes from 1to 1+. Thus, as r changes from 1 to 1+, η* = (A, 0, 0)T changes its character from being a saddle to a center.

4) Case 4: Either A = 0 or r = 1

In this case, λ(1)2 = 0 = λ(1)3.

b. Stability of equilibrium points in ɛF

Let η* = (0, F, 0)TɛF. Then, the Jacobian in (4.1) becomes

 
formula

whose eigenvalues and eigenvectors are given by

 
formula

with λ(2)1 = 0, λ(2)2 = Fac, λ(2)3 = −Fac and

 
formula

The product ac and the ratio c/a are given by

 
formula

and

 
formula

Variation of ac and c/a as a function of r are given in Figs. 7a,b, respectively. From (4.20) and (4.22) it follows that the eigenvalues λ(2)2 and λ(2)3 depend on two factors: the value of the second coordinate F of the equilibrium point η* = (0, F, 0)T and the value of the ratio r. A summary of the properties of these eigenvalues for various combinations of values of F and r are given in Table 3.

Fig. 7.

Variation of ac and c/a with respect to r.

Fig. 7.

Variation of ac and c/a with respect to r.

Table 3.

Properties of eigenvalues of the equilibrium, η* = (0, F, 0)TɛF

Properties of eigenvalues of the equilibrium, η* = (0, F, 0)T ∈ ɛF
Properties of eigenvalues of the equilibrium, η* = (0, F, 0)T ∈ ɛF

Again it can be verified that η* = (0, F, 0)T is an unstable saddle point when r > 1 and F ≠ 0. However, when r < 1 and F ≠ 0, perturbations around this equilibrium will have a local periodic behavior. It is interesting to note that the behavior of the equilibria in ɛA for r < 1 are quite similar to those in ɛF for r > 1 and vice versa.

c. Stability of equilibrium points in ɛG

Let η* = (0, 0, G)TɛG. Then the Jacobian in (4.1) becomes

 
formula

whose complex eigenvalues are given by λ(3)1 = 0, λ(3)2 = iG|ab|, λ(3)3 = −iG|ab|. Using the standard canonical form, the eigen decomposition becomes

 
formula

where

 
formula

and

 
formula

The product ab and ratio b/a are given by

 
formula
 
formula

From (4.25) it follows that the eigenvalues λ(3)2 and λ(3)3 are always complex conjugates (with zero real parts) for all values of r. Variation of ab and b/a as a function of r are given in Figs. 8a,b, respectively. Refer to Table 4 for a summary of the properties of these eigenvalues.

Fig. 8.

Variation of ab and b/a with respect to r.

Fig. 8.

Variation of ab and b/a with respect to r.

Table 4.

Properties of eigenvalues of the equilibrium, η* = (0, 0, G)TɛG

Properties of eigenvalues of the equilibrium, η* = (0, 0, G)T ∈ ɛG
Properties of eigenvalues of the equilibrium, η* = (0, 0, G)T ∈ ɛG

5. Small perturbation of equilibrium

In this section our goal is to quantify the growth of small errors superimposed on equilibrium points.

Let η* be an equilibrium point and ξ ∈ ℝ3 denote a (small) perturbation or error vector. Let η = η* + ξ be the perturbed initial condition for the model (2.6). Since * = g(η*) = 0 by definition, we get

 
formula

as the exact dynamics of the perturbation. Following the developments in section 4, we consider three cases, one for each of the three sets of equilibria.

a. Small perturbation of equilibrium points in ɛA

Let η* = (A, 0, 0)TɛA. Depending on the values of A and r, consistent with section 4, we again consider four subcases.

1) Case 1: Heteroclinic orbits

Let A > 0 and r < 1. In this case, η* = (A, 0, 0)T is a saddle point with v(1)2 = (0, 1, c/b)T denoting the local unstable direction and v(1)3 = (0, 1, −c/b)T denoting the local stable direction. Refer to Tables 1 and 2 for details.

These eigenvectors v(1)2 and v(1)3 partition the FG plane (where A is constant) into four regions R1 to R4. Refer to Fig. 9 for details. It turns out that there is an intrinsic difference in the qualitative behavior of the trajectories starting from each of these four regions but close to the equilibrium point η*.

Fig. 9.

The origin is the equilibrium point η* = (A, 0, 0)T with A = 1. A decomposition of the FG plane (A is constant) by v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/b)]T into four regions ℛ1, ℛ 2, ℛ 3, and ℛ 4. In this figure, r = 0.5, b = 0.1, c = 0.75, and (c/b) = 2.7386.

Fig. 9.

The origin is the equilibrium point η* = (A, 0, 0)T with A = 1. A decomposition of the FG plane (A is constant) by v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/b)]T into four regions ℛ1, ℛ 2, ℛ 3, and ℛ 4. In this figure, r = 0.5, b = 0.1, c = 0.75, and (c/b) = 2.7386.

For concreteness, let A = 1.0, r = 0.5, and ξ(0) = (0.0, 0.001, 0.0)T ∈ ℛ1. Then, η* = (1, 0, 0)T and η = (1.0, 0.001, 0.0)T. A complete picture of the evolution of the model Eq. (5.1) starting from η is given in Fig. 10. A number of observations are in order.

Fig. 10.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 10.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

(i) Heteroclinic orbits

The orbit or the trajectory starting from η ∈ ℛ1 that is close to the equilibrium point η* = (1, 0, 0)T moves away from it along the unstable direction v(1)2. In due course, this orbit bends with A decreasing but F remaining positive. It soon approaches the equilibrium point η** = (−1, 0, 0)T, which is the symmetric dual of η*, along v(1)2, which is the local direction of attraction for η**. As this trajectory gets very close to η**, it gets repelled along v(1)3, which is the local unstable direction for η**. Again it bends with A increasing but F still remaining positive. It now approaches η* along v(1)3, which is its direction of attraction. As it gets closer to η*, it again gets repelled along v(1)2 and this cycle repeats. Figure 10 shows the results of integration over a time period T = 200 units using the adaptive fourth-order Runge–Kutta scheme.

The class of orbits so generated is called a heteroclinic orbit (Sparrow 1982). A distinguishing feature of this orbit is that one part of it connecting η* and η** lies at the intersection of unstable manifold of η* and stable manifold of η**. The other part of the orbit, from η** to η*, lies at the intersection of the unstable manifold of η** and the stable manifold of η*.

(ii) Four distinct heteroclinic orbits

It turns out that each equilibrium point η* ∈ ɛA is endowed with four distinct classes heteroclinic orbits—one for each of the four regions ℛ1 to ℛ4. Figures 11 –13 contain illustrations of heteroclinic orbits starting from initial points in regions ℛ2 to ℛ4, respectively. The key distinguishing features of these four heteroclinic orbits are summarized in Table 5.

Fig. 11.

Heteroclinic orbit starting from η in region ℛ2; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.0, 0.001)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 11.

Heteroclinic orbit starting from η in region ℛ2; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.0, 0.001)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 13.

Heteroclinic orbit starting from η in region ℛ4; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.0, −0.001)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 13.

Heteroclinic orbit starting from η in region ℛ4; η* = (A, 0, 0)T, A = 1, ξ = (0.0, 0.0, −0.001)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Table 5.

Properties of heteroclinic orbits: A = 1.0, r = 0.5 and η* = (1, 0, 0)T.

Properties of heteroclinic orbits: A = 1.0, r = 0.5 and η* = (1, 0, 0)T.
Properties of heteroclinic orbits: A = 1.0, r = 0.5 and η* = (1, 0, 0)T.
(iii) Effect of changing A

From Tables 1 and 2, it is clear that eigenvectors v(1)2 and v(1)3 of 𝗗g(η*) for η* ∈ ɛA are independent of A. But their associated eigenvalues λ(1)2 = Abc and λ(1)3 = −Abc, which are measures of the force of repulsion and attraction are directly proportional to A. Thus, as A increases in magnitude, the resulting increases in the force of repulsion and attraction, while do not alter the basic properties of the heteroclinic orbits, do affect the frequency of oscillations of the components of η. Examples of heteroclinic orbits for A = 5 and 10 with ξ = (0.0, 0.001, 0.0)T and r = 0.5 are given in Figs. 14 and 15.

Fig. 14.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 5, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 14.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 5, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 15.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 10, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 15.

Heteroclinic orbit starting from η in region ℛ1; η* = (A, 0, 0)T, A = 10, ξ = (0.0, 0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

(iv) Effect of changing r

From Fig. 4, it follows that c/b is a decreasing function of r with c/b being positive for 0 < r < 1 and negative for r > 1. We restrict this discussion to the case when 0 < r < 1 (r > 1 is considered in case 3 below). Referring to Fig. 9, it follows that (c/b) determines the slope of the eigenvectors v(1)2 and v(1)3 in the FG plane. To understand the impact r, consider the example of the heteroclinic orbit in Fig. 10. In this example, since F remains positive, combining it Fig. 9, it can be verified that the heteroclinic orbit approaches η* = (1, 0, 0)T along v(1)3 (that separates region ℛ1 and ℛ4) and leaves of η* along v(1)2 (that separates region ℛ1 and ℛ2). That is, this orbit makes a turn (to the right) by an angle ϕ = 180 − 2θ where θ is given by (4.15). Clearly, this angle ϕ is small or large depending on if r is small or large.

(v) Effect of changing L and K

While the solution (2.6) depends only on the ratio r = L/K, from (2.2)(2.5) it follows that L and K affect the values and shape of the streamfunction ψ, vorticity ∇2ψ, and the wind components u and υ.

2) Case 2: A < 0 and r < 1

In this case, the equilibrium point η* = (−A, 0, 0)TɛA is again a saddle point, with v(1)2 being the stable direction and v(1)3 the unstable direction. Hence all the comments and conclusions of case 1 carry over to this case.

3) Case 3: A > 0 and r > 1 or A < 0 and r > 1

In this case b > 0 but c < 0 and recall from Table 1 that the eigenvalues λ(1)2 = iA|bc| and λ(1)3 = −iA|bc|. This in turn implies that η* = (A, 0, 0)T is a center and the components F(t) and G(t) oscillate where the frequency of oscillation depends on both A and r. In verifying this, assume that A is nearly constant. In this case (2.6) becomes

 
formula

Referring to Hirsch and Smale (1974), we get the solution of (5.2) as

 
formula
 
formula

The period T of the harmonic components in (5.3)(5.4) is given by

 
formula

Thus, T is a decreasing function of A and r.

The question now is: when will A be nearly a constant? Recall that when r is large, a < 0 but close to zero (a → 0). Hence is small in magnitude but its sign depends on those of F and G. Refer to Fig. 16. Thus for large r, A is very slowly varying and the orbits of (2.6) look quite similar to a merry-go-round where the cars go round and round in closed circuits but they reach a maximum when they cross the G axis and a minimum when they cross the F axis. Refer to Fig. 17 for an example of the trajectory in this case.

Fig. 16.

Signs of when r is large and a → 0.

Fig. 16.

Signs of when r is large and a → 0.

Fig. 17.

Heteroclinic orbit starting from η = (1.0, 0.01, 0.01)T, L = 10 000, K = 500, r = 20, a = −0.1247 × 10−3, b = 19.95, and c = −9.975.

Fig. 17.

Heteroclinic orbit starting from η = (1.0, 0.01, 0.01)T, L = 10 000, K = 500, r = 20, a = −0.1247 × 10−3, b = 19.95, and c = −9.975.

4) Case 4: A = 0 and r = 1

The special case of r = 1 is considered in detail in section 2. Consider the case when A = 0. In this case η* = (0, 0, 0)T and all the three eigenvalues of 𝗗g(η) are zero. Figure 18 provides examples of the evolution of perturbations for this case.

Fig. 18.

Orbit starting from η = (0.001, 0.001, 0.001)T, where L = 10 000, K = 500, r = 2.0.

Fig. 18.

Orbit starting from η = (0.001, 0.001, 0.001)T, where L = 10 000, K = 500, r = 2.0.

b. Small perturbation of equilibrium points in ɛF

Comparing entries in Table 1, it turns out that there is a one to one correspondence between the stability properties of equilibria in ɛA and ɛF. In particular, the properties of the eigenvalues and eigenvectors of the Jacobian of equilibria in ɛF for r > 1 (r < 1) are very similar to those of ɛA for r < 1 (r > 1). This in turn implies that the qualitative properties of the trajectories starting close to ɛF when r > 1 (r < 1) are quite similar to those starting close to ɛA when r < 1 (r > 1). In view of this similarity, to save space, we observe that all the remarks and conclusion of section 5a carry over to equilibria in ɛF.

c. Small perturbations of equilibrium points in ɛG

From (2.7) and (2.8) we have a < 0 and b > 0 for all r > 0. In view of this, from Table 1 it turns out that the eigenvalues of the Jacobian at the equilibrium point in ɛG are purely imaginary for all r > 0. This in turn implies that the qualitative properties of trajectories starting close to ɛG are quite similar to those starting close to ɛA when r > 1 or ɛF for r < 1. Hence all the relevant conclusions from section 5a carry over to this case as well.

6. Properties of orbits starting from an arbitrary point

In this section we summarize the qualitative properties of orbits starting from an arbitrary (not necessarily close to an equilibrium) point. Based on the conditions (3.2)(3.4), we consider three cases.

a. Case 1

In view of (3.2), from (2.14) and (2.15) consider the ratio

 
formula

This ratio depends on the initial point η = (A, F, G)T and the value of r. Following Lorenz (1960) we define three types of initial conditions: R1(l) = 1 is called critical and R1(l) < (>)1 is called sub (sup) critical cases. From (6.1) it can be verified that

 
formula

Let

 
formula

Since G2/F2 is nonnegative, we consider two subcases.

1) 0 < r < 1

In this case β is positive. Hence, 0 < R1(l) < 1 for all −β < δ < 0 and R1(l) > 1 when δ > 0 with R1(l) = 1 when δ = 0. As an example, if r = ½, then β = 15. Thus, R1(l) is subcritical for all −15 < δ < 0, and supercritical for all δ > 0.

The cross sections of the V and E ellipsoids by the condition (6.3) are given by

 
formula

and

 
formula

Combining these with (6.1) we get

 
formula

Thus, under the condition (6.3), R1(l) = 1 when δ = 0. To check if an equilibrium point in ɛA is reachable, let G = 0 in (6.3). Then F = 0 and from (6.4) and (6.5) we get

 
formula

which from (6.6) is consistent only when δ = 0. This argument leads to the following inescapable conclusion: orbits starting from initial points lying outside of the pair of straight lines given by

 
formula

cannot approach points in ɛA.

A number of observations are in order.

  1. Critical points and eigenvectors of 𝗗g(η*), η* ∈ ɛA. Combining (4.25) with (6.8), the latter can be written as 
    formula
    Thus, the pair of straight lines given by G = ±(c/b)F in the FG plane (refer to Fig. 9) indeed corresponds to the eigenvectors v(1)2 = [0, 1, (c/b)]T and v(1)3 = [0, 1, −(c/b)]T. Thus, points along these two eigenvectors in the FG plane in Fig. 9 correspond to the critical case. Similarly, the region ℛ1 and ℛ3 correspond to the subcritical and region ℛ2 and ℛ4 to the supercritical cases.
  2. Critical points and heteroclinic orbits. Since the orbits starting from points along v(1)2 and v(1)3 approach points in ɛA, these orbits are necessarily heteroclinic orbits corresponding to the unstable case in section 5a. However, orbits starting from the sub and super critical regions can never approach an equilibrium point in ɛA.

2) r > 1

Since

 
formula

it follows from (6.1) that

 
formula

Also, since β < 0, from (6.3) it follows that δ ≥ |β| > 0 for 2G2/F2 to be nonnegative. Hence this corresponds to the supercritical case and in this case orbits starting from any η = (A, F, G)T ∈ ℝ3 cannot reach an equilibrium point in ɛA. This is to be anticipated since when r > 1, any equilibrium point in ɛA is a center.

b. Case 2

In view of (3.3) using (2.14) and (2.15), consider the ratio

 
formula

Analogous to (6.2) we readily see that

 
formula

Let

 
formula

Since the lhs of (6.13) is nonnegative, we consider two subcases.

1) r > 1

In this case, since α is positive, we get 0 < R2(k) < 1 (subcritical) when −α < δ < 0 and R2(k) > 1 (supercritical) when δ > 0. In the critical case (δ = 0), using (4.23), we obtain

 
formula

Thus, the pair of straight lines corresponding to (6.14) are given by G = ±(c/a)A. These straight lines correspond to the eigenvectors v(2)2 = [1, 0, (c/a)]T and v(3)2 = [1, 0, −(c/a)]T. Thus, the orbits starting from points along v(2)2 and v(2)3 approach the equilibrium points η* = (0, F*, 0)T where F* = ±2υ = ±2kE. As in case 1, orbits starting from sub (super) critical regions cannot approach an equilibrium points in ɛF.

2) r < 1

Since

 
formula

it follows that R2(k) > 1 for all 0 < r < 1. Hence this corresponds to the supercritical case and orbits starting from any initial point η ∈ ℝ3 cannot approach an equilibrium point in ɛF. Again this is to be anticipated since η* ∈ ɛF is a center when r < 1.

c. Case 3

From (3.3), (2.14), and (2.15), it follows that

 
formula

Since

 
formula

for all r, R3(l, k) < 1. Again it can be verified by following the above arguments that equilibrium η* ∈ ɛG can never be reached by orbits starting from any point η ∈ ℝ3.

7. Index cycle

Lorenz (1960) demonstrated that the set of three ODEs in (2.1) were capable of portraying the tendency of the average strength of the midlatitude westerly winds to undergo considerable variations in intensity, a phenomenon known as the index cycle (Namias 1950). The index cycle represents a cyclical transformation of the general circulation from a pattern of strong, zonal flow in the midlatitudes along with weak subtropical flow (high index) to strong westerly flow in the subtropical latitudes combined with weak midlatitude flow (low index) and back again, a variation lasting over an extended period of several weeks. The first term in the u component of the wind (2.4) is independent of x, representing the zonal flow. The latitudes of the zonal wind maxima are fixed (y = 3π/2l, 7π/2l, etc.) and the intensity (A/l) varies, therefore Lorenz (1960) noted that the variable A could represent the zonal index. Large positive values of A correspond to strong west-to-east flow at locations of u-wind maxima, while large negative values of A result in a shift in the location of the maximum zonal flow to the south by a distance π/l. An index cycle in this model will therefore be indicated by the value of A changing from Amax to −Amax and back to Amax. Physically plausible values of A are determined by assignment of an arbitrary standard unit of vorticity, for example, Lorenz (1960) set 1 unit of vorticity = 1/3 h−1. Given this, a value of A = 0.5 would correspond to a maximum zonal wind of approximately 75 m s−1.

a. Index cycles associated with heteroclinic orbits

In this work, we define the start of an index cycle to be the time when A first drops below 99% of Amax and the end of the cycle to be when A rises above 99% of Amax. The period of the model index cycle is defined as the length of time between the start and end of the cycle. Given this definition, we can examine the sensitivity of the index cycle period to changing initial conditions and values of r. Since we use the variable A to define the index cycle, we focus our analysis around the A equilibrium. The results of this sensitivity analysis are summarized in Tables 6 –8 

Table 6.

Variation of index cycle with respect to A(0).

Variation of index cycle with respect to A(0).
Variation of index cycle with respect to A(0).
Table 8.

Variation of index cycle with respect to initial perturbation.

Variation of index cycle with respect to initial perturbation.
Variation of index cycle with respect to initial perturbation.

An example of an index cycle portrayed by this model is provided in Fig. 10, given by a small F perturbation imposed upon ɛA. As noted in Fig. 9 this is a perturbation in the ℛ1 region of the FG domain, and results in a heteroclinic orbit with properties summarized in Table 5. For several units of time, A remains near its initial value (A = +1) then at approximately t = 15, the F and G variables quickly increase while the value of A rapidly switches from +1 to −1. As A approaches −1, the values of F and G decrease to zero. The solution remains nearly the same for several units of time (around t = 50). The solution suddenly shifts again at approximately t = 85, when F increases and G decreases rapidly, and the value of A switches from −1 to +1. The period of the index cycle in this case is 69.8 units of time.

A similar example is provided in Fig. 11, given by a small G perturbation imposed upon ɛA. This is a perturbation in the ℛ2 region of the FG domain, which results in a heteroclinic orbit with properties described in Table 5. As with the ℛ1 example, the evolution of A is consistent with the index cycle concept: A remains near its initial value (A = +1) for several units of time, switches quickly from +1 to −1, remains near the minimum value for some time, then quickly changes again from −1 to +1. However, the behavior of the F evolution in this example is similar in character to the G evolution in the previous (ℛ1) example, F changes sign as A switches from +1 to −1 and back again, while G remains positive throughout the evolution. In this case, the period of the index cycle is slightly longer (77.1 units of time) than in the previous example.

1) Effect of changing A

In the unstable situation (r < 1) near ɛA, the index cycle period decreases as the initial value of A increases (Table 6). It appears that relatively strong zonal flow (high initial A) will result in a shorter index cycle. This makes physical sense, in the cases with stronger initial zonal flow, more energy will be available for transfer to the perturbations, allowing for a faster transformation of the flow. This is also consistent with the stability analysis of ɛA described in section 4a(2). Larger values of A produce a larger force of attraction/repulsion about the saddle point, which is consistent with a short index cycle period.

2) Effect of changing r

The sensitivity of the index cycle period to variations in r is summarized in Table 7. Similarly to the previous example, for a given value of A, one expects the minimum index cycle period to occur when the magnitude of the eigenvalue is maximized. Since λ(1)2 = Abc, this would occur when the product bc is maximized. As shown in (4.10), the product bc is maximized for r = 0.6436. This is confirmed experimentally in Table 7 for values of r < 0.6 the index cycle period increases with decreasing r, and for values of r > 0.6 the index cycle increases with increasing r.

Table 7.

Variation of index cycle with respect to r

Variation of index cycle with respect to r
Variation of index cycle with respect to r

3) Effect of changing perturbation size

The size of the initial perturbation also affects the length of the index cycle (Table 8). Larger perturbations result in shorter index cycles. Since the evolution of F in the ℛ1 class of solutions is characterized as an oscillation between the small initial value and some much larger maximum, the initial perturbation fixes the minimum value of F. In order for A to remain near its maximum/minimum value (as in a long period index cycle), F and G must remain near zero for a relatively long time. This will not happen if the initial F perturbation is not near zero. The same is true for G perturbations (ℛ2 class of solutions). It is interesting to note that, in the unstable case, identical index periods are obtained for G perturbation =(c/b) × F perturbation (for r = 0.5, (c/b) = 2.7386).

b. Index cycles in the stable case

Lorenz (1960) stated that cycling of A was possible in the stable (r > 1) situation if the perturbations were large and V/E was supercritical. As shown in section 6, the critical point for r > 1 is determined by G2/A2 = (r4 − 1)/2. However, in both the subcritical and supercritical cases, the evolution of A is not consistent with the definition of an index cycle. An example of a subcritical case is shown in Fig. 3. While A varies considerably in this example, it remains positive and does not remain near its maximum or minimum value for any considerable period of time. On the other hand, the F evolution is consistent with the index cycle definition. In the supercritical situation, the behavior of A and G are switched, A varies from Amax to −Amax and G remains positive for initial G > 0 (negative for initial G < 0). While this is consistent with the index cycle definition, A remains at its maximum or minimum value for only an instant, therefore the A evolution is inconsistent with the index cycle definition in the supercritical case (for r > 1) as well.

8. Summary

In this paper we provide a complete portrayal of the long-term behavior of the solution of the Lorenz’s (1960) maximum simplification equations. While the solution is always periodic, it can be broadly classified into two groups—the unstable case leading to heteroclinic orbits and the stable case of periodic orbits. It turns out that this dichotomy critically depends on the value of the ratio V/E. The unstable case arise when V/E = 2l2 or 2k2 and the stable case for all other values of this ratio. We also examine the effect of various parameters on the index cycle for both the unstable and stable cases.

In the light of this analysis it should be interesting to refine the data assimilation experiment in Lakshmivarahan et al. (2003) and the properties of the moment dynamics in Epstein (1969) to include the unstable heteroclinic orbits and stable periodic orbits.

Fig. 12.

Heteroclinic orbit starting from η in region ℛ3; η* = (A, 0, 0)T, A = 1, ξ = (0.0, −0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Fig. 12.

Heteroclinic orbit starting from η in region ℛ3; η* = (A, 0, 0)T, A = 1, ξ = (0.0, −0.001, 0.0)T, η = η* + ξ, L = 5000, K = 10 000, r = 0.5, b = 0.1, c = 0.75, (c/b) = 2.7386, bc = 0.27386.

Acknowledgments

We wish to thank John Lewis for stimulating our interest in the lower-order models. When this was in progress, Dan Dobrovolschi of the National Meteorological Administration in Bucharest, Romania, suggested that interpretation of heteroclinic orbits is an interesting problem, which provided further impetus for this work. We wish to record our sincere thanks to Alex Gluhovsky, Yunheng Wang, and the two anonymous reviewers whose comments greatly improved the presentation.

REFERENCES

REFERENCES
Abramowitz
,
M.
, and
I. A.
Stegun
,
1965
:
Handbook of Mathematical Functions.
Dover Publications, 1046 pp
.
Arnol’d
,
V. I.
,
1989
:
Mathematical Methods of Classical Mechanics.
Springer-Verlag, 462 pp
.
De Swart
,
H. E.
,
1989
:
Vacillation and predictability properties of low-order atmospheric spectral models. CWI Tract 60, Center for Mathematics and Computer Science, Amsterdam, Netherlands, 121 pp
.
Devaney
,
R. L.
,
1989
:
An Introduction to Chaotic Dynamical Systems.
2d ed. Addison Wesley, 335 pp
.
Epstein
,
E.
,
1969
:
Stochastic dynamic prediction.
Tellus
,
21
,
739
759
.
Gluhovsky
,
A.
, and
E.
Agee
,
1997
:
An interpretation of atmospheric low-order models.
J. Atmos. Sci.
,
54
,
768
773
.
Gluhovsky
,
A.
, and
C.
Tong
,
1999
:
The structure of energy conserving low-order models.
Phys. Fluids
,
11
,
334
343
.
Gluhovsky
,
A.
,
C.
Tong
, and
E.
Agee
,
2002
:
Selection of modes in convective low-order models.
J. Atmos. Sci.
,
59
,
1383
1393
.
Hirsch
,
M. W.
, and
S.
Smale
,
1974
:
Differential Equations, Dynamical Systems and Linear Algebra.
Academic Press, 358 pp
.
Krishnamurti
,
T. N.
,
H. S.
Bedi
, and
V. M.
Hardiker
,
1998
:
An Introduction to Global Spectral Modeling.
Oxford University Press, 253 pp
.
Lakshmivarahan
,
S.
,
Y.
Honda
, and
J. M.
Lewis
,
2003
:
Second-order approximation to the 3DVAR cost function: Application to analysis and forecast.
Tellus
,
55A
,
371
384
.
Landau
,
L. D.
, and
E. M.
Lifshitz
,
1960
:
Mechanics.
Pergamon Press, 165 pp
.
Lewis
,
J. M.
,
S.
Lakshmivarahan
, and
S. K.
Dhall
,
2006
:
Dynamic Data Assimilation: A Least Squares Approach.
Cambridge University Press, 850 pp
.
Lorenz
,
E. N.
,
1960
:
Maximum simplification of the dynamic equations.
Tellus
,
12
,
243
254
.
Lorenz
,
E. N.
,
1963
:
Deterministic nonperiodic flow.
J. Atmos. Sci.
,
20
,
130
141
.
Namias
,
J.
,
1950
:
The index cycle and its role in the general circulation.
J. Meteor.
,
7
,
130
139
.
Obukhov
,
A. M.
,
1973
:
On the problem of nonlinear interactions in fluid dynamics.
Gerlands Beitr. Geophys.
,
82
,
282
290
.
Saltzman
,
B.
,
1962
:
Finite amplitude free convection as an initial value problem—I.
J. Atmos. Sci.
,
19
,
329
341
.
Sparrow
,
C.
,
1982
:
The Lorenz Equations: Bifurcation, Chaos, and Strange Attractors.
Vol. 41, Applied Mathematical Sciences, Springer Verlag, 259 pp
.
Thompson
,
P. D.
,
1957
:
Uncertainty of initial state as a factor in the predictability of large scale atmosphere pattern.
Tellus
,
9
,
275
295
.
Wiin-Nielson
,
A.
,
1992
:
Comparison of low-order atmospheric dynamical systems.
Atmosfera
,
5
,
135
155
.
Wittenburg
,
J.
,
1977
:
Dynamics of Systems of Rigid Bodies.
B. G. Teubner, 224 pp
.

Footnotes

Corresponding author address: S. Lakshmivarahan, School of Computer Science, University of Oklahoma, Norman, OK 73019. Email: E-mail:varahan@ou.edu