We use direct statistical simulation to find the low-order statistics of the well-known dynamical system, the Lorenz63 model. Instead of accumulating statistics from numerical simulation of the dynamical system or solving the Fokker–Planck equation for the full probability distribution of the dynamical system, we directly solve the equations of motion for the low-order statistics after closing them by making several different choices for the truncation. Fixed points of the statistics are obtained either by time evolving or by iterative methods. The stability and statistical realizability of the fixed points of the statistics are analyzed, and the statistics so obtained are compared to those found by the traditional approach. Low-order statistics of the chaotic Lorenz63 system can be obtained from cumulant expansions more efficiently than by accumulation via direct numerical simulation or by solution of the Fokker–Planck equation.

Direct statistical simulation (DSS) is a method for computing the turbulent low-order statistics of a fluid dynamical system using the paradigm of statistical physics. DSS solves directly for the evolution of the cumulants of the probability density function of the fluid motion in statistical space. In this study, we demonstrate the effectiveness of DSS by solving the illustrative Lorenz63 system in a chaotic state utilizing different closure approximations. We find that the chaotic behavior of Lorenz63 can be accurately approximated by the low-order statistics up to third order. The statistical equilibrium of Lorenz63 is efficiently computed either by using the timestepping or iterative method. DSS may, therefore, be an effective scheme for computing the statistics of turbulent flows for cases where it is not possible to solve for the flow directly.

Chaotic dynamical systems can be characterized in a number of different ways, such as Lyapunov exponents, rare events, and low-order statistics. Appropriate methods should be tailored to the particular type of information that is being sought. In this paper, we explore the low-order equal-time statistics of the celebrated dynamical system, the Lorenz63 model. To find the low-order cumulants, we employ the method of direct statistical simulation (DSS) that works directly with the equations of motion for the statistics themselves, a program that Lorenz (1967) himself described in the context of numerical weather prediction. Because nonlinearities couple successive orders of the hierarchy of cumulants together, it is necessary to make closure approximations that truncate the equations of motion. We examine the reliability of several closure approximations by comparison to statistics obtained by the traditional method of time integrating the dynamical system.

The chaotic Lorenz63 model is ideal for a study of the capabilities and limitations of DSS because it is highly non-linear—more nonlinear, in fact, than many anisotropic and inhomogeneous problems in fluid dynamics, such as jets that can be viewed as perturbations to a mean-flow (see, for instance, Marston et al., 2019). The strong quadratic nonlinearity tests the limits of different closure approximations. In addition, the low dimensionality of the Lorenz63 system makes it possible to query all of the low-order equal-time statistics.

The purpose of this work is multifold. We extend the previous study of the Lorenz63 system (Allawala and Marston, 2016) in a more chaotic regime and apply the CE2.5 level of approximation that is closely related to the eddy-damped quasi-normal Markovian (EDQNM) approximation (Marston et al., 2019) to further simplify the statistical representation of the chaotic dynamical system. Furthermore, we attempt to use iterative methods for computing the optimal solution of the truncated cumulant equations that represent the statistical steady state and study the feasibility and the numerical performance of the optimization methods for determining the fixed points. The utility of finding the fixed points directly is important for DSS of more complicated PDE systems. The governing equations of motion for the cumulants are automatically generated using a symbolic manipulation with SymPy (Li et al., 2021), instead of using the general arrays to encode coefficients.

The organization of the rest of the paper is as follows. In Sec. II, we introduce the equal-time cumulants, closure approximations, and methods for finding fixed points. Section III discusses our results for the Lorenz63 attractor. Some conclusions are presented in Sec. IV.

In this section, we briefly summarize how the low-order statistics of the low-order dynamical systems can be approximated using direct statistical simulation.

Low-order dynamical systems are reduced models often used to represent the characteristic behaviors of the fluid dynamical related problems. These models are derived either via a Galerkin truncation of the partial differential equation (PDE) system (e.g., see Holmes et al., 2012 and Maasch and Saltzman, 1990) or via a normal form analysis (see, e.g., Tobias et al., 1995).

We consider a nonlinear system whose governing equation is generically represented by a set of ordinary differential equations with up to quadratic nonlinearities; the ith component reads

(1)

where the coefficient of the quadratic nonlinear interaction is given by Qijk and Lij is for the linear term. Here, we also include the possibility of stochastic forcing, fi, which is assumed to be an independent Gaussian, fiN(μi,σi2), that is introduced to synthesize the unmodeled physical processes, where μi and σi2 are the statistical mean and variance of fi, respectively (see, e.g., Allawala and Marston, 2016).

Direct statistical simulation describes the evolution of dynamical systems in probability space. In this framework, the unknown field, xi, is treated as a random variable with the associated probability density function (PDFs) represented by a series of statistics, namely, cumulants (Kendall et al., 1987). In contrast to the conventional Fokker–Planck approach, which computes the full PDFs of the state variable, xi, of the dynamical system (Pawula, 1967), DSS only solves the low-order statistics of the PDFs of xi. The random variable, xi, which consists of the coherent component, Cxi, and the non-coherent counterpart, δxi, i.e.,

(2)

is required to satisfy the Reynolds averaging rules, e.g.,

(3)

where the statistical average, noted as , is assumed to be the ensemble average in this study. In the cumulant hierarchy, the first three cumulants, Cxi, Cxixj, and Cxixjxk, are identical to the statistical central moments. The fourth and higher cumulants, unlike the second and third terms, are arbitrarily zeros for a Gaussian distribution, while the fourth centered moment does not. In this paper, we explicitly use the fourth cumulant, which is defined as

(4)

The governing equations of the cumulants can be derived directly from the dynamical equation (1) or via the Hopf functional approach (Frisch, 1995). In this study, we choose the first approach so that the first cumulant equation is obtained by taking the statistical average of the dynamical system (1) and reads

(5)

The second cumulant equation, which describes the time evolution of a two-point correlation, is defined by dtCxixm=(dtδxi)δxm+δxi(dtδxm) and written as

(6)

where δi,m is the Kronecker delta. The detailed derivation of the first three cumulant equations for a general dynamical system is documented in Eqs. (2.6)–(2.8) of Li et al. (2021), where for the quadratically nonlinear system (1) considered here, the terms with the cubic coefficients do not appear. In real-time computations, we use the existing software in Python to derive the cumulant equations. This package has been developed by the authors for studying the low-order dynamo problems (Li et al., 2021) and available at https://github.com/Kuan-Li-Math-Geo/dss_low-order.git.

The expansion of the nonlinear dynamical equation (1) leads to an infinite hierarchy of coupled equations; i.e., for the quadratic system, the ith cumulant equation always involves the i+1th cumulant. A proper statistical closure must be chosen to truncate the cumulant expansion at the lowest possible order. We truncate the cumulant equations using three different truncation rules, namely, CE2, CE2.5, and CE3. CE2 is ideal for studying the dynamical systems with the state vector satisfying or close to Gaussian distributions, where the cumulant hierarchy is truncated at second order, and all higher order terms greater than two are neglected (see, e.g., Marston et al., 2019). The CE3 rule is often applied to truncate the dynamical systems with strong asymmetry (skewness) or long tails (flatness) in PDFs, as we will see in Sec. III. In CE3 approximations, the set of cumulant equations is truncated at the third order, while the fourth order cumulant is set to zero Cxixjxkxl=0 (Orszag, 1970). The effects of the fourth order cumulants are further modeled by a diffusion process, Cxixmxn/τd, with the parameter, τd>0, known as the eddy damping parameter (Marston et al., 2019); e.g., refer to Eqs. (2.6), (2.7), and (2.10) in Li et al. (2021). The damping parameter, τd, measures the time scale of the noncoherent components of the flow resulting from their nonlinear interaction. In the highly chaotic regimes of the dynamical system, τd is expected to be a small number. The third cumulant, Cxixjxk, evolves much more rapidly in time as compared with the first and second cumulant. A further simplification of the third order cumulant equations can be made, which leads to the so-called CE2.5 approximation (Marston et al., 2019 and Allawala et al., 2020). In CE2.5, the third cumulant is determined by solving the diagnostic components of the third order equation with dtCxixmxn=0 and with all the first order cumulant setting to zero, Cxi=0, to further speed up the computation; e.g., see Eqs. (2.6), (2.7), and (2.11) in Li et al. (2021).

The low-order statistics governed by the cumulant equations evolves much smoother in space and time than the instantaneous field, xi. As forward evolving the cumulant equations, the solution settles in statistical equilibrium states, which is invariant in time. This approach will determine stable solutions to the cumulant equations. However, many of the time-invariant solutions of the cumulant system are either unstable as we integrate CE2/2.5/3 equations in time or statistically non-realizable. Realizability is satisfied only if the second cumulant satisfies the Cauchy–Schwarz inequality, Cxi2Cxj2Cxixj2. We use three different methods to study the stability and statistical realizability of the fixed points by (1) forward evolving the cumulant equations; (2) directly solving the spatial components of the (algebraic) cumulant equations, where the temporal components are set to zero; and (3) solving an inverse problem using the quasi-Newton method, the conjugate-gradient (CG) method, the trust-region method (Nocedal and Wright, 2006) or sequential quadratic programming (SLSQP) (Kraft, 1988), where the misfit functional, J, is utilized to measure the temporal variation of the cumulant equations; i.e., the misfit, J, is defined as (e.g., see Li et al., 2021)

(7)

for CE2/2.5 and CE3, respectively.

The governing equation of the Lorenz63 system (Lorenz, 1963) is obtained from the truncation of the Galerkin discretization of the atmospheric convection model at the lowest order, where the evolution equations are given by

(8)

The control parameters are the Prandtl number, Pr; the (relative) Rayleigh number, Ra; and the geometric factor, β. The unknown functions x,y, and z represent the velocity, the horizontal temperature variation, and the vertical temperature variation, respectively. The stochastic force, fx,y,z, is introduced to synthesize the unmodeled physical processes (Allawala and Marston, 2016), where fx,y,z is assumed to be independent Gaussian variable, fx,y,zN(μx,y,z,σx,y,z2), with the mean and variance given by μx,y,z and σx,y,z2.

We use the Python package to derive the cumulant equations in (8) and obtain the low-order cumulant approximations of the Lorenz63 system, where the first and the second order equations read

(9)

and

(10)

The third order equations have a complicated form and are detailed in Eq. (A1) of the  Appendix. The cumulant equations are further truncated according to CE2/2.5/3 truncation rules, respectively.

In this study, we vary the Rayleigh number, Ra, and the stochastic force, fx,y,z, to control the dynamics of the Lorenz63 system, while the Prandtl number, Pr=10, and the geometric factor, β=8/3, remain the same as those chosen by Lorenz (1963) for all cases.

We focus our study on the chaotic regime of Lorenz63. We find two physical mechanisms that are able to drive the dynamics of Lorenz63 to the chaotic states. If the Rayleigh number is supercritical, i.e., Ra is greater than a critical value, Rac, the Lorenz63 system spontaneously evolves toward the chaotic state. The chaotic solutions of Lorenz63 can also be obtained in the subcritical branch for Ra<Rac by employing a strong external stochastic force, fx,y,z. We compare the low-order statistics of the Lorenz63 system obtained by direct numerical simulation (DNS) and DSS in the chaotic regime and show the effectiveness of the low-order cumulant approximation for describing the chaos of the Lorenz63 system.

In the absence of external forcing, the Lorenz63 system (8) always evolves toward a steady state for 1<Ra<Rac. Depending on the choice of the initial condition, the steady solution of Lorenz63 converges to one of two stable attractors at

(11)

The critical Rayleigh number is determined by the Prandtl number, Pr, and the geometric factor, β, where Rac is found to be approximately 25 for Pr=10 and β=8/3. In the phase space, two stable attractors are non-connected but accessible with equal probability due to the reflective symmetry, xx and yy.

In the supercritical regime for Ra>Rac, the attractors, XF±, become unstable and the trajectory of Lorenz63 repelled by the unstable (“strange”) attractors oscillates irregularly in the phase space. Shown in Fig. 1 is the typical solution of the Lorenz63 system in the chaotic state for Ra=28. For this parameter set, the characteristic time scales of the Lorenz63 system are found to be tw3.5 and tc0.75, where tw is the average time for the trajectory to stay on one wing of the Lorenz butterfly and tc is the average time for the trajectory to orbit one of the wings. The probability distributions, P(x), P(y), and P(z), that are shown in Figs. 1(c)1(e), are non-Gaussian, which indicate the importance of the higher order statistics in the approximation of P(x), P(y), and P(z), where the green histograms stand for probability distributions, P(x), P(y), and P(z), and the dashed red curves are for the Gaussian distribution with the same mean and variance for comparison purposes. The locations of the strange attractors in x and y coordinates appear as the small peaks on the left and right hand side of P(x) and P(y); the central peak that has the largest amplitude results from the rapid oscillation of the trajectory between two strange attractors. The PDF of z, which is bimodal and unevenly distributed about the statistical mean of the vertical temperature variation, Cz, has two local maxima associated with two strange attractors.

FIG. 1.

The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for Ra=28 and fx,y,z=0, where the green histograms stand for the PDFs of x,y, and z and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).

FIG. 1.

The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for Ra=28 and fx,y,z=0, where the green histograms stand for the PDFs of x,y, and z and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).

Close modal

In the subcritical regime for 1<Ra<Rac, the stable attractor of Lorenz63 may become unstable if we increase the external stochastic force, fx,y,z. For small random forces, i.e., σx,y,z2O(1), the trajectory of Lorenz63 oscillates randomly around one of two steady solutions at XF± with the PDFs, P(x), P(y), and P(z) in Gaussian, where XF±=[±7.1,±7.1,19]; e.g., see Fig. 2. As the stochastic force further increased to σx,y,z2O(1), the attractors, XF±, become unstable. The trajectory of Lorenz63 in this regime exhibits the similar “butterfly” pattern in Figs. 3(a) and 3(b) as we observe for the supercritical case for Ra=28 and fx,y,z=0 in Figs. 1(a) and 1(b). Here, the irregular oscillation of the trajectory about XF+ and XF is due to the combined force of the nonlinear interactions and the stochastic force, fx,y,z. Shown in Fig. 3 is the chaotic solution of Lorenz63 for Ra=20 and fx,y,zN(0,2), where the PDFs, P(x), P(y), and P(z), obtained in DNS are shown in green histograms and the Gaussian distribution with the same mean and variance as for DNS are in dashed red curves. Interestingly, the PDFs of x and y appear to be the superposition of two Gaussian distributions centered at each of the attractor at ±7.1 and the PDF of z is symmetric and close to Gaussian, which is very different from P(z) in the supercritical case with very strong asymmetry. Closely analyzing the data, we find that the mean trajectory of z converges to the value, z=17.2, other than the one in the steady state at Ra1=19, which indicates that the dynamical system is in the chaotic state. The Lorenz63 system comprising two “strange attractors” behaves similarly for different Ra and fx,y,z in the chaotic state.

FIG. 2.

The illustration of the solution of Lorenz63 that oscillates about XF=[7.1,7.1,19] for Ra=20 and fx,y,zN(0,0.1), where the green histograms stand for the PDFs of x,y, and z, and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS (b)–(d).

FIG. 2.

The illustration of the solution of Lorenz63 that oscillates about XF=[7.1,7.1,19] for Ra=20 and fx,y,zN(0,0.1), where the green histograms stand for the PDFs of x,y, and z, and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS (b)–(d).

Close modal
FIG. 3.

The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for Ra=20 and fx,y,zN(0,2), where the green histograms stand for the PDFs of x,y, and z and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).

FIG. 3.

The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for Ra=20 and fx,y,zN(0,2), where the green histograms stand for the PDFs of x,y, and z and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).

Close modal

We integrate the cumulant equations, starting from the random initial conditions, forward in time to obtain the statistical equilibrium of the Lorenz63 system in the chaotic states. The results are summarized and compared with the statistics accumulated from DNS in Tables I and II. We find that for the Lorenz63 system, the skewness (third order cumulant) that quantifies the asymmetry of the probability distribution is important for accurately approximating the statistical equilibrium. Here, the chaoticity of the Lorenz63 system can be very accurately approximated by using the CE2.5 and CE3 equations, where the eddy damping parameter, τd, which parameterizes the correlation time of the interaction among the noncoherent components of the system, must be introduced to stabilize the numerical computation for all cases. The computations of the CE2.5/3 equation are numerically stable for all cases for an eddy damping parameter in the range, τd[103,101], and the optimal τd is found to be O(101), which is approximately ten times smaller than the characteristic time scales of the Lorenz63 system, tw and tc, where tw is the average time for the trajectory to stay on one wing of the Lorenz butterfly and tc is the average time for the trajectory to orbit one of the wings.

TABLE I.

The list of the low-order cumulants of the Lorenz63 system in the chaotic state for the Rayleigh number in the supercritical regime for Ra = 28 for fx,y,z = 0 and fx,y,zN(5,10), where the case labeled as “CE3s” is the solution of the CE3 equation without the temporal components of third order equations.

Rafx,y,zτdCxCyCzCxxCxyCxzCyyCyzCzz
DNS 28  25.10 67.84 66.95 84.02 58.73 
CE2.5   1/40 23.84 63.58 63.58 62.84 75.54 
CE2.5   1/30 23.84 66.21 66.21 65.93 54.01 
CE3   1/30 24.54 65.45 65.45 71.52 58.00 
CE3   1/200 26.38 70.34 70.34 70.80 16.23 
CE3s   1/40 24.92 66.45 66.45 70.99 50.13 
DNS 28 N(5,10)  3.03 2.53 24.84 54.10 53.56 12.04 67.39 4.74 54.82 
CE2.5   1/20.5 3.09 2.59 24.71 53.91 52.91 12.56 57.78 −5.03 63.22 
CE3   1/15 3.00 2.50 24.84 54.76 53.76 11.98 66.10 −2.90 57.68 
Rafx,y,zτdCxCyCzCxxCxyCxzCyyCyzCzz
DNS 28  25.10 67.84 66.95 84.02 58.73 
CE2.5   1/40 23.84 63.58 63.58 62.84 75.54 
CE2.5   1/30 23.84 66.21 66.21 65.93 54.01 
CE3   1/30 24.54 65.45 65.45 71.52 58.00 
CE3   1/200 26.38 70.34 70.34 70.80 16.23 
CE3s   1/40 24.92 66.45 66.45 70.99 50.13 
DNS 28 N(5,10)  3.03 2.53 24.84 54.10 53.56 12.04 67.39 4.74 54.82 
CE2.5   1/20.5 3.09 2.59 24.71 53.91 52.91 12.56 57.78 −5.03 63.22 
CE3   1/15 3.00 2.50 24.84 54.76 53.76 11.98 66.10 −2.90 57.68 
TABLE II.

The list of the low-order cumulants of the Lorenz63 system in the chaotic state with the Rayleigh number in the subcritical regime for Ra = 20 in the presence of noise, fx,y,z ∼ N(0, 2), in DNS and DSS.

Rafx,y,zτdCxCyCzCxxCxyCxzCyyCyzCzz
DNS 20 N(0,2)  16.68 45.06 44.46 53.89 44.74 
CE2  N(0,2)  19.04 50.96 50.76 50.93 0.75 
CE2.5  N(0,2) 1/50 17.25 46.20 46.00 46.02 31.66 
CE3  N(0,2) 1/30 17.14 45.90 45.70 48.30 32.46 
Rafx,y,zτdCxCyCzCxxCxyCxzCyyCyzCzz
DNS 20 N(0,2)  16.68 45.06 44.46 53.89 44.74 
CE2  N(0,2)  19.04 50.96 50.76 50.93 0.75 
CE2.5  N(0,2) 1/50 17.25 46.20 46.00 46.02 31.66 
CE3  N(0,2) 1/30 17.14 45.90 45.70 48.30 32.46 

For all test cases with the Rayleigh number in the supercritical regime, Ra=28, listed in Table I, the CE2.5 system, which assumes the cumulant equation for the third order to be diagnostic in time and as a function of the second cumulant only, gives accurate approximations of the low-order statistics for the Lorenz63 system as compared with the solutions in CE3. We observe that if we set τd to be too small, the second cumulant may become less accurate, while the first cumulants remain very accurate; e.g., see Czz for τd=1/200. We study the dynamics of the Lorenz63 system in the presence of Gaussian noise for fx,y,zN(5,10). The stochastic force has a strong mean value, μi=5, which breaks the symmetry of the probability distribution for the x and y components of the system about zero. We also attempt to regularize the second cumulants of the Lorenz63 system to be closer to Gaussian by setting the variance of fx,y,z to be large; e.g, σ2=10. We find that the CE2.5 and CE3 equations also approximate the low-order cumulants accurately, while the CE2 approximation never converges. Similar observations are made in Allawala and Marston (2016) for the Lorenz63 system with less chaoticity.

For the Lorenz63 system with the Rayleigh number in the subcritical regime, the stochastic force must be applied in order to drive the system to the chaotic state; e.g., see Table II for Ra=20 and fx,y,zN(0,2). We find that the CE2.5/3 approximations also perform very well as compared with the results obtained in DNS. Perhaps, it is of interest to mention that for this parameter set, the CE2 approximation successfully converges to a fixed point but very poorly estimates the statistical equilibrium; e.g., the first cumulant, Cz, fails to converge to the mean trajectory of z in the chaotic state but to the one in the steady state at Ra1=19; the second cumulant, Czz, is inaccurately determined due to the lack of the knowledge of the third cumulant, Cxyz.

Although timestepping allows the access of the stable solutions of the cumulant equations, other fixed points are possible solutions as discussed earlier. Here, we assess the effectiveness of various methods for accessing these fixed points.

The fixed points of Lorenz63 can be solved directly via the symbolic packages of Mathematica or Python, where the fixed points are assumed invariant in time. We obtain 7 statistically realizable fixed points out of 41 roots of CE3 equations of Lorenz63 defined in (9), (10), and (A1) for Ra=28, fx,y,z=0, and τd1=20. The fixed point that corresponds to the strange attractors is the only stable one in time. As increasing the dynamics of Lorenz63 from the steady state to the chaotic state by increasing the Rayleigh number, the stable attractors evolve to the “strange” attractors; see Fig. 4(a) for the continuation of the attractor as a function of Rayleigh number for fx,y,z=0. If τd is small, the third order equations are statistically insignificant in the CE3 equations and the CE3 approximation is reduced to CE2; see Fig. 4(b) for the continuation of the attractor as a function of the eddy damping parameter, τd. We also observe that other fixed points are sensitive to the choice of τd. Shown in Fig. 4(c) is the solution of an unstable attractor of Lorenz63 as a function of τd, where the second order cumulant, Czz>0, is statistically realizable for τd>0.1 and becomes non-realizable (Czz<0) when τd<0.1.

FIG. 4.

The illustration of the continuation of the attractor as a function of Rayleigh number, Ra, from Ra=20 to 30 in (a) and as a function of the eddy damping parameter, τd, in (b); the continuation of an unstable attractor as a function of τd in (c), where the unstable fixed point becomes statistically realizable as τd>0.1. (a) The fixed point of the attractor as a function of Ra. (b) The fixed point of the attractor as a function of τd. (c) An unstable fixed point as a function of τd.

FIG. 4.

The illustration of the continuation of the attractor as a function of Rayleigh number, Ra, from Ra=20 to 30 in (a) and as a function of the eddy damping parameter, τd, in (b); the continuation of an unstable attractor as a function of τd in (c), where the unstable fixed point becomes statistically realizable as τd>0.1. (a) The fixed point of the attractor as a function of Ra. (b) The fixed point of the attractor as a function of τd. (c) An unstable fixed point as a function of τd.

Close modal

We also study gradient based optimization methods for computing the fixed point of the CE2.5 approximations in the presence of noise. For σx,y,z2=2, the dynamical system can be accurately approximated by timestepping the CE2.5 equations, and the optimal solution for J converges to the same fixed point as obtained via timestepping. Shown in Fig. 5 is the convergence to the optimal solution for the cumulants of the CE2.5 approximation of the Lorenz63 in the chaotic state for Ra=38 and fx,y,zN(0,2) via the CG method, where J is normalized by its value at the first iteration. The initial guess is taken from the fixed point for Ra=28 and fx,y,zN(0,2). It is of importance to note that the misfit, J, of Lorenz63 is not convex everywhere. If the initial guess is randomly chosen, the optimal solution of the CE2.5 system may converge to other fixed points. However, we observe that for CE2.5 approximation, the stable fixed point can always be found by continuing the stable fixed point from the solutions of the nearby control parameter. For the first O(300) iterations, the CG method slightly outperforms the timestepping method in terms of the convergence rate, where the misfit, J, reduces by a factor of 107 and the optimal time step, dt102, is used for the numerical integration. After approximately 300 iterations, the convergence rate of J of CG flattens out, but the rate for timestepping remains constant. For this case, the quasi-Newton method performs similarly as CG in terms of the accuracy and the convergence rate. We note that calculating the “downhill” direction for minimizing the misfit, J, that is directly computed via the symbolic differentiation is almost computationally as expensive as for evolving the cumulant equation one time step forward.

FIG. 5.

The stable fixed point is found by minimizing J via the CG and timestepping method for the CE2.5 approximation of the Lorenz63 system for Ra=30 and fx,y,zN(0,2), where (a) and (b) show the path of the cumulants, Cz and Czz, as the solution converges to J=0 and (c) illustrates the convergence of J as a function of the number of iterations, n. The initial guess and the terminal solution of the minimization are shown as green and red dots in (a) and (b). The optimal solution converges to the stable time-invariant solution of CE2.5 via the timestepping method.

FIG. 5.

The stable fixed point is found by minimizing J via the CG and timestepping method for the CE2.5 approximation of the Lorenz63 system for Ra=30 and fx,y,zN(0,2), where (a) and (b) show the path of the cumulants, Cz and Czz, as the solution converges to J=0 and (c) illustrates the convergence of J as a function of the number of iterations, n. The initial guess and the terminal solution of the minimization are shown as green and red dots in (a) and (b). The optimal solution converges to the stable time-invariant solution of CE2.5 via the timestepping method.

Close modal

We find that it is difficult to find the stable fixed point for CE3 approximations via the gradient based method. The optimization either converges slowly to unstable or non-realizable fixed points or is trapped by the local minima that are introduced by the third order cumulants.

As an example, we use CG to continue the stable fixed point of Lorenz63 from Ra=28 (found by timestepping) to Ra=29 with ΔRa=1 for the same stochastic force level [fx,y,zN(0,2)]. The results are in Fig. 6, which shows the path of the low-order cumulant, Cz and Czz, as J converges to zero and the convergence of J as the number of iterations, where the green and purple dots in Figs. 6(a) and 6(b) represent the initial guess and the terminal solution of the optimization and the red dots are for the solution found by timestepping. Disappointingly, the optimization goes into the wrong direction and is trapped by the local minima, where the misfit saturates at the level of JO(104). Similar performance is observed for different Ra and ΔRa with and without the stochastic force, fx,y,z.

FIG. 6.

The migration of Cz and Czz from Ra=28 to 29 via the timestepping and CG method for CE3 equations in (a) and (b) and the convergence of J as the number of time steps/iterations in (c), where the stochastic force is zero; i.e., fx,y,z=0. The green dots in (a) and (b) stand for the initial condition/guess for the timestepping/CG method, the red dots are for the time-invariant solution found by the timestepping method, and the purple ones are for the terminal solution of CG.

FIG. 6.

The migration of Cz and Czz from Ra=28 to 29 via the timestepping and CG method for CE3 equations in (a) and (b) and the convergence of J as the number of time steps/iterations in (c), where the stochastic force is zero; i.e., fx,y,z=0. The green dots in (a) and (b) stand for the initial condition/guess for the timestepping/CG method, the red dots are for the time-invariant solution found by the timestepping method, and the purple ones are for the terminal solution of CG.

Close modal

The stable, time-invariant solution of the CE2.5/3 approximations of Lorenz63 in the chaotic state can always be obtained by the timestepping method. An interesting calculation is to determine the path of approach to a fixed point using a nearby solution as an initial guess. Shown in Fig. 7 are the paths of the low-order cumulants, Cz and Czz, as J approaches zero and the convergence of J as the number of time steps for the CE2.5/3 approximation of Lorenz63. Here, the initial state is that calculated at Ra=28 for fx,y,z=0 and fx,y,zN(0,10), and we show three different cases with ΔRa=1,5, and 10. Interestingly, the path of the cumulant, e.g., Cz and Czz, demonstrates a strong self-similarity as we evolve the CE3 equations in time to reduce J for different ΔRa. The form of the approach to the fixed point suggests that taking reasonably large values of ΔRa is the best strategy for continuation of solutions. For the cumulant system with a non-negligible third order cumulant for fx,y,z=0, the misfit, J, is never found convex as the solution of the cumulant equations converges to the stable fixed point. The same observations have been made for all other test cases. We speculate that this is the reason why the minimization method fails to optimize the CE3 system.

FIG. 7.

The illustration of the fixed points of the CE3 approximation of the Lorenz63 system in the chaotic state and for different Ra found by the timestepping method, where (a) and (b) and (d) and (e) show the path of Cz and Czz as J tends to zero and (c) and (f) show the exponential convergence of J as the number of iterations.

FIG. 7.

The illustration of the fixed points of the CE3 approximation of the Lorenz63 system in the chaotic state and for different Ra found by the timestepping method, where (a) and (b) and (d) and (e) show the path of Cz and Czz as J tends to zero and (c) and (f) show the exponential convergence of J as the number of iterations.

Close modal

We note here that the quasi-Newton method performs similarly to the conjugate-gradient method (CG) in terms of the numerical accuracy and the convergence rate for all test cases, while other minimization methods, such as SLSQP or the trust-region, have some limited power for computing the stable fixed point of the cumulant equations.

In this paper, we implement direct statistical simulation to study the simplified highly nonlinear dynamical system, Lorenz63, and apply DSS to study the statistical behavior of this system in the chaotic regime without and with external random forcing.

We find that the CE2.5 and CE3 approximations are sufficiently accurate to describe the long-term statistical evolution of these systems, though the probability distributions of these systems are either strongly asymmetric or have longer tails than a Gaussian distribution. The mean trajectory is very accurately determined by the cumulant equations with a maximum truncation at third order in the cumulants. There is a relatively small error of less than 1% in the mean statistics compared with those obtained via the traditional approach of long time averaging of direct numerical simulation. The interactions between the coherent and non-coherent components of the dynamics are also accurately quantified by the second order terms with relative errors of less than 20%. For both of the CE2.5 and CE3 closures, a single eddy damping parameter, τd, is introduced to approximate the correlation time among the non-coherent components of the dynamics and to stabilize the numerical integration of DSS equations. The optimal τd is found to be approximately ten times smaller than the characteristic time scale of the Lorenz63 system, tw and tc, for τdO(101). This result is consistent with Allawala and Marston (2016) for the Lorenz63 system in the less chaotic regime.

We also attempt to directly access the fixed points of the cumulant equations of Lorenz63 systems. For this dynamical system, all time-invariant solutions can be solved symbolically. There are seven roots among 41 solutions of Lorenz63 that are statistically realizable for τd=1/20, but the fixed point corresponding to the strange attractor is the only one stable in time. The same symbolic technique cannot be applied to discover the fixed points of systems with large numbers of degree of freedoms due to the rapid increase in computational complexity. We also find that the timestepping method and the gradient methods converge to the statistical equilibrium exponentially at comparable rates, except in the case of CE3 for which the gradient approach never converges.

Cumulant expansions offer a way to find the low-order statistics of dynamical systems directly. By contrast, the Fokker–Planck equation captures all of the equal-time statistics, but it is usually too computationally intensive to solve. In this study, we have demonstrated the effectiveness of the cumulant equations for describing the statistical evolution of the highly nonlinear dynamical system by solving the DSS equations of the Lorenz63 models that represent dimensionally reduced fluids. We note that turbulent fluid dynamical systems with a huge number of freedoms are governed by partial differential equations that are typically less chaotic and, therefore, is expected to be more accurately approximated by the CE2/2.5 equations. Furthermore, we believe that it is of great interest to investigate the use of gradient based methods in conjunction with timestepping the DSS equations to find the steady-state statistics of such fluids.

This work was supported in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. D5S-DLV-786780) and by a grant from the Simons Foundation (Grant No. 662962, GF).

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

The third order cumulant equations for the CE3 approximation of the Lorenz63 system consist of ten equations; i.e.,

(A1)
1
Allawala
,
A.
and
Marston
,
J. B.
, “
Statistics of the stochastically forced Lorenz attractor by the Fokker-Planck equation and cumulant expansions
,”
Phys. Rev. E
94
,
052218
(
2016
).
2
Allawala
,
A.
,
Tobias
,
S. M.
, and
Marston
,
J. B.
, “
Dimensional reduction of direct statistical simulation
,”
J. Fluid Mech.
898
,
A21
(
2020
).
3
Frisch
,
U.
,
Turbulence: The Legacy of A. N. Kolmogorov
(
Cambridge University Press
,
1995
).
4
Holmes
,
P.
,
Lumley
,
J. L.
,
Berkooz
,
G.
, and
Rowley
,
C. W.
, Proper Orthogonal Decomposition, Cambridge Monographs on Mechanics (Cambridge University Press, 2012), 2nd ed., pp. 68–105.
5
Kendall
,
M. G.
,
Stuart
,
A.
, and
Ord
,
J. K.
,
Kendall’s Advanced Theory of Statistics
(
Oxford University Press, Inc.
,
New York
,
1987
), ISBN: 0195205618.
6
Kraft
,
D. A.
, “A software package for sequential quadratic programming,” Technical Report No. DFVLR-FB 88-28, Institute for Flight Mechanics, DLR German Aerospace Center, Cologne, Germany, 1988), pp. 1420–1425.
7
Li
,
K.
,
Marston
,
J. B.
, and
Tobias
,
S. M.
, “
Direct statistical simulation of low-order dynamo systems
,”
Proc. R. Soc. A
477
,
0427
(
2021
).
8
Lorenz
,
E. N.
, “
Deterministic nonperiodic flow
,”
J. Atmos. Sci.
20
(
2
),
130
141
(
1963
).
9
Lorenz
,
E. N.
,
The Nature and Theory of the General Circulation of the Atmosphere
(
World Meteorological Organization
,
Geneva
,
1967
), Vol. 218.
10
Maasch
,
K. A.
and
Saltzman
,
B.
, “
A low-order dynamical model of global climatic variability over the full pleistocene
,”
J. Geophys. Res.: Atmos.
95
(
D2
),
1955
1963
, https://doi.org/10.1029/JD095iD02p01955 (
1990
).
11
Marston
,
J. B.
,
Qi
,
W.
, and
Tobias
,
S. M.
, “Direct statistical simulation of a jet,” in Zonal Jets: Phenomenology, Genesis, and Physics, edited by B. Galperin and P. Read (Cambridge University Press, 2019).
12
Nocedal
,
J.
and
Wright
,
S. J.
,
Numerical Optimization
, 2nd ed. (
Springer
,
New York
,
2006
).
13
Orszag
,
S. A.
, “
Analytical theories of turbulence
,”
J. Fluid Mech.
41
(
2
),
363
386
(
1970
).
14
Pawula
,
R. F.
, “
Approximation of the linear Boltzmann equation by the Fokker-Planck equation
,”
Phys. Rev.
162
,
186
188
(
1967
).
15
Tobias
,
S. M.
,
Weiss
,
N. O.
, and
Kirk
,
V.
, “
Chaotically modulated stellar dynamos
,”
MNRAS
273
,
1150
1166
(
1995
).