Start Submission Become a Reviewer

Reading: On periodic solutions in the non-dissipative Lorenz model: the role of the nonlinear feedbac...


A- A+
Alt. Display

Original Research Papers

On periodic solutions in the non-dissipative Lorenz model: the role of the nonlinear feedback loop


B.-W. Shen

Department of Mathematics and Statistics, San Diego State University, San Diego, CA, US
About B.-W.

X close


In this study, we discuss the role of the linear heating term and nonlinear terms associated with a non-linear feedback loop in the energy cycle of the three-dimensional (X–Y–Z) non-dissipative Lorenz model (3DNLM), where (X, Y, Z) represent the solutions in the phase space. Using trigonometric functions, we first present the closed-form solution of the nonlinear equation d2X/dτ2 + (X2/2) X = 0 without the heating term (i.e. rX), (where τ is a non-dimensional time and r is the normalized Rayleigh number), a solution that has not been previously documented. Since the solution of the simplified 3D-NLM is oscillatory (wave-like) and since the nonlinear term (X3) is associated with the nonlinear feedback loop, here, we suggest that the nonlinear feedback loop may act as a restoring force. When the heating term is considered, the system yields three critical points. A linear analysis suggests that the origin (i.e. the trivial critical point) is a saddle point and that the other two non-trivial critical points (i.e. centers) are stable. Here, we provide an analysis for three types of solutions that are associated with these critical points. Two of the solutions represent closed curves that travel around one non-trivial critical point or all three critical points. The third type of solution, appearing to connect the stable and unstable manifolds of the saddle point, is called the homoclinic orbit. Using the solution that encloses one non-trivial critical point, here, we show that the competing impact of the nonlinear restoring force and the linear (heating) force determines the partition of the average available potential energy from the Y and Z modes. Based on the energy analysis, an energy cycle with four different regimes is identified. The cycle is only half of a ‘large’ cycle, displaying the wing pattern of a glass-winged butterfly. The other half cycle is anti-symmetric with respect to the origin. The two types of oscillatory solutions with either a small cycle or a large cycle are orbitally stable. As compared to the oscillatory solutions, the homoclinic orbit is not periodic because it “takes forever” to reach the origin. Two trajectories with starting points near the homoclinic orbit may be diverged because one moves with a small cycle and the other moves with a large cycle. Therefore, the homoclinic orbit is not orbitally stable. In a future study, dissipation and/or additional nonlinear terms will be further examined in order to determine how their interactions with the original nonlinear feedback loop and the heating term may change the periodic orbits (as well as homoclinic orbits) to become quasi-periodic orbits and chaotic solutions.

How to Cite: Shen, B.-W., 2018. On periodic solutions in the non-dissipative Lorenz model: the role of the nonlinear feedback loop. Tellus A: Dynamic Meteorology and Oceanography, 70(1), p.1471912. DOI:
  Published on 01 Jan 2018
 Accepted on 16 Apr 2018            Submitted on 22 Nov 2017


It has been 50 years since Lorenz published his breakthrough modeling study (Lorenz, 1963) that changed our views on the predictability of weather and climate (Solomon et al., 2007). The model has been extensively investigated by researchers in various fields including earth science, engineering, mathematics, philosophy, and physics (Gleick, 1987; Sprott, 2003; Jordan and Smith, 2007; Anthes, 2011; Hirsch et al., 2013; Strogatz, 2015). Lorenz’s model with three Fourier modes, which represents the solution to the 2-D Rayleigh–Benard equation (Saltzman, 1962; Lorenz, 1963), is known as the three-dimensional Lorenz model (3DLM). In this paper, we use 3D-NLM to refer to the non-dissipative version of the model that is introduced later in the text.

The scientific community agrees that weather is chaotic, with a finite predictability, and that the source of chaos is nonlinearity. Since the degree of nonlinearity is finite within the 3DLM, the impact of increased nonlinearity on system solutions and/or their stability has been studied using generalized LMs with additional Fourier modes (Curry, 1978; Curry et al., 1984; Howard and Krishnamurti, 1986; Hermiz et al., 1995; Thiffeault and Horton, 1996; Musielak et al., 2005; Roy and Musielak, 2007a). As compared to the 3DLM, some of the generalized LMs have suggested larger Rayleigh number values (or heating parameters) for the onset of chaos, while others have indicated smaller values. This discrepancy may be attributed to different mode truncations (Curry et al., 1984; Thiffeault and Horton, 1996; Roy and Musielak, 2007a,b,c; Shen, 2014, 2015, 2016) that lead to different degrees of nonlinearity and different systems whose energy may or may not be conserved (Roy and Musielak, 2007a).

Among studies using generalized LMs, the pioneering study of Prof. Curry (Curry et al., 1984) suggested that chaotic responses disappeared when sufficient modes are included. Recent studies by Prof. Musielak and co-authors (Musielak et al., 2005; Roy and Musielak, 2007a,b,c) have examined the transition to chaos and the fractal dimensions of generalized LMs, and have emphasized the importance of proper mode truncation in energy conservation. More recent studies (Shen, 2014, 2015, 2016) have discussed the importance of proper Fourier mode selection in extending the nonlinear feedback loop of the 3DLM. The feedback loop, which consists of the two nonlinear terms of the 3DLM, includes a pair of down-scale and up-scale transfer processes associated with the Jacobian function of the partial differential equation, discussed in Section 2. As previously suggested, the original feedback loop may help stabilize the solution for 1 < r < 24.74 within the 3DLM. The extended nonlinear feedback loop in high-dimensional LMs can provide negative nonlinear feedback that produces non-trivial stable critical points when r < 42.9 within a five-dimensional LM and when r < 116.9 within a seven-dimensional LM. Based on the above, it has been hypothesized that a system’s stability can be improved further with additional modes that can provide negative nonlinear feedback. While the importance of an increased degree of nonlinearity with more Fourier modes has been discussed in recent studies, the competing role of the nonlinear terms with the linear (heating) term and/or dissipative terms deserves to be examined in order to ascertain the conditions under which nonlinear processes may lead to stable or chaotic solutions.

Roupas (2012) and others have indicated that the 3DLM in the dissipative limit, which is referred to as the 3D-NLM, contains two conserved quantities that represent the conservation of (KE¯+PE¯) and (KE¯+APE¯), respectively. Here, KE¯,PE¯, and APE¯ are the domain-averaged kinetic energy, potential energy, and available potential energy, respectively. These two quantities, (KE¯+PE¯) and (KE¯+APE¯), are related to the Nambu Hamiltonians (Nambu, 1973; Nevir and Blender, 1994; Floratos, 2012; Roupas, 2012; Blender and Lucarini, 2013). As a result of conservation properties, the collective impact of the nonlinear feedback loop and the linear (heating) term may effectively act as a “restoring” force. The simplicity of the 3D-NLM enables an examination of how the nonlinear feedback loop and the linear (heating) term work together to produce oscillatory solutions (in the phase space). In this work, we address this issue in conjunction with how the available potential energy is partitioned amongst two different Fourier modes, Y and Z, where Z is included in order to enable the nonlinear feedback loop.

The paper is organized as follows: In Section 2, we present the governing equations and the 3D-NLM, introduce the nonlinear feedback loop, and derive the energy conservation laws. In Section 3, we illustrate the role of the nonlinear feedback loop (with r = 0) in acting as a restoring force. With inclusion of the heating term, as well as the nonlinear feedback loop, we present a linear stability analysis for the three critical points and then discuss three types of solutions using analytical and numerical methods. The solutions include two types of oscillatory solutions and the so-called homoclinic orbit (Sprott, 2003; Jordan and Smith, 2007). An energy cycle with four regimes is analyzed based on the tendency of KE¯ and the partition of APE¯ at different scales (i.e. either Y or Z). Concluding remarks are provided at the end. Appendix A discusses the derivation of the 3D-NLM, and Appendix B illustrates the relationship between the nonlinear feedback loop and the nonlinear restoring force. Appendix C presents a closed-form solution to the system with r = 0 using elliptic functions (Davis, 1960). Appendix D discusses the equations and initial conditions that are used to obtain the solution of the homoclinic orbit.


Governing equations and the non-dissipative Lorenz model

The following governing equations for a 2D (x, z), Boussinesq flow are introduced in order to derive the non-dissipative Lorenz model (3D-NLM) and to calculate its kinetic and potential energy:

((1) )
((2) )
where ψ is the streamfunction that yields u=-ψz and w=ψx that represent the horizontal and vertical velocity perturbations, respectively, and θ is the temperature perturbation. ΔT represents the temperature difference between the bottom and top boundaries. The constants, g, α, ν, and κ denote the acceleration of gravity, the coefficient of thermal expansion, the kinematic viscosity, and the thermal diffusivity, respectively. The Jacobian of two arbitrary functions is defined as J(A,B)=(A/x)(B/z)-(A/z)(B/x). The crossout symbol indicates the negligence of a term in the dissipationless limit. Equations (1) and (2), with the dissipative terms, were first used in Saltzman (1962) and Lorenz (1963).

The non-dissipative Lorenz model (3D-NLM) is written as:

((3) )
((4) )
((5) )

Here, (X, Y, Z) represent the amplitude of the Fourier modes. τ=κ(1+a2)(π/H)2t is the dimensionless time. a is the ratio of the vertical scale of the convection cell to its horizontal scale. H is the domain height, and 2H/a represents the domain width. σ=ν/κ is the Prandtl number, and r=Ra/Rc is the normalized Rayleigh number or the heating parameter. Ra is the Rayleigh number, Ra=gαH3ΔT/νκ, and Rc is the critical value for the free-slip Rayleigh–Benard problem, Rc=π4(1+a2)3/a2. The ‘forcing’ terms on the right-hand side of Equations (4) and (5) are referred to as the linear force, or the heating term (rX), and the nonlinear force terms (–XZ and XY). Note that as a result of scale selection in the original Lorenz model, the appearance of σY comes from gαθx, which represents a buoyancy term (Hilborn, 2000), and is not associated with the dissipative terms. Appendix A provides detailed derivations.

The 3D-NLM is integrated forward in time using the fourth-order Runge–Kutta scheme. Without a loss of generality, only three different values of the normalized Rayleigh number, r (r = 0, r = 25, and r = 45) are used, while keeping other parameters, including σ = 10 and a=1/2, constant. A dimensionless time interval (τ) of 0.01 is used and the total number of time steps is 10,000, giving a total dimensionless time (τ) of 100. A smaller τ is used to improve accuracy. In this study, unless stated otherwise, the initial conditions (ICs) are as follows:

((6) )

These settings were used in order to examine the stability of the 5DLM and 6DLM in Shen (2014, 2015), who also discussed the dependence of the solution on different r and σ. To illustrate the unique features of solutions, a very small time step and/or different ICs are used, including (X,Y,Z)=(±ϵ,0,0), (0,±ϵ,0), and (2σr,0,2r). The first IC is used for showing solutions with a small cycle, and the second is used for discussing solutions with a large cycle. The third IC is used to perform a numerical simulation for the homoclinic orbit (Sprott, 2003; Jordan and Smith, 2007). In this study, we use a tiny ϵ in Fig. 5 and ϵ = 1 in Fig. 7.


The nonlinear feedback loop and energy conservation laws

Nonlinear terms in the 3D-NLM (and 3DLM) have been shown to result from the Jacobian term, J(ψ,θ), in Equation (2). The nonlinear interaction of two wave modes via the Jacobian term can generate or impact a third wave mode through a downscale (or upscale) transfer process. The subsequent upscale (or downscale) transfer process associated with the third wave mode can provide feedback to the incipient wave mode(s). As illustrated in Shen (2014), XY and –XZ, respectively, represent the downscale and upscale transfer processes that form a nonlinear feedback loop. When new modes are properly included, the feedback loop can be extended. In the following subsections, we discuss the role of the 3D-NLM nonlinear feedback loop in the energy conservation and partition of available potential energy, which, in turn, helps produce oscillatory solutions.

The domain-averaged kinetic energy (KE¯), the potential energy (PE¯), and the available potential energy (APE¯) are defined as follows (Thiffeault and Horton, 1996; Blender and Lucarini, 2013; Shen, 2014, 2015):

((7) )
((8) )
((9) )

Through straightforward derivations, we obtain the following:

((10) )
((11) )
((12) )
where Co=π2κ2(1+a2a)3. APE¯ is non-positive (i.e. APE¯0), since any perturbation reduces the energy transformable to KE¯.

Equations (10) and (11) lead to:

((13) )
while Equations (10) and (12) yield:
((14) )

With Equations (3)–(5), the time derivative for both Equations (13) and (14) is zero, so these two equations describe energy conservation. Both C1 and C2 are constants and are determined by the initial conditions. Thus, if we express Z and Y2 as functions of X, they are single valued. To facilitate our discussions, the contribution to APE from an individual mode is defined as APEI¯=-CoσI2/2r, where I = Y or Z; therefore, APE¯=APEY¯+APEZ¯. Note that Equations (13) and (14) are related to the two Nambu Hamiltonians, C=-X2/2+σZ and H=Y2/2+Z2/2-rZ (Nambu, 1973; Nevir and Blender, 1994; Floratos, 2012; Roupas, 2012; Blender and Lucarini, 2013).

From the initial conditions in Equation (6), we have C1=0 and C2/Co=-σ/2r, the latter of which is –0.2 for r = 25 and –0.11 for r = 45. Figure 1 provides the time evolution of the conserved quantities: (KE¯+PE¯) and (KE¯+APE¯) in Equations (13) and (14) from the 3D-NLM. At a larger r (e.g. r = 45), a finer τ is required to improve the numerical accuracy of simulated total energy (Fig. 1c). In this study, to facilitate discussions and unless stated otherwise, either C1 or C2 is assumed to be zero.

Fig. 1.  

The time evolution of (KE + PE) and (KE + APE) from the 3D-NLM and 5D-NLM. Panels (a) and (b) are provided for r = 25 and r = 45, respectively. Results obtained from the 5D-NLM are presented in order to show the dependence of solutions on the spatial resolutions. In comparison with panel (b), panel (c) shows the dependence of solutions on the temporal resolution using Δτ = 0.0001. All of the fields are normalized using the constant Co (=π2κ2(1+a2a)3).



In this section, we discuss the competing role of the nonlinear terms and the linear forcing term in the transient solutions and the energy cycle of the 3D-NLM. From Equations (3), (4), and (13), we obtain

((15a) )
((15b) )

The three terms on the right-hand side of Equation (15b) represent the impact of nonlinearity, the linear (heating) force, and the initial conditions, respectively. Their competing impacts (i.e. their differences) determine the sign of M2, and thus, the characteristics of the solution. Based on the relative magnitude of the initial state that may lead to (σr+C1/Co)0 or >0, two types of solutions can be identified (Roupas, 2012). Please note that Equation (15) is a special case of Duffing Equation when (σr+C1/Co)0 and is called the double-well potential system when (σr+C1/Co)0. Some exact solutions to the Duffing equation or double-well potential system can be found in other studies (Drazin, 1992; Jordan and Smith, 2007; Elias-Zuniga, 2013; Starossek, 2015). In this study, we focus on the role of nonlinear feedback loop in producing the oscillatory solutions as well as the homoclinic orbit and on the physical interpretation of the solutions, including the mechanism for diverged trajectories, as discussed below. Therefore, we study the case with 0(σr+C1/Co) by using C1=0 in Section 3.2 and C1/Co=Xo2/2 in Section 3.3, where Xo is the initial value of X. To understand the role of the nonlinear terms (i.e. the nonlinear feedback loop), we begin discussions by solving the solution to the equation that contains no nonlinear terms.


Solutions of the linear system without the nonlinear feedback loop

By assuming no nonlinear terms, we can begin with two equations: dX/dτ=σY and dY/dτ=rX and only one conservation law, as follows:

((16) )

This linear case yields M2=-σr, and Equation (15) becomes d2X/dτ2-σrX=0. A local analysis suggests that the origin, the only critical point in the linear system, is a saddle point. After straightforward derivations, the solution is:

((17) )
where X1 and X2 are constant coefficients. The origin (X,Y)=(0,0) is a saddle point, and the trajectory is hyperbolic with solutions exhibiting exponential growth and decay. The initial condition, which determines dX/dY(=σY/rX), can help select the proper mode(s). For example, (X,Y)=(σ/r,1) only provides the growing mode with (X1,X2)=(0,σ/r), while (X,Y)=(σ/r,-1) leads to the decaying mode with (X1,X2)=(σ/r,0). The former and latter display the properties of unstable and stable manifolds, respectively (Ide et al., 2002). For the nonlinear case, a ‘current’ state (i.e. ICs) may vary with time. Therefore, either mode may appear at different stages. For example, as shown with the homoclinic orbit in Section 3.3, a trajectory with an initial point of (X,Y,Z)=(2σr,0,2r) may approach the origin at a decay rate of σr, while a trajectory beginning near the origin may depart at a growth rate of σr. Based on Equations (16) and (17), although the time change of (KE¯+APE¯) remains zero, the KE¯ produced using only the linear (heating) force has no upper limit. Such an outcome could violate the linear assumption, and, thus, nonlinearity should be included.


Nonlinear solutions and the nonlinear restoring force for r = 0

Here, we examine the impact of nonlinear terms using a special case with r = 0 and (X,Y,Z)=(0,1,0), leading to M2=X2/2. Thus, Equation (15) and its solution are written as follows:

((18a) )
((18b) )
((18c) )

In Appendix B, we show that the cubic term (X3) in Equation (18a) is derived from the nonlinear feedback loop (–XZ and XY). The solution in Equations (18b) and (18c) is periodic and it has time varying frequencies. More importantly, as discussed later, when the approach is applied to the system with r0, we can also obtain periodic and bounded solutions, which is different from the unbounded solution of the linear system in Equation (17). Thus, we suggest that the nonlinear feedback loop (i.e. the cubic term) acts as a restoring force. Solutions for other components (Y and Z) can be obtained as follows.

As compared to the case with r0 in Equations (13) and (14), the energy conservation laws obtained using r = 0 are:

((19) )
((20) )
which, in turn, lead to:
((21) )

Since X2 is defined in Equation (18b), solutions for Y and Z can be obtained as follows:

((22a) )
((22b) )
where the phase function, ϕ in Equation (18c), can also be displayed as:
((23) )

To illustrate the solution’s characteristics, Equations (22a) and (23) are solved using the following iterative method:

((24a) )
((24b) )
where N is the number of iterations. Over a period of time, an initial guess for the phase function is given as ϕ0(τ)=τ. We insert the first phase function, ϕ0, into Equation (24a) to obtain Y0. We then calculate the next phase function, ϕ1, using Y0 and Equation (24b). The integral in Equation (24b) is calculated using the trapezoidal rule. We then repeat the above calculations N times. Numerical results using N = 100 are provided in Fig. 2. The phase function (Fig. 2a) oscillates with time and varies between 0 and π, consistent with Equation (18b) as a result of sin(ϕ)0. Therefore, the solution to Equation (18a) is oscillatory instead of growing or decaying exponentially (as shown in Figs. 2b and 2c). The nonlinear term in Equation (18a) may be viewed as a nonlinear restoring force. Such a suggestion is consistent with the view (Shen, 2014) that the pair of nonlinear terms (–XZ and XY), leading to the nonlinear term (X3/2), can form a nonlinear feedback loop within the 3DLM. Appendix B provides a brief summary. When σ = 10 is replaced by σ = 20, an oscillatory solution with a different period is obtained, as shown in Fig. 2d. As a result of the simple method for integral calculation, a fine τ may be required in order to obtain accurate solutions, as indicated by the red and green lines for the results obtained using τ=0.0001 and 0.01, respectively (Fig. 2a).

Fig. 2.  

Closed-form solutions obtained from Equations (24a,b) using r=0 and the iterative method. (a) The phase function(ϕ) with Δτ=0.01 (green) and Δτ=0.0001 (red) (Equation (24b)). (b) Y=cos(ϕ) (Equation (24a)). (c) X2=2σsin(ϕ) (Equation 18b). (d) The phase function using a different value of σ (=20).

To verify the integral form of the solutions in Equations (24a) and (24b), the numerical solutions of the 3D-NLM using r = 0 (Equations (3)–(5)) are provided in Fig. 3. In panel (a), the blue ‘dot’ indicates the initial temporal evolution of the phase function that is calculated by performing a time integration of X using Equation (18c), where X is obtained from the 3D-NLM; and the red line indicates the phase function calculated using Equations (24a) and (24b). Both are in good agreement. The simulated trajectories in the XY and XZ sections are elliptic and parabolic (Figs. 3b and 3c), respectively, consistent with the analytical relationships in Equations (21) and (19), respectively. Figure 3d provides the time evolution of oscillatory Y (red) and X (blue), consistent with the results provided in Figs. 2b and 2c, respectively. As indicated in Appendix C, the oscillatory characteristics of the solution can also be illustrated using elliptic functions.

Fig. 3.  

Numerical solutions obtained from the 3D-NLM (Equations (3)–(5)) using r=0. (a) A comparison of the phase functions that are calculated from numerical solutions of the 3D-NLM (blue dot) and closed-form solutions of Equation (24a,b) (red). (b) X–Y plot. (c) X–Z plot. (d) The time evolution of X (blue) and Y (red).


Closed form solutions near a non-trivial critical point for r0

In the previous sections, we used Equation (15) to illustrate the individual impact of the linear (heating) force and the nonlinear feedback loop using cases with M2 < 0 and M20, respectively, neither of which changes the sign during the integration. Here, using the initial condition (X,Y,Z)=(Xo,0,0) and 0 < Xo, we consider a more general case with M2=(X2/2-σr-Xo2/2), whose sign may vary during the time integration depending on the relative magnitude of the nonlinearity and the linear (heating) force. M2 has two zeros at X=±Xt, where Xt is defined as 2σr+Xo2. These are called turning points. Based on an analysis using the WKB approximation (Bender and Orszag, 1978), there appears to be a growing or decaying solution for |X|<Xt and an oscillatory (wave-like) solution for |X|>Xt. The former is largely impacted by the linear (heating) force while the latter is impacted by the nonlinear restoring force. Additional analyses are provided below.

To understand the stability of solutions, below, we present a local analysis by linearizing the system with respect to a non-trivial critical point. Using Equations (3)–(5), we first solve for the non-trivial critical point, which is (Xc,Yc,Zc)=(Xc,0,r). Here, Xc can be any constant and is selected as the value of X at the turning point of Equation (15) using the justification presented below. From Equations (3) and (15), the 3D-NLM system can be written as:

((25) )
((26) )

The two equations, with an initial condition of (X, Y) =  (X0,0), lead to three critical points, (Xc, Yc) at (0, 0) and (±2σr+X02, 0). Note that the non-trivial critical points ±Xc are determined by the condition M = 0, which defines the turning points in Equation (15). Therefore, we have Xc = Xt.

Next, each of the total fields is decomposed into its basic part and perturbation, as follows:

((27) )
where V indicates the total fields (X, Y, Z), Vc represents the basic state obtained from the solution of the critical point, and V is a perturbation that measures the departure from the critical point. Using the above equation, the linearized system corresponding to the nonlinear system from Equations (3)–(5) is:
((28) )
where s=(X,Y,Z) and the matrix A that is evaluated at the non-trivial critical point is written as:
((29) )

The above system with Equations (28)–(29) yields eigenvalues of 0 and ±iXc, suggesting that the non-trivial critical point (Xt,0,r) for the linearized system is a stable node as a center. The other non-trivial critical point (−Xt,0,r) shares the same features. One important feature is the dependence of the solution’s period (2π/Xc) on the initial condition (Xo). Next, to illustrate its periodicity, we present a closed form solution for the nonlinear system.

Equations (13)–(14), with the initial conditions (X, Y, Z) = (Xo,0,0), lead to the following two equations:

((30) )
((31) )

Using the same procedures as those provided in Section 3.2 (for r = 0), the following closed-form solutions (for r0) are obtained:

((32a) )
((32b) )
((32c) )
((32d) )

Equations (32c) and (32d) are iteratively computed in order to determine the value of ϕ that is used to calculate the solutions of X, Y, and Z. Note that because Xt2=2σr+Xo2, Equations (32a) and (32c) lead to:

((33) )

Since rrcos(ϕ), X2 is always positive in Equation (32c). Without a loss of generality, X=+2σ(r-rcos(ϕ))+Xo2 is first discussed, yielding a minimum of Xmin=Xo and a maximum of Xmax=4σr+Xo2. In a similar manner, given an initial condition of (X,Y,Z)=(-Xo,0,0), X has a maximum of -Xo and a minimum of -4σr+Xo2. Note that the average of Xmin2 and Xmax2 is equal to 2σr+Xo2, which is equal to Xt2.

Applying the initial conditions (X,Y,Z)=(2σr,0,0) that yield Xc=Xt=2σr (31.62), Fig. 4 shows the closed form solutions obtained using Equation (32), as well as numerical solutions of the 3D-NLM obtained using Equations (3)–(5). The results are in good agreement. For any positive Xo, a trajectory beginning with (X,Y,Z)=(Xo,0,0) is a closed curve. An oscillatory (closed) trajectory associated with a center is orbitally stable (Jordan and Smith, 2007). In Section 3.4, the time varying energy cycle in Equation (32) will be analyzed. Next, we discuss the solutions in (32) for the special case of Xo = 0.

Fig. 4.  

Closed-form solutions obtained from Equations (32c,d) using r=25 and Xo=2σr and the iterative method. (a) A comparison of the closed-form solution (r*sin(ϕ)) from Equation (32a) and the numerical solution of Y from the 3D-NLM. (b) X–Y plot. Closed-form solutions are displayed using black lines and numerical solutions obtained from the 3D-NLM are shown with blue dots. The green line passes through the corresponding critical points at (Xc,Yc)=(2σr+Xo2,0).

When considering a homoclinic solution that begins with and ends at the saddle point, initial conditions of (X,Y,Z)=(0,0,0) are required. At the saddle point, the time derivatives of the 3D-NLM in Equations (3)–(5) are zero. Therefore, numerically, only zero solutions for X, Y, and Z are possible. Obtaining a non-trivial solution can be attempted by initially adding a small perturbation (i.e. (X,Y,Z)=(ϵ,0,0)) and then performing numerical integrations. A similar approach can be applied in Equation (32) to obtain an approximate solution. Here, we provide an initial guess for ϕ,ϕ(τ)=τ, over a target period in Equation (32d). We then calculate X using Equation (32c) and ϕ using Equation (32d) through iterations. Figure 5 displays the solutions determined using Equation (32). The time evolution of Y is provided in Fig. 5a and suggests that Y increases slowly with time, reaches a maximum at τ ≈ 2.3, and then decreases with time to its minimum. The X–Y plot in Fig. 5b appears reasonable and is later compared to the analytical solution. In Fig. 5a, the evolution seems ‘periodic’ but displays some irregularities that are different from those obtained using Xo0. This case is further analyzed using the following analytical solution, as well as the numerical solutions in Section 3.5.

Fig. 5.  

Closed-form solutions with Xo = 0 and the analytical solution of the homoclinic orbit. (a) The closed-form solution (r*sin(ϕ) of Equation (32a). (b) The X–Y plot for the closed-form solutions of Equation (32) (black) and the analytical solutions of the homoclinic orbit using Equation (34) (lightblue). The green line passes through the corresponding critical points at (Xc,Yc)=(2σr+Xo2,0).

From Equation (D3) of Appendix D, the ‘second’ half of the homoclinic orbital solution (e.g. the lower part of the orbit in Fig. 5b) can be obtained as follows:

((34a) )
((34b) )
((34c) )

When t, the above solution approaches the origin. Note that a special point on the solution trajectory is (X,Y,Z)=(2σr,0,2r). Therefore, the ‘second’ half of the solution in Equation (34) describes the trajectory that begins at (X,Y,Z)=(2σr,0,2r) and that moves forward in time toward the origin, as shown in light blue in Fig. 5b. The special point provides an alternative IC for the numerical integration of the 3D-NLM required to simulate (qualitatively) the homoclinic orbit, discussed with numerical results later. When t, X(τ) ≈ 4σre-σrτ and Y(τ) ≈ -4re-σrτ, yielding X/Y=-σ/r. These features are the same as those for the linear (stable) solution with a heating term, as outlined in Section 3.1. As the system is invariant under t-t and y-y (Strogatz, 2015), the solution (X(τ),-Y(τ),Z(τ)) in backward time, which represents the ‘first’ half of the solution, displays the trajectory that begins at the origin and then moves to the point (X,Y,Z)=(2σr,0,2r) (e.g. the upper part of the orbit in Fig. 5b). The total solution that consists of the first and second half connects the stable and unstable manifolds of the saddle point and is called the homoclinic orbit. Since the trajectory ‘takes forever’ to reach the saddle point (as t), it is, therefore, not a periodic solution. As a result, the ‘periodicity’ of closed-form solutions with Xo = 0, as provided in Fig. 5, is spurious, suggesting a numerical error. Unique features of the homoclinic orbit within the analytical solution, as well as the closed-form solution, will be further explored using numerical solutions below.


An energy cycle of solutions with r0

Here, we analyze the time evolution of energy along the homoclinic orbit and apply it to nearby trajectories with a small energy cycle. In Equation (D3), Y20 leads to |X|2σr, which yields a maximum of X (Xm=2σr) at Y = 0, denoted by Ym = 0. Taking the partial derivative of Equation (D3) with respect to X indicates that Y has extrema when X=±2σr=±Xt. At X=±Xt, Z = r, and Y=±r from Equations (13) and (14), respectively, and, thus, APEY¯ = APEZ¯, in other words equal contributions to APE¯ from the Y and Z modes. Furthermore, Equation (D3) suggests that Y2 initially increases in association with the increase in X2 but later decreases in association with an increase in X4. The former is consistent with the linear case in Equation (16), while the latter is consistent with the simplified nonlinear case in Equations (19) and (20). The distribution of Y as a function of X Equation ((D3)), with the aforementioned four points, is provided in Fig. 6. The energy cycle: (1) begins at point A; (2) goes through B, C, and D; and (3) returns back to A. The segment from point P to point Q is designated as Leg PQ, where P and Q are either one of the following: A(X,Y)=(0,0), B=(Xt,Yt),C=(Xm,Ym), or D=(Xt,-Yt). Point A is a saddle point as discussed in Section 3.1, (Xt,Yt)=(2σr,r), and (Xm,Ym)=(2σr,0). The energy cycle is referred to as a small (energy) cycle. From a perspective of potential energy partitioning, the linear force dominates in Legs A–B and D–A where |APEY||APEZ|, while the nonlinear restoring force dominates, in Legs B–C and C–D where |APEY||APEZ|. This also illustrates the importance of the inclusion of the Z mode and thus the nonlinear terms. Here, it should be noted that the homoclinic orbit is not periodic, and that a small energy cycle may appear within the oscillatory trajectories that move close to, but do not pass through, the saddle point. In the discussion provided in the next section, the analysis discussed above is inter-compared to the numerical solutions obtained from the 3D-NLM.

Fig. 6.  

An energy cycle in the 3D non-dissipative Lorenz model (3D-NLM), as shown in the X–Y plot. The four points for the selected ICs are identified, as follows: A(X,Y) = (0,0), B = (2σr,r), C = (2σr,0), and D = (2σr,-r). (σ,r) = (10, 25). The energy cycle: (1) begins at the saddle point A=(0,0); (2) goes through B, C, and D; and (3) returns near A. KE increases as APE decreases along the upper curve (Legs A–B and B–C), and KE decreases as APE increases along the lower curve (Legs C–D and D–A). From a perspective of potential energy partitioning, |APEY||APEZ| in Legs A–B and D–A where the linear force dominates, and |APEY||APEZ| in Legs B–C and C–D where the nonlinear restoring force dominates.


Numerical solutions with r0

In this section, we first discuss the numerical solutions with small cycles that were analyzed in the previous section, and then present the numerical solutions with large cycles that enclose three critical points. The former are obtained using ICs of (X,Y,Z)=(±1,0,0), while the latter are simulated with ICs of (X,Y,Z)=(0,±1,0). Since the features of solutions (e.g. homoclinic orbits) require a high level of numerical accuracy, for Fig. 7, we perform numerical integrations using a very small time step of τ=10-5 but for a shorter period of time (τ = 4). Figure 7a displays oscillatory solutions with small cycles, each of which moves around one of the non-trivial critical points. Figure 7b displays the two oscillatory solutions with a large cycle, moving around all three critical points. The two solutions share the same orbit. The numerical result of the ‘homoclinic orbit’ is shown in Fig. 7c. Figure 7d shows the evolution of these solutions. While the evolution of energy described by the small or large (energy) cycle will be analyzed later, we first discuss the results of ‘homoclinic orbit’.

Fig. 7.  

Three types of solutions associated with different ICs from the 3D-NLM (Equations (3)–(5)) with Δτ = 10−5 and τ = 4. (a) Oscillatory solutions with small cycles. (b) Oscillatory solutions with large cycles. (c) A numerical solution for the homoclinic orbit. (d) The time evolution of the simulated homoclinic orbit (green), and oscillatory solutions with a small cycle (red) and a large cycle (lightblue). Panels (c–d) indicate that the homoclinic trajectory begins at (X, Y, Z) = (2σr 0, 2r), approaches the origin, and remains near the saddle point for a much longer period of time as compared to the oscillatory solutions. Note that as a result of accumulated numerical error, the ‘simulated orbit’ travels around the saddle point, moves to the other side in panel (c), and is no longer the original homoclinic orbit.

The special features of the ‘homoclinic orbit’ are illustrated using numerical simulations with ICs of (X,Y,Z)=(2σr,0,2r). Here, the ‘second’ half of the orbital solution with X > 0 and Y < 0 is a focus. As indicated by the green line in Figs. 7c and 7d, the solution Y decreases and then increases along the homoclinic orbit. Nearly zero values of Y in Fig. 7d suggest that the homoclinic trajectory approaches the saddle point and ‘remains’ close to the saddle point for a long period of time, as indicated by a comparison between the homoclinic orbit solution (in green) and the oscillatory solutions (in light blue and red). Between τ = 0 and τ=1.7 in Fig. 7d, the numerical solution in green represents the ‘second’ half of the homoclinic orbit. In Fig. 7c, the corresponding trajectory in the X–Y phase space that appears in the region with positive values of X compares well with the analytical solution in Fig. 5b. Near the saddle point, (X, Y, Z) are near zero. Therefore, it is easy for numerical errors (e.g. associated with rounding or truncation) to introduce non-zero forcing to drive the flow. When this occurs, the trajectory moves into a region that has negative values of X, and moves around a non-trivial critical point, as shown in Fig. 7c. This move is inconsistent with the analytical solution. However, it is challenging to avoid this type of error since τ is already very small, 10-5. Avoiding this type of error may be more complicated in a higher-dimensional space (e.g. the 5DLM), which is beyond the scope of this study.

The above discussion suggests that the homoclinic orbit separates some trajectories that‘encircle’ a non-trivial critical point from other trajectories that move around the three critical points. Therefore, the homoclinic trajectory that connects the stable and unstable manifolds of the saddle point may be viewed as a separatrix. Two different points that are initially near but not on the homoclinic path will be on periodic orbits that depart from one another, suggesting diverged trajectories. Next, we use the results of Fig. 7a,b as an example to illustrate the divergence of two trajectories starting from (ϵ, 0, 0) and (0, −ϵ, 0). Here, ϵ is a small number. Applying initial conditions (X,Y,Z)=(1,0,0) and (0,-1,0), two trajectories shown with red and light blue lines display ‘repeated’ divergence and convergence. (The two trajectories may give rise to ‘periodic’ divergence and convergence if the ratio of their periods is a rational number.) Therefore, while each of the oscillatory trajectories is orbitally stable, the homoclinic orbit is not. Note that this type of solution dependence on ICs is different from that of a chaotic attractor in dissipative systems (e.g. the 3DLM). Next, we analyze the energy cycle of the oscillatory solutions.

As shown in Fig. 8a using the selected ICs, the large cycle includes two ‘'small cycles’ that are anti-symmetric with respect to the saddle point A. The large cycle resembles the wing of a glasswinged butterfly. The right-hand side of the wing displays the same (or comparable) characteristics as those of the cycle in Fig. 6, while the left-hand side wing is anti-symmetric with respect to the origin (i.e. the saddle point). While the trajectories in XY and XZ are shown in Figs. 8a and 8b, respectively, the distribution of APEY¯ and APEZ¯ as a function of X is provided in Fig. 8c. Here, APEY¯ and APEZ¯ represent averaged available potential energy from the Y and Z modes, respectively. From the perspective of the APE partition, the APEY¯ (red) dominates in Leg AB (DA) while APEZ¯ dominates in Leg BC (CD) where nonlinearity is stronger. Alternatively, KE¯ is predominantly converted from APEY¯ when X2 is small, and is largely converted from APEZ¯ when X2 is large. Therefore, inclusion of the Z mode can lead to the oscillatory solution by enabling the partition of APE on different scales at different stages (i.e. linear and nonlinear stages). More detailed analysis of the energy cycle is provided later. Here, we briefly discuss the impact of different initial conditions.

Fig. 8.  

Solutions obtained from the 3D-NLM (Equations (3)–(5)). (a) X–Y plot. (b) X–Z plot. (c) X-APE plot. The black, red, and blue lines show normalized −KE, APEY, and APEZ, respectively. Green lines are plotted at X=±Xt=±2σr where Y2 = Z2. (d) The X–Y plot using different initial conditions, (X, Y, Z) = (X0, 0, 0). The black, red, and blue curves represent results using X0=2σr,2σr+15, and 2σr+30, respectively. The three green lines pass through the corresponding critical points at (Xc, Yc) = (2σr+X02,0).

In Fig. 8a, a solution may move ‘around’ the saddle point and the two stable critical points, and its trajectory defines the large cycle. With a different initial condition, a solution may go around one of the non-trivial stable points (±2σr+X02, 0), and its trajectory forms a small cycle. For example, Fig. 8d displays results using three different ICs including X0 = 2σr, 2σr + 15, and 2σr + 30 as shown with black, blue, and red lines, respectively. Note that the numerical solution with the first IC was previously compared with the closed-form solution in Fig. 4. The three vertical green lines indicate the locations of the corresponding stable critical points at Xc = 2σr+X02. Solutions in Figs. 8a and 8d suggest that the larger the value of (initial) Xo is, the further the trajectory may move away from the saddle point. When an initial Xo is distant from zero (i.e. the saddle point), nonlinearity may become dominant as compared to the (linear) heating term. The corresponding solutions (e.g. the blue curve in Fig. 8d) become more symmetric with respect to the green line that passes through the nontrivial stable critical point.

Figure 9 provides the time evolution of APEY¯ and APEZ¯ for the solutions in Figs. 8a–8c. The energy cycle: (1) begins near point A, at A+ = (0,0+) to be specific; (2) goes through B, C, and D; and 3) returns back to A, (i.e., A = (0,0)). In Leg AB, where the linear (heating) force dominates, the solution X grows gradually with {KE¯,APEY¯, APEZ¯} and |APEY¯|>|APEZ¯|. The upper arrow (↑) and the down arrow (↓) are, respectively, used to indicate an increase and decrease. In Leg BC (or CD), where the nonlinear restoring force becomes dominant, the solution X increases (or decreases) rapidly with {KE¯, APEY¯,APEZ¯} (or {KE¯, APEY¯,APEZ¯}) and |APEY¯|<|APEZ¯|. In Leg DA, solution X decays slowly with {KE¯,APEY¯, APEZ¯} and |APEY¯|>|APEZ¯|. Once the trajectory returns back to point A, it may experience another small cycle, going to B, C, and D and returning back to point A+. Here, B-=(-Xt,-Yt),C-=(-Xm,Ym), and D-=(-Xt,Yt). The two ‘small cycles’ form the large cycle, and the time evolution of energy is the same for both wings.

Fig. 9.  

Time evolution for the KE and APE in the 3D-NLM (Equations (13)–(14)). Black open circles display KE. Red, blue, and green lines indicate −APEY, −APEZ, and KE +APE, respectively. All of the fields are normalized by Co. The gray line is plotted at a value of σr/2, which is equal to −APEY/Co and −APEZ /Co at X = Xt. These panels show that the solutions are oscillatory and periodic.


Concluding remarks

More than 50 years ago using his forced, dissipative model, Lorenz showed that chaos may appear in the presence of nonlinearity, suggesting that nonlinearity may be the source of chaos. In this study, using the non-dissipative Lorenz model (3D-NLM), we discussed how nonlinearity may act as a restoring force to produce oscillatory solutions. We first presented the closed-form solution of the nonlinear equation d2X/dτ2+(X2/2)X=0, which is derived from the 3D-NLM with r = 0. The corresponding solution is oscillatory (wave-like) and the nonlinear term (X3) is associated with the nonlinear feedback loop. Therefore, we suggest that the nonlinear feedback loop may act as a restoring force. The closed-form solution obtained using a trigonometric function (Equation (18b)) compares well with the closed-form solution obtained using an elliptic function (Equation (C9)); and the simplicity of the former effectively illustrates the fundamental role of nonlinearity in producing oscillatory solutions.

While a heating term is considered, a linear analysis shows a saddle point at the origin and two stable non-trivial points (i.e. two centers) at (Xc,Yc,Zc)=(±Xt,0,r). Here ±Xt depends on the product of σ and r, as well as the initial conditions. In addition, three types of solutions could be obtained. Two of the solutions display periodic movement around one non-trivial critical point or all three critical points. The corresponding trajectories are orbitally stable. The third type of solution is the so-called homoclinic orbit. Two of the trajectories, beginning at points near the homoclinic orbit, may be diverged, displaying the dependence of the solutions on the initial conditions. The homoclinic orbit is not orbitally stable. Note that this solution dependence is different from that associated with a chaotic attractor in a dissipative system that includes two unstable non-trivial critical points.

The homoclinic orbit connects the stable and unstable manifolds of the saddle point, but it is not a periodic solution. However, the time evolution of energy along the homoclinic orbit can qualitatively depict the energy evolution for trajectories that that begin with initial points close to the saddle point and move around a non-trivial critical point. From the perspective of an energy partition, inclusion of the Z mode within the 3D-NLM not only introduces additional APE to be transferred to KE, but it also limits KE to be finite, as compared to the linear system. We illustrated that the relative impact of the nonlinear restoring force and the linear (heating) force determines the partition of the averaged available potential energy associated with the Y and Z modes, denoted as APEY¯ and APEZ¯, respectively. Based on the energy analysis, an energy cycle with four different regimes could be identified with the following four points: A(X,Y)=(0,0),B=(Xt,Yt), C=(Xm,Ym), and D=(Xt,-Yt). With the initial perturbation, (X,Y,Z)=(0,1,0),(Xt,Yt)=(2σr,r), and (Xm,Ym)=(2σr,0). As shown in Fig. 6, the energy cycle may: (1) begin at (near) point A; (2) go through B, C, and D; and (3) return back to A. |APEY¯|<|APEZ¯| appears when solutions are strongly nonlinear (i.e. in Legs B–C and C–D), suggesting that inclusion of the Z mode can enable the partition of APE on different scales, leading to nonlinear solutions. Since point A is a saddle point, the ‘cycle’ is only half of a large cycle, representing one wing of the glass-winged butterfly. More detailed discussions of the energy cycle and the large cycle are provided near the end of Section 3.

In summary, the nonlinear feedback loop alone may act as a restoring force and the heating term alone can produce a saddle point within the 3D-NLM. The nonlinear feedback loop and the heating term collectively lead to three critical points and three types of solutions. The existence of oscillatory solutions reveals the role of the nonlinearity (i.e. the nonlinear feedback loop) in producing recurrence and the appearance of the homoclinic orbit suggests the divergence of two nearby trajectories. This type of solution dependence on ICs will be examined and compared with the one associated with a chaotic attractor in the 3DLM and high-dimensional LMs in a future study.

Fig. C1.  

A comparison of the nonlinear oscillatory solutions using r = 0 in Equations (18b) and (C9), expressed using an elementary trigonometric function and an elliptic function, respectively. (a) The solution of Equation (18b) in red. (b) The solution of Equation (C9) in blue. (c) Both of the solutions in (a) and (b). (d) Differences of the two solutions.


We thank Prof. Zdzislaw Musielak, Dr. Kayo Ide, Dr. Tzi-Cheng Lai, Prof. Xubin Zeng, Prof. Ricardo Carretero, and Prof. Juan Serna for valuable comments and/or discussions, and Jill Dunbar for proofreading the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author.


  1. Anthes, R.2011. Turning the tables on chaos: is the atmosphere more predictable than we assume?, UCAR Magazine, spring/summer. Online at: (last access: 6 January 2016). 

  2. Bender, C. M. and Orszag, S. A.1978. Advanced Mathematical Methods for Scientists and Engineers. McGraw-Hill, New York, pp. 593. 

  3. Blender, R. and Lucarini, V.2013. Nambu representation of an extended Lorenz model with viscous heating. Physica D. 243, 86–91. DOI: 

  4. Curry, J. H.1978. Generalized Lorenz systems. Commun. Math. Phys. 60, 193–204. 

  5. Curry, J. H., Herring, J. R., Loncaric, J. and Orszag, S. A.1984. Order and disorder in two- and three-dimensional Benard convection. J. Fluid Mech. 147, 1–38. DOI: 

  6. David, D. and HolmD. D.1992. Multiple Lie-Poisson structures, reductions, and geometric phases for the Maxwell-Bloch travelling wave equations. J. Nonlinear Sci. 2, 241–262. DOI: 

  7. Davis, H. T.1960. Introduction to Bonlinear Differential and Integral Equations. Dover, New York, pp. 566.  

  8. Drazin, P. G.1992. Nonlinear Systems. Cambridge University Press, New York, pp. 333. 

  9. Elias-Zuniga, A.2013. Exact solution of the cubic-quintic Duffing oscillator. Appl. Math. Model. 37, 2574–2579. DOI: 

  10. Floratos, E.2012. Matrix quantization of turbulence. Int. J. Bifurcat. Chaos. 22, 1250213. DOI: 

  11. Gleick, J.1987. Chaos: Making a New Science. Penguin, New York, pp. 360. 

  12. Hermiz, K. B., Guzdar, P. N. and Finn, J. M.1995. Improved low-order model for shear flow driven by Rayleigh–Benard convection. Phys. Rev. E. 51, 325–331. DOI: 

  13. Hilborn, R. C.2000. Chaos and Nonlinear Dynamics. An Introduction for Scientists and Engineers. 2nd ed. Oxford University Press, New York, pp. 650. 

  14. Hirsch, M., Smale, S. and Devaney, R. L.2013. Differential Equations, Dynamical Systems, and an Introduction to Chaos. 3rd ed. Academic Press, San Diego, pp. 432. 

  15. Howard, L. N. and Krishnamurti, R. K.1986. Large-scale flow in turbulent convection: a mathematical model. J. Fluid Mech. 170, 385–410. DOI: 

  16. Ide, K., Small, D. and Wiggins, S.2002. Distinguished hyperbolic trajectories in time-dependent fluid flows: analytical and computational approach for velocity fields defined as data sets. Nonlin. Processes Geophys. 9, 237–263. DOI: 

  17. Jordan, D. W. and Smith, P.2007. Nonlinear Ordinary Differential Equations. An Introduction for Scientists and Engineers. 4th ed. Oxford University Press, New York, pp. 560. 

  18. Lorenz, E.1963. Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130–141. DOI:<0130:DNF>2.0.CO;2. 

  19. Musielak, Z. E., Musielak, D. E. and Kennamer, K. S.2005. The onset of chaos in nonlinear dynamical systems determined with a new fractal technique. Fractals. 13, 19–31. DOI: 

  20. Nambu, Y.1973. Generalized Hamiltonian dynamics. Phys. Rev. D. 7, 2403. DOI: 

  21. Nevir, P. and Blender, R.1994. Hamiltonian and Nambu representation of the non-dissipative Lorenz equations. Contrib. to Atmosph. Phys. 67, 133–140. DOI: 

  22. Roupas, Z.2012. Phase space geometry and chaotic attractors in dissipative Nambu mechanics. J. Phys. A-Math. Theor. 45, 195101. DOI: 

  23. Roy, D. and Musielak, Z. E.2007a. Generalized Lorenz models and their routes to chaos. I. energy-conserving vertical mode truncations. Chaos Soliton. Fract. 32, 1038–1052. DOI: 

  24. Roy, D. and Musielak, Z. E.2007b. Generalized Lorenz models and their routes to chaos. II. Energy-conserving horizontal mode truncations. Chaos Soliton. Fract. 31, 747–756. DOI: 

  25. Roy, D. and Musielak, Z. E.2007c. Generalized Lorenz models and their routes to chaos. III. Energy-conserving horizontal and vertical mode truncations. Chaos Soliton. Fract. 33, 1064–1070. DOI: 

  26. Saltzman, B.1962. Finite amplitude free convection as an initial value problem. J. Atmos. Sci. 19, 329–341. DOI:<0329:FAFCAA>2.0.CO;2. 

  27. Shen, B.-W.2014. Nonlinear feedback in a five-dimensional Lorenz model. J. Atmos. Sci. 71, 1701–1723. DOI: 

  28. Shen, B.-W.2015. Nonlinear feedback in a six-dimensional Lorenz Model. Impact of an additional heating term. Nonlin. Processes Geophys. 22, 749–764. DOI: 

  29. Shen, B.-W.2016. Hierarchical scale dependence associated with the extension of the nonlinear feedback loop in a seven-dimensional Lorenz model. Nonlin. Processes Geophys. 23, 189–203. DOI: 

  30. Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M. and co-authors. 2007. Climate Change 2007: The Physical Science Basis. Cambridge University Press, New York, pp. 996.  

  31. Sparrow, C.1982. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. Springer, New York, pp. 41. 

  32. Sprott, J. C.2003. Chaos and Time-Series Analysis. Oxford University Press. New York, pp. 528. The numerical method is briefly discussed on 

  33. Starossek, U.2015. A low-frequency pendulum mechanism. Mech. Mach. Theory. 83, 81–90. DOI: 

  34. Strogatz, S. H.2015. Nonlinear Dynamics and Chaos. With Applications to Physics, Biology, Chemistry, and Engineering. Westpress view, Boulder, CO, pp. 513.  

  35. Thiffeault, J.-L. and Horton, W.1996. Energy-conserving truncations for convection with shear flow. Phys. Fluids. 8, 1715–1719. DOI: 

Appendix A

Derivations of Equation (3)

In this section, we describe the procedures used to derive Equation (3). As discussed in Shen (2014), the nonlinear advection term in Equation (1) does not explicitly appear within the Lorenz model. Therefore, Equation (1) can be written as follows:

((A1) )

Using the M1-M3 of Shen (2014), ψ and θ can be represented as:

((A2) )
((A3) )

Here, C1 and C2 are defined in Equation (9) of Shen (2014). Note that the two constants C1 and C2 are different from those in Equations (13) and (14). Next, the left- and right-hand side terms of Equation (A1) are expressed using M1-M3, respectively, as follows:

((A4) )
((A5) )

Here, l and m, as defined in Shen (2014), represent the horizontal and vertical wavenumbers, respectively. Substituting Equations (A4) and (A5) into Equation (A1), we have:

((A6) )

We define a time scale (To) as:

((A7) )
and obtain:
((A8) )

The term H represents the domain height, and a is the ratio of the vertical scale of the convection cell to its horizontal scale, (i.e. a=l/m). With Equation (9) in Shen (2014), the right-hand side of Equation (A6) becomes:

((A9) )

With Equation (A9), Equation (A6) becomes:

((A10) )

The above equation then becomes Equation (3), dX/dτ=σY, when τ=κTot is assumed.

Appendix B

Nonlinear feedback loop and nonlinear restoring force

The non-dissipative Lorenz model (3D-NLM) with no heating term is written as:

((B1) )
((B2) )
((B3) )

Note that the above equations represent a special system with a very large r for the original 3DLM (e.g. Sparrow, 1982; Strogatz, 2015). Equations (B1)–(B2) lead to

((B4) )

From Equation (B1) and (B3) (i.e. X * Equation (B1)− σ*Equation (B3)), we have

((B5) )
((B6) )
where C is a constant. The initial condition of (X,Y,Z)=(0,1,0) leads to C = 0. Therefore, Equations (B4) and (B5) with the above initial condition gives us
((B7) )

The above derivations illustrate the relationship between the X3 and the nonlinear feedback loop (i.e. –XZ and XY), which are associated with a pair of up-scaling and down-scaling processes as discussed in Shen (2014). Note that replacing (-XZ) and (XY) by (XZ) and (-XY) in Equations (B2) and (B3), respectively, we obtain the real-valued Maxwell-Bloch equations (Equations (1.20)) of David and Holm, 1992). Therefore, it is suggested that the pair of (XZ) and (-XY) serves as a restoring force in the Maxwell–Bloch equations.

Appendix C

A closed-form solution for r = 0 using elliptic functions

Nonlinear dynamics have been widely studied using solutions to the nonlinear Duffing equation. The Duffing equation that governs certain damped and driven (i.e. dissipative and forced) oscillatory motions (Bender and Orszag, 1978; Jordan and Smith, 2007; Wikipedia) is given by:

((C1) )
where δ, α, β, γ, and ω are constants. When δ=γ=0 and α=-σr and β=1/2, Equation (C1) appears in the form of the 3D-NLM (Equation (15)). While an exact solution to the general form of the Duffing equation may not have been determined, closed-form, approximated, or numerical solutions have been documented in the literature. In this section, we discuss how closed-form solutions to the special case of the non-dissipative and unforced Duffing equation (i.e. δ=α=γ=0 and β=1/2) can be obtained and expressed in terms of elliptic functions (Davis, 1960).

The incomplete elliptic integral of the first kind, F(ϕ,k), is written as:

((C2) )
u=F(ϕ,k)=0ϕdθ1-k2sin2θ.  0<k2<1.

The elliptic functions of the Jacobi are defined as the inverse of the elliptic integral, as follows:

((C3) )
sn(u,k)=sinϕ,  cn(u,k)=cosϕ.

The elliptic functions have the following properties:

((C4) )
sn(0)=0,  cn(0)=1,  sn2u+cn2=1,
and their derivatives can be obtained as follows:
((C5) )
((C6) )

Equations (C4) and (C5) can be helpful for obtaining a solution to the equation d2X/dτ2+X3/2=0 (Equation (18a)). Using Equation (C5) and k2=-1, the solution can be written as:

((C7) )

Through Equation (C6) and k2=1/2, the solution becomes:

((C8) )

Here, D1, D2, θ1, or θ2 are constants and are determined from the initial conditions. θ1 or θ2 can be expressed in terms of complete elliptic integrals of the first kind, K(k)=F(π/2,k2=1/2)1.8541, which is calculated using Matlab ( In addition, as a result of the definition of the elliptic integral (Equation (C2)) and/or validity of the numerical algorithms for the elliptic functions, k2 may not be equal to –1. Therefore, only Equation (C8) is discussed below. In this study, the IC with X(0)=0 and Y(0) = 1 (i.e., ddτX(0)=σ)) leads to the following solution with θ2=3K and D2=2σ:

((C9) )

Figure C1 displays solutions obtained from Equation (18b) (Fig. C1-a) and Equation (C9) (Fig. C1-b) and their differences (Fig. C1-d). The solutions are very close. This comparison confirms the validity of the method for obtaining closed-form solutions for the 3D-NLM with r = 0 using elementary trigonometric functions (Fig. 2) that are used to illustrate the role of nonlinearity in producing oscillatory solutions. For the case of α≠0 (i.e. r≠0), similar procedures can be applied using Equation (C5) or (C6) in order to obtain closed-form solutions. The procedures are beyond the scope of the present study.

Appendix D

Equations for the solution of the homoclinic orbit

Here, we present equations that are used for obtaining an analytical solution that begins and ends at the saddle point. Equations (13)–(14) with the ICs (X,Y,Z)=(0,0,0) become:

((D1) )
((D2) )
which lead to:
((D3) )

Using the above with Equation (3) yields:

((D4) )

The term X4 in Equation (D4) can be shown to be equivalent to the term X3 in Equation (15) that is associated with the nonlinear feedback loop. Without this term, the system in Equation (D4) or (15) becomes linear. Applying the method described in Chapter 3 of Jordan and Smith (2007), we obtain an analytical solution to the nonlinear system in Equation (D3) that begins and ends at the saddle point (0,0,0), which are shown in Equations (34a)–(34c).

comments powered by Disqus