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 ($\overline{\text{KE}}+\overline{\text{PE}}$) and ($\overline{\text{KE}}+\overline{\text{APE}}$), respectively. Here, $\overline{\text{KE}},\overline{\text{PE}}$, and $\overline{\text{APE}}$ are the domain-averaged kinetic energy, potential energy, and available potential energy, respectively. These two quantities, ($\overline{\text{KE}}+\overline{\text{PE}}$) and ($\overline{\text{KE}}+\overline{\text{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 $\overline{\text{KE}}$ and the partition of $\overline{\text{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.
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:
The non-dissipative Lorenz model (3D-NLM) is written as:
Here, (X, Y, Z) represent the amplitude of the Fourier modes. $\tau =\kappa (1+{a}^{2}){(\pi /H)}^{2}t$ 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. $\sigma =\nu /\kappa $ is the Prandtl number, and $r={R}_{\text{a}}/{R}_{\text{c}}$ is the normalized Rayleigh number or the heating parameter. ${R}_{\text{a}}$ is the Rayleigh number, ${R}_{\text{a}}=g\alpha {H}^{3}\Delta T/\nu \kappa $, and ${R}_{\text{c}}$ is the critical value for the free-slip Rayleigh–Benard problem, ${R}_{\text{c}}={\pi}^{4}{(1+{a}^{2})}^{3}/{a}^{2}$. 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\alpha \frac{\partial \theta}{\partial 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/\sqrt{2}$, constant. A dimensionless time interval ($\u25b3\tau $) of 0.01 is used and the total number of time steps is 10,000, giving a total dimensionless time (τ) of 100. A smaller $\u25b3\tau $ is used to improve accuracy. In this study, unless stated otherwise, the initial conditions (ICs) are as follows:
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)=(\pm \mathrm{\u03f5},0,0)$, $(0,\pm \mathrm{\u03f5},0)$, and $(2\sqrt{\sigma 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.
Nonlinear terms in the 3D-NLM (and 3DLM) have been shown to result from the Jacobian term, $J(\psi ,\theta )$, 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 ($\overline{\text{KE}}$), the potential energy ($\overline{\text{PE}}$), and the available potential energy ($\overline{\text{APE}}$) are defined as follows (Thiffeault and Horton, 1996; Blender and Lucarini, 2013; Shen, 2014, 2015):
Through straightforward derivations, we obtain the following:
Equations (10) and (11) lead to:
With Equations (3)–(5), the time derivative for both Equations (13) and (14) is zero, so these two equations describe energy conservation. Both C_{1} and C_{2} are constants and are determined by the initial conditions. Thus, if we express Z and Y^{2} as functions of X, they are single valued. To facilitate our discussions, the contribution to APE from an individual mode is defined as $\overline{{\text{APE}}_{I}}=-{C}_{o}\sigma {I}^{2}/2r$, where I = Y or Z; therefore, $\overline{\text{APE}}=\overline{{\text{APE}}_{Y}}+\overline{{\text{APE}}_{Z}}$. Note that Equations (13) and (14) are related to the two Nambu Hamiltonians, $C=-{X}^{2}/2+\sigma Z$ and $H={Y}^{2}/2+{Z}^{2}/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 ${C}_{1}=0$ and ${C}_{2}/{C}_{o}=-\sigma /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: ($\overline{\text{KE}}+\overline{\text{PE}}$) and ($\overline{\text{KE}}+\overline{\text{APE}}$) in Equations (13) and (14) from the 3D-NLM. At a larger r (e.g. r = 45), a finer $\u25b3\tau $ is required to improve the numerical accuracy of simulated total energy (Fig. 1c). In this study, to facilitate discussions and unless stated otherwise, either C_{1} or C_{2} is assumed to be zero.
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
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 M^{2}, and thus, the characteristics of the solution. Based on the relative magnitude of the initial state that may lead to $(\sigma r+{C}_{1}/{C}_{o})\le 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 $(\sigma r+{C}_{1}/{C}_{o})\le 0$ and is called the double-well potential system when $(\sigma r+{C}_{1}/{C}_{o})\ge 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\le (\sigma r+{C}_{1}/{C}_{o})$ by using ${C}_{1}=0$ in Section 3.2 and ${C}_{1}/{C}_{o}={X}_{o}^{2}/2$ in Section 3.3, where X_{o} 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.
By assuming no nonlinear terms, we can begin with two equations: $\text{d}X/\text{d}\tau =\sigma Y$ and $\text{d}Y/\text{d}\tau =rX$ and only one conservation law, as follows:
This linear case yields ${M}^{2}=-\sigma r$, and Equation (15) becomes ${\text{d}}^{2}X/\text{d}{\tau}^{2}-\sigma 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:
Here, we examine the impact of nonlinear terms using a special case with r = 0 and $(X,Y,Z)=(0,1,0)$, leading to ${M}^{2}={X}^{2}/2$. Thus, Equation (15) and its solution are written as follows:
In Appendix B, we show that the cubic term (X^{3}) 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 $r\ne 0$, 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 $r\ne 0$ in Equations (13) and (14), the energy conservation laws obtained using r = 0 are:
Since X^{2} is defined in Equation (18b), solutions for Y and Z can be obtained as follows:
To illustrate the solution’s characteristics, Equations (22a) and (23) are solved using the following iterative method:
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 X–Y and X–Z 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.
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 M^{2} < 0 and ${M}^{2}\ge 0$, respectively, neither of which changes the sign during the integration. Here, using the initial condition $(X,Y,Z)=({X}_{o},0,0)$ and $0\hspace{0.17em}<\hspace{0.17em}{X}_{o}$, we consider a more general case with ${M}^{2}=({X}^{2}/2-\sigma r-{X}_{o}^{2}/2)$, whose sign may vary during the time integration depending on the relative magnitude of the nonlinearity and the linear (heating) force. M^{2} has two zeros at $X=\pm {X}_{\text{t}}$, where ${X}_{\text{t}}$ is defined as $\sqrt{2\sigma r+{X}_{o}^{2}}$. 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 $\left|X\right|\hspace{0.17em}<\hspace{0.17em}{X}_{\text{t}}$ and an oscillatory (wave-like) solution for $\left|X\right|>{X}_{\text{t}}$. 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 $({X}_{c},{Y}_{c},{Z}_{c})=({X}_{c},0,r)$. Here, X_{c} 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:
The two equations, with an initial condition of (X, Y) = $({X}_{0},0)$, lead to three critical points, (X_{c}, Y_{c}) at (0, 0) and $(\pm \sqrt{2\sigma r+{X}_{0}^{2}}\mathrm{,}0)$. Note that the non-trivial critical points $\pm {X}_{c}$ are determined by the condition M = 0, which defines the turning points in Equation (15). Therefore, we have X_{c} = X_{t}.
Next, each of the total fields is decomposed into its basic part and perturbation, as follows:
The above system with Equations (28)–(29) yields eigenvalues of 0 and $\pm i{X}_{c}$, suggesting that the non-trivial critical point (${X}_{t},0,r$) for the linearized system is a stable node as a center. The other non-trivial critical point (−${X}_{t},0,r$) shares the same features. One important feature is the dependence of the solution’s period ($2\pi /{X}_{c}$) on the initial condition (X_{o}). 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) = (${X}_{o},0,0$), lead to the following two equations:
Using the same procedures as those provided in Section 3.2 (for r = 0), the following closed-form solutions (for $r\ne 0)$ are obtained:
Equations (32c) and (32d) are iteratively computed in order to determine the value of $\varphi $ that is used to calculate the solutions of X, Y, and Z. Note that because ${X}_{t}^{2}=2\sigma r+{X}_{o}^{2}$, Equations (32a) and (32c) lead to:
Since $r\ge rcos(\varphi )$, X^{2} is always positive in Equation (32c). Without a loss of generality, $X=+\sqrt{2\sigma (r-rcos(\varphi ))+{X}_{o}^{2}}$ is first discussed, yielding a minimum of ${X}_{min}={X}_{o}$ and a maximum of ${X}_{max}=\sqrt{4\sigma r+{X}_{o}^{2}}$. In a similar manner, given an initial condition of $(X,Y,Z)=(-{X}_{o},0,0)$, X has a maximum of $-{X}_{o}$ and a minimum of $-\sqrt{4\sigma r+{X}_{o}^{2}}$. Note that the average of ${X}_{min}^{2}$ and ${X}_{max}^{2}$ is equal to $2\sigma r+{X}_{o}^{2}$, which is equal to ${X}_{t}^{2}$.
Applying the initial conditions $(X,Y,Z)=(\sqrt{2\sigma r},0,0)$ that yield ${X}_{c}={X}_{t}=2\sqrt{\sigma r}$ ($\approx 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 X_{o}, a trajectory beginning with $(X,Y,Z)=({X}_{o},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 X_{o} = 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)=(\mathrm{\u03f5},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 $\varphi ,\varphi (\tau )=\tau $, over a target period in Equation (32d). We then calculate X using Equation (32c) and $\varphi $ 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 $\tau \hspace{0.17em}\approx \hspace{0.17em}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 ${X}_{o}\ne 0$. This case is further analyzed using the following analytical solution, as well as the numerical solutions in Section 3.5.
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:
When $t\to \infty $, the above solution approaches the origin. Note that a special point on the solution trajectory is $(X,Y,Z)=(2\sqrt{\sigma r},0,2r)$. Therefore, the ‘second’ half of the solution in Equation (34) describes the trajectory that begins at $(X,Y,Z)=(2\sqrt{\sigma 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\to \infty $, $X(\tau )\hspace{0.17em}\approx \hspace{0.17em}4\sqrt{\sigma r}{e}^{-\sqrt{\sigma r}\tau}$ and $Y(\tau )\hspace{0.17em}\approx \hspace{0.17em}-4r{e}^{-\sqrt{\sigma r}\tau}$, yielding $X/Y=-\sqrt{\sigma /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\to -t$ and $y\to -y$ (Strogatz, 2015), the solution ($X(\tau ),-Y(\tau ),Z(\tau ))$ 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\sqrt{\sigma 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\to \infty $), it is, therefore, not a periodic solution. As a result, the ‘periodicity’ of closed-form solutions with X_{o} = 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.
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), ${Y}^{2}\ge 0$ leads to $\left|X\right|\le 2\sqrt{\sigma r}$, which yields a maximum of X (${X}_{\text{m}}=2\sqrt{\sigma r}$) at Y = 0, denoted by Y_{m} = 0. Taking the partial derivative of Equation (D3) with respect to X indicates that Y has extrema when $X=\pm \sqrt{2\sigma r}=\pm {X}_{\text{t}}$. At $X=\pm {X}_{\text{t}}$, Z = r, and $Y=\pm r$ from Equations (13) and (14), respectively, and, thus, $\overline{{\text{APE}}_{Y}}$ = $\overline{{\text{APE}}_{Z}}$, in other words equal contributions to $\overline{\text{APE}}$ from the Y and Z modes. Furthermore, Equation (D3) suggests that Y^{2} initially increases in association with the increase in X^{2} but later decreases in association with an increase in X^{4}. 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 P–Q, where P and Q are either one of the following: $A(X,Y)=(0,0)$, $B=({X}_{\text{t}},{Y}_{\text{t}}),C=({X}_{\text{m}},{Y}_{\text{m}})$, or $D=({X}_{\text{t}},-{Y}_{\text{t}})$. Point A is a saddle point as discussed in Section 3.1, $({X}_{\text{t}},{Y}_{\text{t}})=(\sqrt{2\sigma r},r)$, and $({X}_{\text{m}},{Y}_{\text{m}})=(2\sqrt{\sigma 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 $\left|\mathrm{A}\mathrm{P}{\mathrm{E}}_{Y}\right|\ge \left|\mathrm{A}\mathrm{P}{\mathrm{E}}_{Z}\right|$, while the nonlinear restoring force dominates, in Legs B–C and C–D where $\left|\mathrm{A}\mathrm{P}{\mathrm{E}}_{Y}\right|\le \left|\mathrm{A}\mathrm{P}{\mathrm{E}}_{Z}\right|$. 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.
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)=(\pm 1,0,0)$, while the latter are simulated with ICs of $(X,Y,Z)=(0,\pm 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 $\u25b3\tau ={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’.
The special features of the ‘homoclinic orbit’ are illustrated using numerical simulations with ICs of $(X,Y,Z)=(2\sqrt{\sigma 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 $\tau =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 $\u25b3\tau $ 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 X–Y and X–Z are shown in Figs. 8a and 8b, respectively, the distribution of $\overline{{\text{APE}}_{Y}}$ and $\overline{{\text{APE}}_{Z}}$ as a function of X is provided in Fig. 8c. Here, $\overline{{\text{APE}}_{Y}}$ and $\overline{{\text{APE}}_{Z}}$ represent averaged available potential energy from the Y and Z modes, respectively. From the perspective of the APE partition, the $\overline{{\text{APE}}_{Y}}$ (red) dominates in Leg A–B (D–A) while $\overline{{\text{APE}}_{Z}}$ dominates in Leg B–C (C–D) where nonlinearity is stronger. Alternatively, $\overline{\text{KE}}$ is predominantly converted from $\overline{{\text{APE}}_{Y}}$ when X^{2} is small, and is largely converted from $\overline{{\text{APE}}_{Z}}$ when X^{2} 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.
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 ($\pm \sqrt{2\sigma r+{X}_{0}^{2}}$, 0), and its trajectory forms a small cycle. For example, Fig. 8d displays results using three different ICs including X_{0} = $\sqrt{2\sigma r}$, $\sqrt{2\sigma r}$ + 15, and $\sqrt{2\sigma 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 X_{c} = $\sqrt{2\sigma r+{X}_{0}^{2}}$. Solutions in Figs. 8a and 8d suggest that the larger the value of (initial) X_{o} is, the further the trajectory may move away from the saddle point. When an initial X_{o} 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 $\overline{{\text{APE}}_{Y}}$ and $\overline{{\text{APE}}_{Z}}$ 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 A–B, where the linear (heating) force dominates, the solution X grows gradually with {$\overline{\text{KE}}\uparrow ,\overline{{\text{APE}}_{Y}}\downarrow $, $\overline{{\text{APE}}_{Z}}\downarrow $} and $\left|\overline{{\text{APE}}_{Y}}\right|>\left|\overline{{\text{APE}}_{Z}}\right|$. The upper arrow (↑) and the down arrow (↓) are, respectively, used to indicate an increase and decrease. In Leg B–C (or C–D), where the nonlinear restoring force becomes dominant, the solution X increases (or decreases) rapidly with {$\overline{\text{KE}}\uparrow $, $\overline{{\text{APE}}_{Y}}\uparrow ,\hspace{0.17em}\overline{{\text{APE}}_{Z}}\downarrow $} (or {$\overline{\text{KE}}\downarrow $, $\overline{{\text{APE}}_{Y}}\downarrow ,\hspace{0.17em}\overline{{\text{APE}}_{Z}}\uparrow $}) and $\left|\overline{{\text{APE}}_{Y}}\right|<\left|\overline{{\text{APE}}_{Z}}\right|$. In Leg D–A, solution X decays slowly with {$\overline{\text{KE}}\downarrow ,\hspace{0.17em}\overline{{\text{APE}}_{Y}}\uparrow $, $\overline{{\text{APE}}_{Z}}\uparrow $} and $\left|\overline{{\text{APE}}_{Y}}\right|>\left|\overline{{\text{APE}}_{Z}}\right|$. 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}^{-}=(-{X}_{\text{t}},-{Y}_{\text{t}}),{C}^{-}=(-{X}_{\text{m}},{Y}_{\text{m}})$, and ${D}^{-}=(-{X}_{\text{t}},{Y}_{\text{t}})$. The two ‘small cycles’ form the large cycle, and the time evolution of energy is the same for both wings.
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 ${\mathrm{d}}^{2}X/\mathrm{d}{\tau}^{2}+({X}^{2}/2)X=0$, which is derived from the 3D-NLM with r = 0. The corresponding solution is oscillatory (wave-like) and the nonlinear term (X^{3}) 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 $({X}_{c},{Y}_{c},{Z}_{c})=(\pm {X}_{t},0,r)$. Here $\pm {X}_{t}$ 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 $\overline{{\text{APE}}_{Y}}$ and $\overline{{\text{APE}}_{Z}}$, 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=({X}_{\text{t}},{Y}_{\text{t}})$, $C=({X}_{\text{m}},{Y}_{\text{m}})$, and $D=({X}_{\text{t}},-{Y}_{\text{t}})$. With the initial perturbation, $(X,Y,Z)=(0,1,0),({X}_{\text{t}},{Y}_{\text{t}})=(\sqrt{2\sigma r},r)$, and $({X}_{\text{m}},{Y}_{\text{m}})=(2\sqrt{\sigma 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. $\left|\overline{{\text{APE}}_{Y}}\right|<\left|\overline{{\text{APE}}_{Z}}\right|$ 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.
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.
No potential conflict of interest was reported by the author.
Anthes, R.2011. Turning the tables on chaos: is the atmosphere more predictable than we assume?, UCAR Magazine, spring/summer. Online at: https://www2.ucar.edu/atmosnews/opinion/turning-tables-chaos-atmosphere-more-predictable-we-assume-0 (last access: 6 January 2016).
Bender, C. M. and Orszag, S. A.1978. Advanced Mathematical Methods for Scientists and Engineers. McGraw-Hill, New York, pp. 593.
Blender, R. and Lucarini, V.2013. Nambu representation of an extended Lorenz model with viscous heating. Physica D. 243, 86–91. DOI:https://doi.org/10.1016/j.physd.2012.09.007.
Curry, J. H.1978. Generalized Lorenz systems. Commun. Math. Phys. 60, 193–204. https://projecteuclid.org/euclid.cmp/1103904126
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:https://doi.org/10.1017/S0022112084001968.
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:https://doi.org/10.1007/BF02429857.
Davis, H. T.1960. Introduction to Bonlinear Differential and Integral Equations. Dover, New York, pp. 566.
Drazin, P. G.1992. Nonlinear Systems. Cambridge University Press, New York, pp. 333.
Elias-Zuniga, A.2013. Exact solution of the cubic-quintic Duffing oscillator. Appl. Math. Model. 37, 2574–2579. DOI:https://doi.org/10.1016/j.apm.2012.04.005.
Floratos, E.2012. Matrix quantization of turbulence. Int. J. Bifurcat. Chaos. 22, 1250213. DOI:https://doi.org/10.1142/S0218127412502136.
Gleick, J.1987. Chaos: Making a New Science. Penguin, New York, pp. 360.
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:https://doi.org/10.1103/PhysRevE.51.325.
Hilborn, R. C.2000. Chaos and Nonlinear Dynamics. An Introduction for Scientists and Engineers. 2nd ed. Oxford University Press, New York, pp. 650.
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.
Howard, L. N. and Krishnamurti, R. K.1986. Large-scale flow in turbulent convection: a mathematical model. J. Fluid Mech. 170, 385–410. DOI:https://doi.org/10.1017/S0022112086000940.
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:https://doi.org/10.5194/npg-9-237-2002.
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.
Lorenz, E.1963. Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130–141. DOI:https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.
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:https://doi.org/10.1142/S0218348X0500274X.
Nambu, Y.1973. Generalized Hamiltonian dynamics. Phys. Rev. D. 7, 2403. DOI:https://doi.org/10.1103/PhysRevD.7.2405.
Nevir, P. and Blender, R.1994. Hamiltonian and Nambu representation of the non-dissipative Lorenz equations. Contrib. to Atmosph. Phys. 67, 133–140. DOI:https://doi.org/10.1088/0305-4470/26/22/010.
Roupas, Z.2012. Phase space geometry and chaotic attractors in dissipative Nambu mechanics. J. Phys. A-Math. Theor. 45, 195101. DOI:https://doi.org/10.1088/1751-8113/45/19/195101.
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:https://doi.org/10.1016/j.chaos.2006.02.013.
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:https://doi.org/10.1016/j.chaos.2006.03.082.
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:https://doi.org/10.1016/j.chaos.2006.05.084.
Saltzman, B.1962. Finite amplitude free convection as an initial value problem. J. Atmos. Sci. 19, 329–341. DOI:https://doi.org/10.1175/1520-0469(1962)019<0329:FAFCAA>2.0.CO;2.
Shen, B.-W.2014. Nonlinear feedback in a five-dimensional Lorenz model. J. Atmos. Sci. 71, 1701–1723. DOI:https://doi.org/10.1175/JAS-D-13-0223.1.
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:https://doi.org/10.5194/npg-22-749-2015.
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:https://doi.org/10.5194/npg-23-189-2016.
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.
Sparrow, C.1982. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. Springer, New York, pp. 41.
Sprott, J. C.2003. Chaos and Time-Series Analysis. Oxford University Press. New York, pp. 528. The numerical method is briefly discussed on http://sprott.physics.wisc.edu/chaos/lyapexp.htm
Starossek, U.2015. A low-frequency pendulum mechanism. Mech. Mach. Theory. 83, 81–90. DOI:https://doi.org/10.1016/j.mechmachtheory.2014.08.010.
Strogatz, S. H.2015. Nonlinear Dynamics and Chaos. With Applications to Physics, Biology, Chemistry, and Engineering. Westpress view, Boulder, CO, pp. 513.
Thiffeault, J.-L. and Horton, W.1996. Energy-conserving truncations for convection with shear flow. Phys. Fluids. 8, 1715–1719. DOI:https://doi.org/10.1063/1.868956.
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:
Using the ${M}_{1}-{M}_{3}$ of Shen (2014), ψ and θ can be represented as:
Here, C_{1} and C_{2} 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 ${M}_{1}-{M}_{3}$, respectively, as follows:
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:
We define a time scale (T_{o}) as:
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:
With Equation (A9), Equation (A6) becomes:
The above equation then becomes Equation (3), $\mathrm{d}X/\mathrm{d}\tau =\sigma Y,$ when $\tau =\kappa {T}_{o}t$ is assumed.
The non-dissipative Lorenz model (3D-NLM) with no heating term is written as:
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
From Equation (B1) and (B3) (i.e. X * Equation (B1)− σ*Equation (B3)), we have
The above derivations illustrate the relationship between the X^{3} 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.
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:
The incomplete elliptic integral of the first kind, $F(\varphi ,k)$, is written as:
The elliptic functions of the Jacobi are defined as the inverse of the elliptic integral, as follows:
The elliptic functions have the following properties:
Equations (C4) and (C5) can be helpful for obtaining a solution to the equation ${\mathrm{d}}^{2}X/\mathrm{d}{\tau}^{2}+{X}^{3}/2=0$ (Equation (18a)). Using Equation (C5) and ${k}^{2}=-1$, the solution can be written as:
Through Equation (C6) and ${k}^{2}=1/2$, the solution becomes:
Here, D_{1}, D_{2}, θ_{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(\pi /2,{k}^{2}=1/2)\sim 1.8541$, which is calculated using Matlab (http://www.mathworks.com/help/matlab/ref/ellipke.html). 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, k^{2} 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., $\frac{\mathrm{d}}{\mathrm{d}\tau}X(0)=\sigma $)) leads to the following solution with θ_{2}=3K and ${D}_{2}=\sqrt{2\sigma}$:
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.
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:
Using the above with Equation (3) yields:
The term X^{4} in Equation (D4) can be shown to be equivalent to the term X^{3} 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).