Equation-free implementation of statistical moment closures

TL;DR

This paper introduces a numerical scheme that leverages equation-free methods to implement statistical moment closures for complex nonlinear systems, enabling efficient simulation without explicit closure equations.

Contribution

It develops a general, flexible approach for numerically integrating closure dynamics in large-scale nonlinear systems without requiring closed-form equations.

Findings

Successfully applied to a nonlinear stochastic PDE

Demonstrates computational advantages over direct simulation

Applicable to systems without scale separation

Abstract

We present a general numerical scheme for the practical implementation of statistical moment closures suitable for modeling complex, large-scale, nonlinear systems. Building on recently developed equation-free methods, this approach numerically integrates the closure dynamics, the equations of which may not even be available in closed form. Although closure dynamics introduce statistical assumptions of unknown validity, they can have significant computational advantages as they typically have fewer degrees of freedom and may be much less stiff than the original detailed model. The closure method can in principle be applied to a wide class of nonlinear problems, including strongly-coupled systems (either deterministic or stochastic) for which there may be no scale separation. We demonstrate the equation-free approach for implementing entropy-based Eyink-Levermore closures on a nonlinear stochastic partial differential equation.

Figures

Equations

$$\dot{\mathbf{X}}=\mathbf{U}(\mathbf{X}(t),\mathbf{N}(t),t)$$

$${\mathbf{X}}_{t+1}={\mathbf{U}}_t({\mathbf{X}}_t,{\mathbf{N}}_t)$$

$$d\mathbf{X}=\mathbf{U}(\mathbf{X},t)dt+\sqrt{2}\mathbf{S}(\mathbf{X},t)d\mathbf{W}(t).$$

$${\partial }_{t}P={\mathcal{L}}^{\ast }(t)P,$$

$${\mathcal{L}}^{\ast }(t)\psi (\mathbf{X})=-{\nabla }_{\mathbf{X}}\cdot (\mathbf{U}(\mathbf{X},t)\psi (\mathbf{X}))+{\nabla }_{\mathbf{X}}^2:(\mathbf{D}(\mathbf{X},t)\psi (\mathbf{X})).$$

$${\partial }_{t}P+{\nabla }_{\mathbf{X}}\cdot (\mathbf{U}P)={\nabla }_{\mathbf{X}}^2:(\mathbf{D}P)$$

$$\displaystyle{\mbox{\boldmath$\mu$}}(t)=\int{\mbox{\boldmath$\xi$}}({\bf X},t)P({\bf X},t)d{\bf X}$$

$$\displaystyle{\dot{\mbox{\boldmath$\mu$}}(t)}=\int{\dot{\mbox{\boldmath$\xi$}}}({\bf X},t)P({\bf X},t)d{\bf X},$$

$$\displaystyle{\dot{\mbox{\boldmath$\xi$}}}({\bf X},t)=\partial_{t}{\mbox{\boldmath$\xi$}}({\bf X},t)+{\cal L}(t){\mbox{\boldmath$\xi$}}({\bf X},t)$$

$$\displaystyle{\dot{\mbox{\boldmath$\mu$}}}(t)={\bf V}({\mbox{\boldmath$\mu$}},t)\equiv\int{\dot{\mbox{\boldmath$\xi$}}}({\bf X},t)P({\bf X},t,{\mbox{\boldmath$\mu$}})d{\bf X}.$$

$$\displaystyle{\mbox{\boldmath$\mu$}}(t+\Delta t)={\mbox{\boldmath$\mu$}}(t)+\Delta t\left[\frac{{\mbox{\boldmath$\mu$}}(t+\delta t)-{\mbox{\boldmath$\mu$}}(t)}{\delta t}\right].$$

$$\frac{\partial \varphi (x,t)}{\partial t}=D\Delta \varphi (x,t)-V^{\prime }(\varphi (x,t))+\eta (x,t)$$

$$V(\varphi )=\frac{1}{2}{\varphi }^2+\frac{1}{4}{\varphi }^4$$

$$H[\varphi ]=\int [\frac{1}{2}D|\nabla \varphi (x){|}^2+V(\varphi (x))]dx.$$

$$\frac{\partial ⟨\varphi (x,t)⟩}{\partial t}=\Delta ⟨\varphi (x,t)⟩-⟨\varphi (x,t)⟩-⟨{\varphi }^3(x,t)⟩.$$

$$H[\varphi ;\alpha ]=H[\varphi ]+\alpha \int \varphi (x)dx$$

$$\alpha ={\mathrm{argmax}}_{\alpha }[\alpha \mu -F(\alpha )],$$

$$(\partial /\partial \alpha )[\alpha \mu -F(\alpha )]=\mu -\mu (\alpha ),$$

$$\varphi (x,t+\delta t)=\varphi (x,t)-\delta t[\varphi (x,t)+{\varphi }^3(x,t)]+$$

$$\frac{D\delta t}{(\delta x{)}^2}[\varphi (x+\delta x,t)-2\varphi (x,t)+\varphi (x-\delta x,t)]+$$

$$\sqrt{2kT(\delta t/\delta x)}N(x,t)$$

$$H_{\delta }=\frac{D}{2\delta x}\sum _{⟨x,x^{\prime }⟩}(\varphi (x)-\varphi (x^{\prime }){)}^2$$

$$+\sum _x\delta x[\frac{1}{2}{\varphi }^2(x)+\frac{1}{4}{\varphi }^4(x)]$$

$$H_{\delta }[\varphi ;\alpha ]=H_{\delta }[\varphi ]+\alpha \sum _x\delta x\varphi (x).$$

$${\dot{\mu }}_{app}(t)=[\mu (t+\delta t)-\mu (t)]/\delta t.$$

$$\mu (t+\Delta t)=\mu (t)+\Delta t{\dot{\mu }}_{app}(t).$$

Full text

I Introduction

Accurate, fast simulations of complex, large-scale, nonlinear systems remain a challenge for computational science and engineering, despite extraordinary advances in computing power. Examples range from molecular dynamics simulations of proteins MDProteinsI , MDProteinsII and glasses MDGlasses , to stochastic simulations of cellular biochemistry  Gillespie ; Wilkinson , to global-scale, geophysical fluid dynamics MajdaWang . Often for the systems under consideration there is no obvious scale separation, and their many degrees of freedom are strongly coupled. The complex and multiscale nature of these processes therefore makes them extremely difficult to model numerically. To make matters worse, one is often interested not in a single, time-dependent solution of the equations governing these processes, but rather in ensembles of solutions consisting of multiple realizations (e.g., sampling noise, initial conditions, and/or uncertain parameters). Often real-time answers are needed (e.g., for control, tracking, filtering). These demands can easily exceed the computational resources available not only now but also for the foreseeable future.

In principle, all statistical information for the problem under investigation is contained in solutions to the Liouville (if deterministic)/Kolmogorov (if stochastic) equations. These are partial differential equations in a state space of high (possibly infinite) dimension. A straightforward discretization of the Liouville / Kolmogorov equations is therefore impractical. An ensemble approach to solving these equations can be taken; however, quite often, the practical application of the ensemble approach is also problematic. Generating a sufficient number of independent samples for statistical convergence can be a challenge. For some problems, computing even one realization may be prohibitive.

The traditional approach to making these problems computationally tractable is to replace the Liouville/Kolmogorov equation by a (small) set of equations (PDEs or ODEs) for a few, low order statistical moments of its solution. When taking this approach for nonlinear systems, one must make an approximation, a closure, for the dependence of higher order moments on lower order moments. Typically the form of the closure equation is based on expert knowledge, empirical data, and/or physical insight. For example, in the superposition approximation and its extensions Superposition for dense liquids and plasmas, both quantum or classical, one approximates third order moments as functions of second order moments. Moment closure methods of this type have been applied to a number of areas including fluid turbulence (see  TurbClosure and references therein, and also the work of Chorin et. al.). Of course, as with any approximation strategy, the quality of the resulting reduced description depends on the approximations made – poor closures lead to poor answers/predictions. In addition to replacing the ensemble with a small set of equations for low order moments, these equations are typically easier to solve. They are deterministic and generally far less stiff than the original equations.

A less exploited variant of this approximation scheme is the probability density function (PDF) based moment-closure approach. For PDF moment closures one makes an ansatz for the system statistics guided by available information (e.g., symmetries). One then uses this ansatz in conjunction with the original dynamical equations to derive moment equations. Such PDF-based closures have been developed for reacting scalars advected by turbulence CCK , phase-ordering dynamics Oono and a variety of other systems. This approach to moment-closure is a close analogue of the Rayleigh-Ritz method frequently used in solving the quantum-mechanical Schroedinger equation, by exploiting an ansatz for the wave-function. For a formal development of this point of view, see Eyink96 .

One of the obstacles to applying moment closures is that often the closure equations are too complicated to write down explicitly, even with the availability of computer algebra / symbolic computation systems. This is especially true for large-scale, complex systems, e.g. global climate models. Because of their great complexity, even if one could in principle derive the closure equations analytically, this procedure would be extremely difficult and time-intensive. Moreover, each time a model is updated, as climate and ocean models regularly are, the closure equations would have to be rederived. In other cases it may simply be impossible to determine the closure equations analytically. This is especially likely when PDF’s are not Gaussian, which is the case for most useful closures. Monte Carlo or other numerical methods may be needed in order to evaluate integrals for the moments Levermore96 . In addition, there may be situations where neither analytic nor numerical/MC integration will yield the closure equations due to the black-box nature of the available numerical simulator such as a compiled numerical code with an inaccessible source. Clearly, a need exists for a robust approach to the general closure protocol which circumvents analytical difficulties.

We address that need here by combining PDF closures with equation free modeling EQNFREE  EqnFree2 . The basic premise of the equation-free method is to use an ensemble of short bursts of simulation of the original dynamical system to estimate, on demand, the time-evolution of the the closure equations that we may not explicitly have. The equation-free approach extends the applicability of statistical closures beyond the rare cases where they can be expressed in closed form. This hybrid strategy may be faster than the brute-force solution of a large ensemble of realizations of the dynamical equations since the closure version is generally smoother than the original problem.

This paper is organized as follows. In Section 2 we describe the general features of PDF-based moment closures. In Section 3 we explain how to implement the equation-free approach with these closures. We then, in Section 4, apply these ideas for a specific dynamical system, the stochastic Ginzburg-Landau (GL) equations using a particular PDF-based closure scheme, the entropy method of Eyink and Levermore EL . We conclude with a discussion of closure quality, computational issues, and the application of our approach to large-scale systems.

II PDF-Based Moment Closures

We consider the very general class of dynamical systems, including maps, formally represented by

[TABLE]

or

[TABLE]

where ${\bf N}(t)$ is a stochastic process with prescribed statistics. The stochastic component arises from unknown parameters, random forcing, neglected degrees of freedom and/or random initial conditions. This class includes both deterministic and stochastic systems with discrete and/or continuous states. Queueing systems, molecular dynamics, and stochastic PDEs are just some of the many examples that fall into this category.

For concreteness in this paper we restrict ourselves to a special case of equation (2), namely, situations where ${\bf N}(t)$ is a Markov process (Brownian motion, Poisson process, etc.) and—more specifically still—Itô stochastic differential equations of the form:

[TABLE]

The deterministic component of the state, ${\bf X}$, is governed by the continuously differentiable vector field, ${\bf U}:{\mathbb{R}}^{N}\times{\mathbb{R}}\to{\mathbb{R}}^{N}$. For many problems of interest (e.g., climate) ${\bf U}$ is a highly nonlinear function. The noise component is modeled by the standard mean [math], covariance matrix ${\bf I}$ Wiener process, ${\bf W}\in{\mathbb{R}}^{N}$, possibly modulated by a state-dependent matrix ${\bf S}:{\mathbb{R}}^{N}\times{\mathbb{R}}\to{\mathbb{R}}^{N\times N}$. Equation (3) encompasses a wide class of systems including deterministic $({\bf S}=0)$ ones.

In many cases one is interested in knowing the low order statistics of equation (3), for example an instantaneous mean value or possibly multi-point covariance of ${\bf X}$. These statistics can be obtained by averaging over an ensemble of stochastic systems, solving equation (3). They can also be obtained via the forward Kolmogorov equation for the probability density function $P({\bf X},t)$:

[TABLE]

where $P$ satisfies the conditions: $P({\bf X},t)\geq 0$, and $\int P({\bf X},t)\,d{\bf X}=1$, and where ${\cal L}^{*}$ is the generator of the Markov process. In the case of equation (3) this operator takes the form

[TABLE]

The forward Kolmogorov equation then becomes a Fokker-Planck equation

[TABLE]

where ${\bf D}({\bf X},t)={\bf S}({\bf X},t){\bf S}({\bf X},t)^{T}$ is the nonnegative-definite diffusion matrix arising from the noise term. Unlike the original dynamical equation (3), the forward Kolmorogov equation (FKE) is both linear and deterministic. Dealing with it, therefore, has apparent advantages over the original ensemble of stochastic systems simulations. The price to pay for these advantages is that the FKE lives in a typically high, potentially infinite-dimensional, space. When equation (3) is a nonlinear PDE, numerical solution to the FKE is usually ruled out.

For computational purposes, we would therefore like to reduce the FKE (if possible and useful) to a small system of ordinary differential equations. This reduction should simplify the computation as much as possible while retaining fidelity to the original dynamical processes. The reduction proceeds by taking moments of the FKE with respect to a vector-valued function {\mbox{\boldmath\xi}}({\bf X},t) from ${\mathbb{R}}^{N}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{M}$. The $\xi$ selected should include the relevant variables in the system (slow modes, conserved quantities, etc.). The moments {\mbox{\boldmath\mu}}(t) of {\mbox{\boldmath\xi}}({\bf X},t) are defined by

[TABLE]

and give rise to

[TABLE]

where

[TABLE]

and ${\cal L}$ is the adjoint of ${\cal L}^{*}$ or the backward Kolmogorov operator. The result (8) can be obtained by averaging over an ensemble of realizations of the stochastic dynamics (3). In general, however, (8) is not a closed equation for the moments, $\mu$. One can close this equation by choosing a PDF, P({\bf X},t,{{\mbox{\boldmath\mu}}}), which itself is a function of the moments $\mu$.

[TABLE]

Alternatively, one can select a family of probability densities P({\bf X},t,{\mbox{\boldmath\alpha}}), specified by parameters {\mbox{\boldmath\alpha}}={\mbox{\boldmath\alpha}}({\mbox{\boldmath\mu}},t) rather than directly by the moments $\mu$. This is analogous to specifying the temperature in the canonical ensemble as opposed to the average energy. The equivalence of these approaches is guaranteed provided that the parameters and moments can be determined uniquely from one another. The translation between the parameters and their corresponding moments can be carried out by one of several methods. In some cases one may require Monte Carlo evaluation of the resulting integrals.

If the moments and/or parameters are selected judiciously, one hopes that the approximate PDF P({\bf X},t,{\mbox{\boldmath\alpha}}({\mbox{\boldmath\mu}})(t)) will be close to the exact solution of the Liouville/Kolmogorov equation (4). The mapping closure approach of Chen et al CCK and the Gaussian mapping method of Yeung et al. Oono are based on this type of parametric PDF closure 111In the case of  CCK the dynamics is an advection-reaction-diffusion equation for a scalar concentration field ${\bf X}(t)=\{\theta({\bf x},t):\,{\bf x}\in{\mathbb{R}}^{d}\}$. The moment functions are the “fine-grained PDF” $\xi_{\vartheta,{\bf x}}[{\bf X},t]=\delta(\theta({\bf x},t)-\vartheta),$ labelled by space point ${\bf x}$ and scalar value $\vartheta.$ The moment average $\mu_{\vartheta,{\bf x}}(t)=\langle\delta(\theta({\bf x},t)-\vartheta)\rangle$ is the 1-point PDF $p(\vartheta;{\bf x},t)$ which gives the distribution of scalar values $\vartheta$ at space-time point $({\bf x},t).$ The parametric model P[{\bf X};{\mbox{\boldmath\alpha}},t] is the distribution over scalar fields obtained by the ansatz $\theta({\bf x},t)=X(\theta_{0}({\bf x},t),{\bf x},t)$ where $\theta_{0}({\bf x},t)$ is a reference random field of known (Gaussian) statistics and $X(\cdot,{\bf x},t):{\mathbb{R}}\rightarrow{\mathbb{R}}$ is a “mapping function”. The latter function is the “parameter” $\alpha_{\vartheta_{0},{\bf x}}(t)=X(\vartheta_{0},{\bf x},t)$ which determines (and is determined by) the “moment” $\mu_{\vartheta,{\bf x}}(t)$ from the relation $p(X(\vartheta_{0},{\bf x},t);{\bf x},t)|\partial X/\partial\vartheta_{0}|=p_{0}(\vartheta_{0},{\bf x},t).$ Here $p_{0}$ is the 1-point PDF of the reference Gaussian field $\theta_{0}({\bf x},t).$

       The approach of Oono is similar. The problem is phase-ordering dynamics as given, for example, by our equation (12) and ${\bf X}(t)=\{\phi({\bf x},t):\,{\bf x}\in{\mathbb{R}}^{d}\}.$ The moment functions are the quadratic products \xi_{\bf r}[{\bf X},t]=\phi({\bf r},t)\phi({\mbox{\boldmath0}},t), labelled by the displacement ${\bf r}\in{\mathbb{R}}^{d}$ and the moment averages $\mu_{\bf r}(t)$ are the spatial correlation function $C({\bf r},t)$. The parametric model P[{\bf X};{\mbox{\boldmath\alpha}},t] is the distribution obtained by the ansatz $\phi({\bf x},t)=f(u({\bf x},t))$ where $u({\bf x},t)$ is a homogeneous Gaussian random field with mean zero and covariance G({\bf r},t)=\langle u({\bf r},t)u({\mbox{\boldmath0}},t)\rangle and $f(z)$ is the stationary planar interface solution of the TDGL equation (12). In this case, it is the auxiliary correlation function $G({\bf r},t)$ which plays the role of the “parameter” $\alpha_{\bf r}(t).$ It is shown in Oono for various cases how this function may be uniquely related to the “moment” $\mu_{\bf r}(t)=C({\bf r},t).$. In fact, perhaps the most familiar application of the parametric approach is the use of the Rayleigh-Ritz method in quantum mechanical calculations. This is the essential approach of our paper.

III Equation-Free Computation

Although we now have obtained a closed moment equation (equation  10), we still need to determine the dynamical vector field ${\bf V}$. As explained above, this step can be a serious obstacle to the practical implementation of PDF-based moment-closure (PDFMC). A method to calculate ${\bf V}$ is desirable that (i) does not require a radical revision each time the underlying code or model changes, and (ii) is relatively insensitive to the complexity of the PDFMC. The equation-free approach of Kevrekidis and collaborators EQNFREE meets those requirements. It permits one to work with much more sophisticated, physically realistic closures.

Equation-free computation is motivated by the simple observation that numerical computations involving the closure equations ultimately do not require closed formulae for the closure equations. Instead, one must only be able to sample an ensemble of system states ${\bf X}$ distributed according to the closure ansatz P({\bf X},t;{\mbox{\boldmath\alpha}}) and then evolve each of these via equation (3) for short intervals of time. Such sampling and subsequent dynamical evolution would be necessary to calculate the statistics of interest even when not using a closure strategy. It is sufficient to have a (possibly black-box) subroutine available which, given a specific state variable ${\bf X}(t)$ as input, returns the value of the state ${\bf X}(t+\delta t)$ after a short time $\delta t$. The ensemble of systems, each of which satisfies equation (3), is evolved over a time interval $\delta t$. The moments/parameters $\mu$ or $\alpha$ are determined at the beginning and end of this interval and the time derivative {\dot{{\mbox{\boldmath\mu}}}} is estimated from the results of these short ensemble runs. This “coarse timestepper” can be used to estimate locally the right hand side of the closure evolution equations, namely {\bf V}({{\mbox{\boldmath\mu}}},t).

Coarse projective forward Euler (arguably the simplest of equation-free algorithms) which we will use below illustrates the approach succinctly: Starting from a set of coarse-grained initial conditions specified by moments {\mbox{\boldmath\mu}}(t) we first (a) lift to a consistent fine scale description, that is, sample the PDF ansatz P({\bf X},t;{\mbox{\boldmath\alpha}}(t)) to generate ensembles of initial conditions ${\bf X}$ for equation (3) consistent with the set {\mbox{\boldmath\mu}}(t); (b) starting with these consistent initial conditions we evolve the fine scale description for a (relatively short) time $\delta t$; we subsequently restrict back to coarse observables by evaluating the moments {\mbox{\boldmath\mu}}(t+\delta t) as ensemble-averages and (d) use the results to estimate locally the time derivative d{\mbox{\boldmath\mu}}/dt. This is precisely the right hand-side of the explicitly unavailable closure, obtained not through a closed form formula, but rather through short, judicious computational experiments with the original fine scale dynamics/code. Given this local estimate of the coarse-grained observable time derivatives, we can now exploit the smoothness of their evolution in time (in the form of Taylor series) and take a single long projective forward Euler step:

[TABLE]

The procedure then repeats itself: lifting, fine scale evolution, restriction, estimation, and then (connecting with continuum traditional numerical analysis) a new forward Euler step. Beyond coarse projective forward Euler, many other coarse initial-value solvers (e.g. coarse projective Adams-Bashforth, and even implicit coarse solvers) have been implemented; the stability and accuracy study of such algorithms is progressing EQNFREE . These developments allow us to construct a nonintrusive implementation of PDF moment closures, nonintrusive in the sense that we compute with the closures without explicitly obtaining them, but rather by intelligently chosen computational experiments with the original, fine-scale problem.

There is, however, an obvious objection to the equation-free implementation of moment-closures. Using the same ingredients, one can clearly obtain an estimate of any statistics of interest (for example, the moment-averages {\mbox{\boldmath\mu}}(t)) without the need of making any closure assumptions whatsoever. This can be done by the much simpler method of direct ensemble averaging. That is, one can sample an ensemble of initial conditions ${\bf X}$ from any chosen distribution $P_{0}({\bf X})$, evolve each of these realizations according to the fine-scale dynamics of equation (3), and then evaluate any statistics of interest at time $t$ by averaging over the ensemble of solutions ${\bf X}(t).$ It would seem that this direct ensemble approach is much more straightforward and accurate than the equation-free implementation of a moment-closure, which introduces additional statistical hypotheses.

The response to this important objection is that the fine-scale dynamics (3) is often very stiff for the applications considered, in which the system contains many-degrees-of-freedom interacting on a huge range of length- and time-scales. In contrast, the closure equation (10) is much less stiff, because of statistical-averaging, and its solutions {\mbox{\boldmath\mu}}(t) are much smoother in time (and space). Thus, to evolve an ensemble of solutions of the fine-scale dynamics (3) from an initial time $t_{0}$ to a final time $t_{0}+T$ would require $O(T/\delta t)$ integration steps, where the time-step $\delta t$ is required to be very small by the intrinsic stiffness of the micro-dynamics. In the closure approach, the evolution of the moment equations (10) from time $t_{0}$ to time $t_{0}+T$ requires only $O(T/\Delta t)$ integration steps, with (hopefully) $\Delta t\gg\delta t.$ Each of these closure integration steps by an increment $\Delta t$ requires in the equation-free approach just one (or just a few) fine-scale integration step by an increment $\delta t.$ Thus, there is an over-all savings by a (hopefully) large factor $O(\Delta t/\delta t).$ This crude estimate is based on a single step coarse projective forward Euler algorithm; clearly, more sophisticated projective integration algorithms can be used.

In all of them, however, the computational savings are predicated on the smoothness of the closure equations, and are governed by the ratio of the time that it takes to obtain a good local estimate of d{\mbox{\boldmath\mu}}/dt from full direct simulation to the time that we can (linearly or even polynomially) extrapolate {\mbox{\boldmath\mu}}(t) in time. It is also worth noting that a variety of additional computational tasks, beyond projective integration (e.g. accelerated fixed point computation) can be performed within the equation-free framework

In the next section we show by a concrete example how significant computational economy can be achieved with statistical moment closures implemented in the equation-free framework.

IV A Numerical Example

We illustrate here the equation-free implementation of moment-closures for a canonical equation of phase-ordering kinetics Bray , the stochastic time-dependent Ginzburg-Landau (TDGL) equation in one spatial dimension. This is written as

[TABLE]

where $\phi(x,t)$ represents a local order parameter, e.g. a magnetization. The noise has mean zero and covariance $\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle=2kT\delta(x-x^{\prime})\delta(t-t^{\prime})$. The potential $V$ shall be chosen as

[TABLE]

to represent a single quartic/quadratic well. This stochastic dynamics has an invariant measure which is formally of Hamiltonian form $P_{*}[\phi]\propto\exp(-H[\phi]/kT)$ where

[TABLE]

The Gibbsian measure $P_{*}[\phi]$ is approached at long times for any random distribution $P_{0}[\phi]$ of initial states.

One of the simplest dynamical quantities of interest is the bulk magnetization $\overline{\phi}(t)=(1/V)\int\phi(x,t)dx,$ where $V$ is the total volume. If the initial statistics are space-homogeneous, then the ensemble average $\mu(t)=\langle\overline{\phi}(t)\rangle$ is also given by $\mu(t)=\langle\phi(x,t)\rangle$ for any space point $x.$ Equation (12) leads to a hierarchy of equations for statistical moments of $\phi(x,t)$. For example, the first moment satisfies the equation

[TABLE]

The evolution of the mean total magnetization is thus a function of the mean cubic total magnetization. One could write a time evolution equation for $\langle\phi^{3}\rangle$, but it would involve a higher order term $\langle\phi^{5}\rangle$, and so on. Each equation contains higher moments and therefore the hierarchy does not close.

To close the equation for $\mu(t)$ we assume a parametric PDF of the form $P[\phi;\alpha]\propto\exp(-H[\phi;\alpha]/kT)$ where

[TABLE]

is a perturbation of the Hamiltonian (13) by a term proportional to the moment variable $\xi[\phi]=(1/V)\int\phi(x)\,dx.$ This is a special case of a general “entropy-based” closure prescription proposed by Eyink and Levermore EL . This closure scheme guarantees that $\alpha(t)\to 0$ at long times and therefore the PDF ansatz $P[\phi;\alpha(t)]$ relaxes to the correct stationary distribution $P_{*}[\phi]$ of the stochastic process. The determination of the parameter $\alpha$ given the moment $\mu$ is here accomplished by Legendre transform

[TABLE]

where the “moment-generating function” $F(\alpha)=\log\langle\exp[\alpha\int\phi(x)\,dx]\rangle_{*}$ and $\langle\cdot\rangle_{*}$ denotes average with respect to the invariant measure $P_{*}[\phi]$. The numerical optimization required for the Legendre transform is well-suited to gradient descent algorithms such as the conjugate gradient method, since

[TABLE]

where $\mu(\alpha)=\langle\xi\rangle_{\alpha}$ is the average of the moment-function in the PDF ansatz $P[\phi;\alpha].$ In simple cases, $F(\alpha)$ and $\mu(\alpha)=F^{\prime}(\alpha)$ may be given by closed analytical expressions. If not, then both of these averages may be determined together by Monte Carlo sampling techniques.

In the numerical calculations below, we discretize equation (12) using a forward Euler-Maruyama stochastic integrator and 3-point stencil for the Laplacian (other discretizations are possible).

[TABLE]

where $N(x,t)$ are independent, identically distributed standard normal random variables for each space-time point $(x,t).$ The invariant distribution of the stochastic dynamics space-discretized in this manner has a Gibbsian form $\propto\exp(-H_{\delta}/kT)$ with discrete Hamiltonian

[TABLE]

where $\langle x,x^{\prime}\rangle$ are nearest-neighbor pairs. The closure ansatz can be adopted in the consistently discretized form $P_{\delta}[\phi;\alpha]\propto\exp(-H_{\delta}[\phi;\alpha]/kT)$ where

[TABLE]

In this numerical experiment, we integrate an $N=1000$ member ensemble of solutions of equation (IV), and measure the ensemble-averaged, global magnetization $\mu(t)=\langle\overline{\phi}(t)\rangle=(1/V)\sum_{x}\langle\phi(x,t)\rangle$ at each time-step. With this we compare the results of the entropy-based closure simulation implemented by the equation-free framework using also an ensemble with $N=1000$ samples. In this concrete example, the projective integration scheme works as follows: Suppose we are given the parameter $\alpha(t)$ at time $t$. The mean $\mu(t)$ is first calculated from the parametric ensemble at time $t$ by Monte Carlo sampling. Next all $N$ samples are integrated over a short time-step $\delta t$ to create a time-advanced ensemble. From this ensemble $\mu(t+\delta t)$ is calculated, which yields an estimate of the local time derivative.

[TABLE]

A large, projective Euler time-step of the moment average is then taken via

[TABLE]

The parameter is finally updated by using the Legendre transform inversion to obtain $\alpha(t+\Delta t)$ from the known value $\mu(t+\Delta t).$ The cycle may now be repeated to integrate the closure equations by successive time-steps of length $\Delta t.$

A critical issue in general application of projective integration is the criterion to determine the projective time-step $\Delta t.$ For stiff problems with time-scale separation, the projective time step for stability purposes is of the order of (1/fastest “slow group” eigenvalues), while the “preparatory” simulation time is of the order of (1/slowest “fast group” eigenvalue). Variants of the approach have been developed for problems with several gaps in their spectrum GK1 . Accuracy considerations in real-time projective step selection can, in principle, be dealt with in the traditional way for integrators with adaptive step-size selection and error control: through on-line a posteriori error estimates. An additional “twist” arises from the error inherent in the estimation of the (unavailable) reduced time derivatives from the ensemble simulations; issues of variance reduction and even on-line hypothesis testing (are the data consistent with a local linear model?) must be considered. These are important research issues that are currently explored by several research groups. Nevertheless, the main factor in computational savings comes from the effective smoothness of the unavailable closed equation: the separation of time scales between the low-order statistics we follow and the higher order statistics whose effect we model (and, eventually, the time scales of the direct simulation of the original model).

Figure 1 is a plot comparing Projective Integration with Entropy Closure and direct Ensemble Integration with equation (12) for diffusion constant $D=1000.0$ We have selected both the “fine-scale” integration step $\delta t$ and the “coarse-scale” projective integration step $\Delta t$ to be as large as possible, consistent with stability and accuracy. Thus, only steps small enough to avoid numerical blow-ups were considered. Then, values were selected both for $\delta t$ and for $\Delta t$ so that the numerical integrations with those time-steps differed by at most a few percent from fully converged integrations with very small steps. In this manner, the time step required for the Euler-Maruyama integration of (12) was determined to be $\delta t=0.0004$. On the other hand, for projective integration of the closure equation a time step $\Delta t=0.01$ could be taken. This indicates a gain in time step by a factor of 25, which is also roughly the speed-up in the algorithm or savings in CPU time. The present example is not as stiff as equations that appear in more realistic applications, with a very broad range of length- and time-scales, where even greater computational economies might be expected.

In general, the moment-closure results need not agree so well with those of the direct ensemble approach, even when both are converged. In the example presented here, there is good agreement because the closure effectively captures the one-point PDF (see Fig.2). This one-point PDF is the only statistical quantity that enters into Equation (14) as long as the statistics are homogeneous and the Laplacian term vanishes.

V Conclusions

In this paper, we have described how one can combine recently developed equation-free methods with statistical moment closures to model nonlinear problems. With this method we can numerically integrate complex nonlinear systems, for which closure equations may not be available in closed form. In the example presented here the specific entropy-based closure we selected has an H-theorem which guarantees relaxation to the equilibrium state of the original dissipative dynamics. However, we stress that the general approach outlined above can be used with a variety of closure methods.

The equation-free method has the potential to enhance the flexibility, power, and applications set of the statistical moment closure approach. Since little or no analytic work is required, the sophistication of statistical moment closures can greatly enhanced beyond Gaussian PDF ansätze. The “practical usefulness” criterion for parametric PDF models that they permit analytical calculations is replaced by the criterion that they can be efficiently sampled. We believe that this approach can significantly increase the usefulness of closure methods.

In order to model systems like global climate, oceans, and reaction diffusion processes in systems biology, one will have to construct more complex closures. These will likely include higher order moments, correlation functions of the relevant variables, highly non-Gaussian statistics, etc. As the closures become more complex, the lifting step will require more efficient sampling approaches. One will likely have to use nonlocal, accelerated sampling methods. One will also likely employ the latest in adaptive time and adaptive mesh methods to optimize performance for large-scale problems.

VI Acknowledgements

This work, LA-UR-07-2218, was carried out in part at Los Alamos National Laboratory under the auspices of the US National Nuclear Security Administration of the US Department of Energy. It was supported under contract number DE-AC52-06NA25396. The work of IGK was partially supported by DARPA and by and US DOE(CMPD). G. Eyink was supported by NSF-ITR grant, DMS-0113649.