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.
I. INTRODUCTION
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.
II. CUMULANT STATISTICS OF LOW-ORDER DYNAMICAL SYSTEMS
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 th component reads
where the coefficient of the quadratic nonlinear interaction is given by and is for the linear term. Here, we also include the possibility of stochastic forcing, , which is assumed to be an independent Gaussian, , that is introduced to synthesize the unmodeled physical processes, where and are the statistical mean and variance of , 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, , 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, , of the dynamical system (Pawula, 1967), DSS only solves the low-order statistics of the PDFs of . The random variable, , which consists of the coherent component, , and the non-coherent counterpart, , i.e.,
is required to satisfy the Reynolds averaging rules, e.g.,
where the statistical average, noted as , is assumed to be the ensemble average in this study. In the cumulant hierarchy, the first three cumulants, , , and , 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
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
The second cumulant equation, which describes the time evolution of a two-point correlation, is defined by and written as
where 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 th cumulant equation always involves the th 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 (Orszag, 1970). The effects of the fourth order cumulants are further modeled by a diffusion process, , with the parameter, , 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, , 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, is expected to be a small number. The third cumulant, , 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 and with all the first order cumulant setting to zero, , 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, . 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, . 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, , is utilized to measure the temporal variation of the cumulant equations; i.e., the misfit, , is defined as (e.g., see Li et al., 2021)
for CE2/2.5 and CE3, respectively.
III. LORENZ63 SYSTEM
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
The control parameters are the Prandtl number, ; the (relative) Rayleigh number, ; and the geometric factor, . The unknown functions , and represent the velocity, the horizontal temperature variation, and the vertical temperature variation, respectively. The stochastic force, , is introduced to synthesize the unmodeled physical processes (Allawala and Marston, 2016), where is assumed to be independent Gaussian variable, , with the mean and variance given by and .
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
and
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, , and the stochastic force, , to control the dynamics of the Lorenz63 system, while the Prandtl number, , and the geometric factor, , 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., is greater than a critical value, , the Lorenz63 system spontaneously evolves toward the chaotic state. The chaotic solutions of Lorenz63 can also be obtained in the subcritical branch for by employing a strong external stochastic force, . 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.
A. Lorenz63 system in the chaotic regime
In the absence of external forcing, the Lorenz63 system (8) always evolves toward a steady state for . Depending on the choice of the initial condition, the steady solution of Lorenz63 converges to one of two stable attractors at
The critical Rayleigh number is determined by the Prandtl number, , and the geometric factor, , where is found to be approximately for and . In the phase space, two stable attractors are non-connected but accessible with equal probability due to the reflective symmetry, and .
In the supercritical regime for , the attractors, , 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 . For this parameter set, the characteristic time scales of the Lorenz63 system are found to be and , where is the average time for the trajectory to stay on one wing of the Lorenz butterfly and is the average time for the trajectory to orbit one of the wings. The probability distributions, , , and , 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 , , and , where the green histograms stand for probability distributions, , , and , 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 and coordinates appear as the small peaks on the left and right hand side of and ; the central peak that has the largest amplitude results from the rapid oscillation of the trajectory between two strange attractors. The PDF of , which is bimodal and unevenly distributed about the statistical mean of the vertical temperature variation, , has two local maxima associated with two strange attractors.
The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for and , where the green histograms stand for the PDFs of , and and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).
The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for and , where the green histograms stand for the PDFs of , and and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).
In the subcritical regime for , the stable attractor of Lorenz63 may become unstable if we increase the external stochastic force, . For small random forces, i.e., , the trajectory of Lorenz63 oscillates randomly around one of two steady solutions at with the PDFs, , , and in Gaussian, where ; e.g., see Fig. 2. As the stochastic force further increased to , the attractors, , 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 and in Figs. 1(a) and 1(b). Here, the irregular oscillation of the trajectory about and is due to the combined force of the nonlinear interactions and the stochastic force, . Shown in Fig. 3 is the chaotic solution of Lorenz63 for and , where the PDFs, , , and , 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 and appear to be the superposition of two Gaussian distributions centered at each of the attractor at and the PDF of is symmetric and close to Gaussian, which is very different from in the supercritical case with very strong asymmetry. Closely analyzing the data, we find that the mean trajectory of converges to the value, , other than the one in the steady state at , which indicates that the dynamical system is in the chaotic state. The Lorenz63 system comprising two “strange attractors” behaves similarly for different and in the chaotic state.
The illustration of the solution of Lorenz63 that oscillates about for and , where the green histograms stand for the PDFs of , and , and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS (b)–(d).
The illustration of the solution of Lorenz63 that oscillates about for and , where the green histograms stand for the PDFs of , and , and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS (b)–(d).
The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for and , where the green histograms stand for the PDFs of , and and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).
The illustration of the trajectory, the time series, and the PDFs of Lorenz63 in the chaotic state for and , where the green histograms stand for the PDFs of , and and the dashed red curves are for the Gaussian distribution with the same mean and variance as for DNS in (c)–(e).
B. Direct statistical simulation of Lorenz63 in the chaotic states
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, , 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, , and the optimal is found to be , which is approximately ten times smaller than the characteristic time scales of the Lorenz63 system, and , where is the average time for the trajectory to stay on one wing of the Lorenz butterfly and is the average time for the trajectory to orbit one of the wings.
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 , where the case labeled as “CE3s” is the solution of the CE3 equation without the temporal components of third order equations.
. | Ra . | fx,y,z . | τd . | Cx . | Cy . | Cz . | Cxx . | Cxy . | Cxz . | Cyy . | Cyz . | Czz . |
---|---|---|---|---|---|---|---|---|---|---|---|---|
DNS | 28 | 0 | 0 | 0 | 25.10 | 67.84 | 66.95 | 0 | 84.02 | 0 | 58.73 | |
CE2.5 | 1/40 | 0 | 0 | 23.84 | 63.58 | 63.58 | 0 | 62.84 | 0 | 75.54 | ||
CE2.5 | 1/30 | 0 | 0 | 23.84 | 66.21 | 66.21 | 0 | 65.93 | 0 | 54.01 | ||
CE3 | 1/30 | 0 | 0 | 24.54 | 65.45 | 65.45 | 0 | 71.52 | 0 | 58.00 | ||
CE3 | 1/200 | 0 | 0 | 26.38 | 70.34 | 70.34 | 0 | 70.80 | 0 | 16.23 | ||
CE3s | 1/40 | 0 | 0 | 24.92 | 66.45 | 66.45 | 0 | 70.99 | 0 | 50.13 | ||
DNS | 28 | 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 |
. | Ra . | fx,y,z . | τd . | Cx . | Cy . | Cz . | Cxx . | Cxy . | Cxz . | Cyy . | Cyz . | Czz . |
---|---|---|---|---|---|---|---|---|---|---|---|---|
DNS | 28 | 0 | 0 | 0 | 25.10 | 67.84 | 66.95 | 0 | 84.02 | 0 | 58.73 | |
CE2.5 | 1/40 | 0 | 0 | 23.84 | 63.58 | 63.58 | 0 | 62.84 | 0 | 75.54 | ||
CE2.5 | 1/30 | 0 | 0 | 23.84 | 66.21 | 66.21 | 0 | 65.93 | 0 | 54.01 | ||
CE3 | 1/30 | 0 | 0 | 24.54 | 65.45 | 65.45 | 0 | 71.52 | 0 | 58.00 | ||
CE3 | 1/200 | 0 | 0 | 26.38 | 70.34 | 70.34 | 0 | 70.80 | 0 | 16.23 | ||
CE3s | 1/40 | 0 | 0 | 24.92 | 66.45 | 66.45 | 0 | 70.99 | 0 | 50.13 | ||
DNS | 28 | 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 |
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.
. | Ra . | fx,y,z . | τd . | Cx . | Cy . | Cz . | Cxx . | Cxy . | Cxz . | Cyy . | Cyz . | Czz . |
---|---|---|---|---|---|---|---|---|---|---|---|---|
DNS | 20 | 0 | 0 | 16.68 | 45.06 | 44.46 | 0 | 53.89 | 0 | 44.74 | ||
CE2 | 0 | 0 | 19.04 | 50.96 | 50.76 | 0 | 50.93 | 0 | 0.75 | |||
CE2.5 | 1/50 | 0 | 0 | 17.25 | 46.20 | 46.00 | 0 | 46.02 | 0 | 31.66 | ||
CE3 | 1/30 | 0 | 0 | 17.14 | 45.90 | 45.70 | 0 | 48.30 | 0 | 32.46 |
. | Ra . | fx,y,z . | τd . | Cx . | Cy . | Cz . | Cxx . | Cxy . | Cxz . | Cyy . | Cyz . | Czz . |
---|---|---|---|---|---|---|---|---|---|---|---|---|
DNS | 20 | 0 | 0 | 16.68 | 45.06 | 44.46 | 0 | 53.89 | 0 | 44.74 | ||
CE2 | 0 | 0 | 19.04 | 50.96 | 50.76 | 0 | 50.93 | 0 | 0.75 | |||
CE2.5 | 1/50 | 0 | 0 | 17.25 | 46.20 | 46.00 | 0 | 46.02 | 0 | 31.66 | ||
CE3 | 1/30 | 0 | 0 | 17.14 | 45.90 | 45.70 | 0 | 48.30 | 0 | 32.46 |
For all test cases with the Rayleigh number in the supercritical regime, , 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 to be too small, the second cumulant may become less accurate, while the first cumulants remain very accurate; e.g., see for . We study the dynamics of the Lorenz63 system in the presence of Gaussian noise for . The stochastic force has a strong mean value, , which breaks the symmetry of the probability distribution for the and 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 to be large; e.g, . 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 and . 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, , fails to converge to the mean trajectory of in the chaotic state but to the one in the steady state at ; the second cumulant, , is inaccurately determined due to the lack of the knowledge of the third cumulant, .
C. Fixed points for Lorenz63 in the chaotic state
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 statistically realizable fixed points out of roots of CE3 equations of Lorenz63 defined in (9), (10), and (A1) for , , and . 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 . If 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, . We also observe that other fixed points are sensitive to the choice of . Shown in Fig. 4(c) is the solution of an unstable attractor of Lorenz63 as a function of , where the second order cumulant, , is statistically realizable for and becomes non-realizable () when .
The illustration of the continuation of the attractor as a function of Rayleigh number, , from to in (a) and as a function of the eddy damping parameter, , in (b); the continuation of an unstable attractor as a function of in (c), where the unstable fixed point becomes statistically realizable as . (a) The fixed point of the attractor as a function of . (b) The fixed point of the attractor as a function of . (c) An unstable fixed point as a function of .
The illustration of the continuation of the attractor as a function of Rayleigh number, , from to in (a) and as a function of the eddy damping parameter, , in (b); the continuation of an unstable attractor as a function of in (c), where the unstable fixed point becomes statistically realizable as . (a) The fixed point of the attractor as a function of . (b) The fixed point of the attractor as a function of . (c) An unstable fixed point as a function of .
We also study gradient based optimization methods for computing the fixed point of the CE2.5 approximations in the presence of noise. For , the dynamical system can be accurately approximated by timestepping the CE2.5 equations, and the optimal solution for 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 and via the CG method, where is normalized by its value at the first iteration. The initial guess is taken from the fixed point for and . It is of importance to note that the misfit, , 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 iterations, the CG method slightly outperforms the timestepping method in terms of the convergence rate, where the misfit, , reduces by a factor of and the optimal time step, , is used for the numerical integration. After approximately iterations, the convergence rate of 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, , that is directly computed via the symbolic differentiation is almost computationally as expensive as for evolving the cumulant equation one time step forward.
The stable fixed point is found by minimizing via the CG and timestepping method for the CE2.5 approximation of the Lorenz63 system for and , where (a) and (b) show the path of the cumulants, and , as the solution converges to and (c) illustrates the convergence of as a function of the number of iterations, . 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.
The stable fixed point is found by minimizing via the CG and timestepping method for the CE2.5 approximation of the Lorenz63 system for and , where (a) and (b) show the path of the cumulants, and , as the solution converges to and (c) illustrates the convergence of as a function of the number of iterations, . 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.
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 (found by timestepping) to with for the same stochastic force level []. The results are in Fig. 6, which shows the path of the low-order cumulant, and , as converges to zero and the convergence of 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 . Similar performance is observed for different and with and without the stochastic force, .
The migration of and from to via the timestepping and CG method for CE3 equations in (a) and (b) and the convergence of as the number of time steps/iterations in (c), where the stochastic force is zero; i.e., . 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.
The migration of and from to via the timestepping and CG method for CE3 equations in (a) and (b) and the convergence of as the number of time steps/iterations in (c), where the stochastic force is zero; i.e., . 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.
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, and , as approaches zero and the convergence of as the number of time steps for the CE2.5/3 approximation of Lorenz63. Here, the initial state is that calculated at for and , and we show three different cases with , and . Interestingly, the path of the cumulant, e.g., and , demonstrates a strong self-similarity as we evolve the CE3 equations in time to reduce for different . The form of the approach to the fixed point suggests that taking reasonably large values of is the best strategy for continuation of solutions. For the cumulant system with a non-negligible third order cumulant for , the misfit, , 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.
The illustration of the fixed points of the CE3 approximation of the Lorenz63 system in the chaotic state and for different found by the timestepping method, where (a) and (b) and (d) and (e) show the path of and as tends to zero and (c) and (f) show the exponential convergence of as the number of iterations.
The illustration of the fixed points of the CE3 approximation of the Lorenz63 system in the chaotic state and for different found by the timestepping method, where (a) and (b) and (d) and (e) show the path of and as tends to zero and (c) and (f) show the exponential convergence of as the number of iterations.
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.
IV. CONCLUSION
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, , 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 is found to be approximately ten times smaller than the characteristic time scale of the Lorenz63 system, and , for . 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 , 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.
ACKNOWLEDGMENTS
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).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
DATA AVAILABILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.
APPENDIX: THE THIRD ORDER CUMULANT EQUATIONS OF LORENZ63
The third order cumulant equations for the CE3 approximation of the Lorenz63 system consist of ten equations; i.e.,