Equation-Free Implementation of Statistical Moment Closures
Francis J. Alexander and Gregory Johnson
Los Alamos National Laboratory, P.O.Box 1663,
Los Alamos, NM, 87545.
Gregory L. Eyink
Department of Mathematical Sciences
Johns Hopkins University
Baltimore, MD 21218
Ioannis G. Kevrekidis
Department of Chemical Engineering and PACM
Princeton University
Princeton, NJ 08544
We present a general numerical scheme for the practical implementation of statistical moment clo-
sures 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 typi-
cally 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.
PACS numbers:
INTRODUCTION
Accurate, fast simulations of complex, large-scale, non-
linear systems remain a challenge for computational sci-
ence and engineering, despite extraordinary advances in
computing power. Examples range from molecular dy-
namics simulations of proteins [1], [2] and glasses [3], to
stochastic simulations of cellular biochemistry [4, 5], to
global-scale, geophysical fluid dynamics [6]. Often for the
systems under consideration there is no obvious scale sep-
aration, and their many degrees of freedom are strongly
coupled. The complex and multiscale nature of these pro-
cesses therefore makes them extremely difficult to model
numerically. To make matters worse, one is often in-
terested not in a single, time-dependent solution of the
equations governing these processes, but rather in ensem-
bles of solutions consisting of multiple realizations (e.g.,
sampling noise, initial conditions, and/or uncertain pa-
rameters). 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 Liou-
ville (if deterministic)/Kolmogorov (if stochastic) equa-
tions. These are partial differential equations in a state
space of high (possibly infinite) dimension. A straightfor-
ward discretization of the Liouville / Kolmogorov equa-
tions is therefore impractical. An ensemble approach to
solving these equations can be taken; however, quite of-
ten, 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 real-
ization may be prohibitive.
The traditional approach to making these prob-
lems computationally tractable is to replace the Liou-
ville/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 phys-
ical insight. For example, in the superposition approxi-
mation and its extensions [7] for dense liquids and plas-
mas, 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 [8]
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 ensem-
ble 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
http://arxiv.org/abs/0704.0804v1
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 in-
formation (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 turbu-
lence [10], phase-ordering dynamics [11] 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 [12].
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 com-
puter 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 up-
dated, 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 clo-
sure equations analytically. This is especially likely when
PDF’s are not Gaussian, which is the case for most use-
ful closures. Monte Carlo or other numerical methods
may be needed in order to evaluate integrals for the mo-
ments [13]. 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 numer-
ical code with an inaccessible source. Clearly, a need ex-
ists 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 [14] [15]. 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 clo-
sure equations that we may not explicitly have. The
equation-free approach extends the applicability of sta-
tistical 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 clo-
sures. 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 sys-
tem, the stochastic Ginzburg-Landau (GL) equations us-
ing a particular PDF-based closure scheme, the entropy
method of Eyink and Levermore [16]. We conclude with
a discussion of closure quality, computational issues, and
the application of our approach to large-scale systems.
PDF-BASED MOMENT CLOSURES
We consider the very general class of dynamical sys-
tems, including maps, formally represented by
Ẋ = U(X(t),N(t), t) (1)
Xt+1 = Ut(Xt,Nt) (2)
where N(t) is a stochastic process with prescribed statis-
tics. The stochastic component arises from unknown pa-
rameters, 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
N(t) is a Markov process (Brownian motion, Poisson pro-
cess, etc.) and—more specifically still—Itô stochastic dif-
ferential equations of the form:
dX = U(X, t)dt +
2S(X, t)dW(t). (3)
The deterministic component of the state, X, is gov-
erned by the continuously differentiable vector field, U :
N × R → RN . For many problems of interest (e.g., cli-
mate) U is a highly nonlinear function. The noise com-
ponent is modeled by the standard mean 0, covariance
matrix I Wiener process, W ∈ RN , possibly modulated
by a state-dependent matrix S : RN×R → RN×N . Equa-
tion (3) encompasses a wide class of systems including
deterministic (S = 0) ones.
In many cases one is interested in knowing the low
order statistics of equation (3), for example an instanta-
neous mean value or possibly multi-point covariance of
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 (X, t):
∂tP = L∗(t)P, (4)
where P satisfies the conditions: P (X, t) ≥ 0, and
P (X, t) dX = 1, and where L∗ is the generator of the
Markov process. In the case of equation (3) this operator
takes the form
L∗(t)ψ(X) = −∇X·(U(X, t)ψ(X))+∇2X : (D(X, t)ψ(X)).
The forward Kolmogorov equation then becomes a
Fokker-Planck equation
∂tP +∇X · (UP ) = ∇2X : (DP ) (6)
where D(X, t) = S(X, t)S(X, t)T is the nonnegative-
definite diffusion matrix arising from the noise term. Un-
like the original dynamical equation (3), the forward Kol-
morogov equation (FKE) is both linear and determinis-
tic. Dealing with it, therefore, has apparent advantages
over the original ensemble of stochastic systems simu-
lations. 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 reduc-
tion should simplify the computation as much as pos-
sible while retaining fidelity to the original dynamical
processes. The reduction proceeds by taking moments of
the FKE with respect to a vector-valued function ξ(X, t)
from RN × R+ → RM . The ξ selected should include
the relevant variables in the system (slow modes, con-
served quantities, etc.). The moments µ(t) of ξ(X, t) are
defined by
µ(t) =
ξ(X, t)P (X, t)dX (7)
and give rise to
µ̇(t) =
ξ̇(X, t)P (X, t)dX, (8)
where
ξ̇(X, t) = ∂tξ(X, t) + L(t)ξ(X, t) (9)
and L is the adjoint of L∗ or the backward Kolmogorov
operator. The result (8) can be obtained by averaging
over an ensemble of realizations of the stochastic dynam-
ics (3). In general, however, (8) is not a closed equa-
tion for the moments, µ. One can close this equation by
choosing a PDF, P (X, t,µ), which itself is a function of
the moments µ.
µ̇(t) = V(µ, t) ≡
ξ̇(X, t)P (X, t,µ)dX. (10)
Alternatively, one can select a family of probability den-
sities P (X, t,α), specified by parameters α = α(µ, t)
rather than directly by the moments µ. 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 param-
eters 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 (X, t,α(µ)(t)) will be close to the exact solution of
the Liouville/Kolmogorov equation (4). The mapping
closure approach of Chen et al [10] and the Gaussian
mapping method of Yeung et al. [11] are based on this
type of parametric PDF closure [19]. In fact, perhaps
the most familiar application of the parametric approach
is the use of the Rayleigh-Ritz method in quantum me-
chanical calculations. This is the essential approach of
our paper.
EQUATION-FREE COMPUTATION
Although we now have obtained a closed moment equa-
tion (equation 10), we still need to determine the dynam-
ical vector field 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
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 collabora-
tors [14] 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 for-
mulae for the closure equations. Instead, one must only
be able to sample an ensemble of system states X dis-
tributed according to the closure ansatz P (X, t;α) and
then evolve each of these via equation (3) for short inter-
vals 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 avail-
able which, given a specific state variable X(t) as input,
returns the value of the state X(t + δt) after a short
time δt. The ensemble of systems, each of which satis-
fies equation (3), is evolved over a time interval δt. The
moments/parameters µ or α are determined at the be-
ginning and end of this interval and the time derivative
µ̇ 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 equa-
tions, namely V(µ, t).
Coarse projective forward Euler (arguably the simplest
of equation-free algorithms) which we will use below il-
lustrates the approach succinctly: Starting from a set
of coarse-grained initial conditions specified by moments
µ(t) we first (a) lift to a consistent fine scale descrip-
tion, that is, sample the PDF ansatz P (X, t;α(t)) to
generate ensembles of initial conditions X for equation
(3) consistent with the set µ(t); (b) starting with these
consistent initial conditions we evolve the fine scale de-
scription for a (relatively short) time δt; we subsequently
restrict back to coarse observables by evaluating the mo-
ments µ(t+ δt) as ensemble-averages and (d) use the re-
sults to estimate locally the time derivative dµ/dt. This
is precisely the right hand-side of the explicitly unavail-
able closure, obtained not through a closed form formula,
but rather through short, judicious computational exper-
iments 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:
µ(t+∆t) = µ(t) + ∆t
µ(t+ δt)− µ(t)
. (11)
The procedure then repeats itself: lifting, fine scale evo-
lution, restriction, estimation, and then (connecting with
continuum traditional numerical analysis) a new for-
ward Euler step. Beyond coarse projective forward Eu-
ler, many other coarse initial-value solvers (e.g. coarse
projective Adams-Bashforth, and even implicit coarse
solvers) have been implemented; the stability and accu-
racy study of such algorithms is progressing [14]. These
developments allow us to construct a nonintrusive imple-
mentation of PDF moment closures, nonintrusive in the
sense that we compute with the closures without explic-
itly 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 µ(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 X from
any chosen distribution P0(X), evolve each of these real-
izations 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 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 appli-
cations 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 µ(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 t0 to a final time t0+T
would require O(T/δt) integration steps, where the time-
step δt is required to be very small by the intrinsic stiff-
ness of the micro-dynamics. In the closure approach, the
evolution of the moment equations (10) from time t0 to
time t0+T requires only O(T/∆t) integration steps, with
(hopefully) ∆t ≫ δt. Each of these closure integration
steps by an increment ∆t requires in the equation-free
approach just one (or just a few) fine-scale integration
step by an increment δt. Thus, there is an over-all savings
by a (hopefully) large factor O(∆t/δt). This crude esti-
mate 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µ/dt from full direct
simulation to the time that we can (linearly or even poly-
nomially) extrapolate µ(t) in time. It is also worth not-
ing that a variety of additional computational tasks, be-
yond 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.
A NUMERICAL EXAMPLE
We illustrate here the equation-free implementation
of moment-closures for a canonical equation of phase-
ordering kinetics [17], the stochastic time-dependent
Ginzburg-Landau (TDGL) equation in one spatial di-
mension. This is written as
∂φ(x, t)
= D∆φ(x, t) − V ′(φ(x, t)) + η(x, t) (12)
where φ(x, t) represents a local order parameter, e.g. a
magnetization. The noise has mean zero and covariance
〈η(x, t)η(x′, t′)〉 = 2kT δ(x−x′)δ(t− t′). The potential V
shall be chosen as
V (φ) =
to represent a single quartic/quadratic well. This
stochastic dynamics has an invariant measure which is
formally of Hamiltonian form P∗[φ] ∝ exp(−H [φ]/kT )
where
H [φ] =
D|∇φ(x)|2 + V (φ(x))] dx. (13)
The Gibbsian measure P∗[φ] is approached at long times
for any random distribution P0[φ] of initial states.
One of the simplest dynamical quantities of interest is
the bulk magnetization φ(t) = (1/V )
φ(x, t)dx, where
V is the total volume. If the initial statistics are space-
homogeneous, then the ensemble average µ(t) = 〈φ(t)〉 is
also given by µ(t) = 〈φ(x, t)〉 for any space point x. Equa-
tion (12) leads to a hierarchy of equations for statistical
moments of φ(x, t). For example, the first moment satis-
fies the equation
∂〈φ(x, t)〉
= ∆〈φ(x, t)〉 − 〈φ(x, t)〉 − 〈φ3(x, t)〉. (14)
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 〈φ3〉, but it
would involve a higher order term 〈φ5〉, and so on. Each
equation contains higher moments and therefore the hi-
erarchy does not close.
To close the equation for µ(t) we assume a parametric
PDF of the form P [φ;α] ∝ exp(−H [φ;α]/kT ) where
H [φ;α] = H [φ] + α
φ(x) dx
is a perturbation of the Hamiltonian (13) by a term pro-
portional to the moment variable ξ[φ] = (1/V )
φ(x) dx.
This is a special case of a general “entropy-based” clo-
sure prescription proposed by Eyink and Levermore [16].
This closure scheme guarantees that α(t) → 0 at long
times and therefore the PDF ansatz P [φ;α(t)] relaxes to
the correct stationary distribution P∗[φ] of the stochastic
process. The determination of the parameter α given the
moment µ is here accomplished by Legendre transform
α = argmaxα[αµ− F (α)], (15)
where the “moment-generating function” F (α) =
log〈exp[α
φ(x) dx]〉∗ and 〈·〉∗ denotes average with re-
spect to the invariant measure P∗[φ]. The numerical op-
timization required for the Legendre transform is well-
suited to gradient descent algorithms such as the conju-
gate gradient method, since
(∂/∂α)[αµ− F (α)] = µ− µ(α),
where µ(α) = 〈ξ〉α is the average of the moment-function
in the PDF ansatz P [φ;α]. In simple cases, F (α) and
µ(α) = F ′(α) may be given by closed analytical expres-
sions. If not, then both of these averages may be deter-
mined together by Monte Carlo sampling techniques.
In the numerical calculations below, we discretize equa-
tion (12) using a forward Euler-Maruyama stochastic in-
tegrator and 3-point stencil for the Laplacian (other dis-
cretizations are possible).
φ(x, t + δt) = φ(x, t)− δt[φ(x, t) + φ3(x, t)] + (16)
(δx)2
[φ(x + δx, t)− 2φ(x, t) + φ(x− δx, t)] +
2kT (δt/δx)N(x, t)
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 ∝ exp(−Hδ/kT ) with discrete Hamiltonian
〈x,x′〉
(φ(x) − φ(x′))2 (17)
φ2(x) +
φ4(x)]
where 〈x, x′〉 are nearest-neighbor pairs. The closure
ansatz can be adopted in the consistently discretized form
Pδ[φ;α] ∝ exp(−Hδ[φ;α]/kT ) where
Hδ[φ;α] = Hδ[φ] + α
δxφ(x).
In this numerical experiment, we integrate an N =
1000 member ensemble of solutions of equation (17),
and measure the ensemble-averaged, global magnetiza-
tion µ(t) = 〈φ(t)〉 = (1/V )
〈φ(x, t)〉 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 inte-
gration scheme works as follows: Suppose we are given
the parameter α(t) at time t. The mean µ(t) is first cal-
culated from the parametric ensemble at time t by Monte
Carlo sampling. Next all N samples are integrated over
a short time-step δt to create a time-advanced ensemble.
From this ensemble µ(t + δt) is calculated, which yields
an estimate of the local time derivative.
µ̇app(t) = [µ(t+ δt)− µ(t)]/δt.
A large, projective Euler time-step of the moment aver-
age is then taken via
µ(t+∆t) = µ(t) + ∆t µ̇app(t).
The parameter is finally updated by using the Legendre
transform inversion to obtain α(t+∆t) from the known
value µ(t + ∆t). The cycle may now be repeated to in-
tegrate the closure equations by successive time-steps of
length ∆t.
A critical issue in general application of projective inte-
gration is the criterion to determine the projective time-
step ∆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/slow-
est “fast group” eigenvalue). Variants of the approach
have been developed for problems with several gaps in
their spectrum [18]. 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 poste-
riori 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 hypothe-
sis testing (are the data consistent with a local linear
model?) must be considered. These are important re-
search issues that are currently explored by several re-
search groups. Nevertheless, the main factor in computa-
tional 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 δt and
the “coarse-scale” projective integration step ∆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 δt and for ∆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 δt = 0.0004. On
the other hand, for projective integration of the closure
equation a time step ∆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 equa-
tions 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 Equa-
tion (14) as long as the statistics are homogeneous and
the Laplacian term vanishes.
CONCLUSIONS
In this paper, we have described how one can combine
recently developed equation-free methods with statisti-
cal moment closures to model nonlinear problems. With
this method we can numerically integrate complex non-
linear 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-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
FIG. 1: Mean total field as a function of time. Line (sym-
bols): traditional (coarse projective) integration, respectively.
See the text for a description of the stepsize selection.
FIG. 2: Comparison of the time dependent PDF’s of the local
field φ(x, t) for the exact solution (blue) and for the projective
integration / closure solution (red).
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 statis-
tical moment closure approach. Since little or no an-
alytic work is required, the sophistication of statistical
moment closures can greatly enhanced beyond Gaussian
PDF ansätze. The “practical usefulness” criterion for
parametric PDF models that they permit analytical cal-
culations 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 ap-
proaches. One will likely have to use nonlocal, acceler-
ated sampling methods. One will also likely employ the
latest in adaptive time and adaptive mesh methods to
optimize performance for large-scale problems.
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 un-
der 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.
[1] T. Schlick, R. D. Skeel, A. T. Brunger, L. V. Kale, J. A.
Board, Jr. , J. Hermans, and K. Schulten J. Comp. Phys.
151, 9, (1999).
[2] M. Karplus and J. A. McCammon, Nature, Structural
and Molecular Biology, 9 , 646, (2002).
[3] P. G. Debenedetti and F. H. Stillinger, Nature, 410, 259,
(2001).
[4] D. T. Gillespie, J. Phys. Chem., 81 , 2340, (1977).
[5] D. J. Wilkinson, Stochastic Modeling for Systems Biol-
ogy, Chapman & Hall / CRC Press, Boca Raton, (2006).
[6] A. J. Majda and X. Wang, Nonlinear Dynamics and Sta-
tistical Theories for Basic Geophysical Flows, Cambridge
Univeristy Press, Cambridge UK, 2006
[7] J. P. Hansen and I. R. MacDonald, Theory of Simple
Liquids, Academic, New York, (1986).
[8] S. B. Pope, Turbulent Flows, Cambridge University
Press, Cambridge, UK, (2000).
[9] A. J. Chorin, O. H. Hald, and R. Kupferman, Proceedings
of the National Academy of Sciences of the United States
of America, 97, 2968, (2000).
[10] H. Chen, S. Chen, and R. H. Kraichnan, Phys. Rev. Lett.,
63, 2657–2660, 1989.
[11] C. Yeung, Y. Oono, and A. Shinozaki Phys. Rev. E, 49,
2693 (1994)
[12] G. L. Eyink, Phys. Rev. E 54 (1996) 3419–3435.
[13] C.D. Levermore, J. Stat. Phys., 86 (1996), 1021–1065.
[14] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G.
Kevrekidis, O. Runborg and K. Theodoropoulos, Comm.
Math. Sciences 1(4) pp.715-762 (2003); S. L. Lee and C.
W. Gear, J. Comp. App. Math., 201, 258, (2007).
[15] I. G. Kevrekidis, C. William Gear and G. Hummer,
A.I.Ch.E Journal, 50(7) pp.1346-1354 (2004)
[16] G. L. Eyink and C. D. Levermore (preprint) Entropy-
Based Closures of Nonlinear Stochastic Dynamics. sub-
mitted to ”Communications in Mathematical Sciences”
(2006).
[17] A. J. Bray, Adv. in Phys., 43, 357, (1994)
[18] C. W. Gear and I. G. Kevrekidis, J. Comp. Phys., 187,
95, (2003)
[19] In the case of [10] the dynamics is an advection-reaction-
diffusion equation for a scalar concentration field X(t) =
{θ(x, t) : x ∈ Rd}. The moment functions are the “fine-
grained PDF” ξϑ,x[X, t] = δ(θ(x, t) − ϑ), labelled by
space point x and scalar value ϑ. The moment average
µϑ,x(t) = 〈δ(θ(x, t) − ϑ)〉 is the 1-point PDF p(ϑ;x, t)
which gives the distribution of scalar values ϑ at space-
time point (x, t). The parametric model P [X;α, t] is
the distribution over scalar fields obtained by the ansatz
θ(x, t) = X(θ0(x, t),x, t) where θ0(x, t) is a reference ran-
dom field of known (Gaussian) statistics and X(·,x, t) :
R → R is a “mapping function”. The latter function is
the “parameter” αϑ0,x(t) = X(ϑ0,x, t) which determines
(and is determined by) the “moment” µϑ,x(t) from the re-
lation p(X(ϑ0,x, t);x, t)|∂X/∂ϑ0| = p0(ϑ0,x, t). Here p0
is the 1-point PDF of the reference Gaussian field θ0(x, t).
The approach of [11] is similar. The problem is phase-
ordering dynamics as given, for example, by our equation
(12) and X(t) = {φ(x, t) : x ∈ Rd}. The moment func-
tions are the quadratic products ξr[X, t] = φ(r, t)φ(0, t),
labelled by the displacement r ∈ Rd and the moment av-
erages µr(t) are the spatial correlation function C(r, t).
The parametric model P [X;α, t] is the distribution ob-
tained by the ansatz φ(x, t) = f(u(x, t)) where u(x, t) is a
homogeneous Gaussian random field with mean zero and
covariance G(r, t) = 〈u(r, t)u(0, t)〉 and f(z) is the sta-
tionary planar interface solution of the TDGL equation
(12). In this case, it is the auxiliary correlation function
G(r, t) which plays the role of the “parameter” αr(t). It
is shown in [11] for various cases how this function may
be uniquely related to the “moment” µr(t) = C(r, t).

  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.

Equation-Free Implementation of Statistical Moment Closures
Francis J. Alexander and Gregory Johnson
Los Alamos National Laboratory, P.O.Box 1663,
Los Alamos, NM, 87545.
Gregory L. Eyink
Department of Mathematical Sciences
Johns Hopkins University
Baltimore, MD 21218
Ioannis G. Kevrekidis
Department of Chemical Engineering and PACM
Princeton University
Princeton, NJ 08544
We present a general numerical scheme for the practical implementation of statistical moment clo-
sures 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 typi-
cally 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.
PACS numbers:
INTRODUCTION
Accurate, fast simulations of complex, large-scale, non-
linear systems remain a challenge for computational sci-
ence and engineering, despite extraordinary advances in
computing power. Examples range from molecular dy-
namics simulations of proteins [1], [2] and glasses [3], to
stochastic simulations of cellular biochemistry [4, 5], to
global-scale, geophysical fluid dynamics [6]. Often for the
systems under consideration there is no obvious scale sep-
aration, and their many degrees of freedom are strongly
coupled. The complex and multiscale nature of these pro-
cesses therefore makes them extremely difficult to model
numerically. To make matters worse, one is often in-
terested not in a single, time-dependent solution of the
equations governing these processes, but rather in ensem-
bles of solutions consisting of multiple realizations (e.g.,
sampling noise, initial conditions, and/or uncertain pa-
rameters). 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 Liou-
ville (if deterministic)/Kolmogorov (if stochastic) equa-
tions. These are partial differential equations in a state
space of high (possibly infinite) dimension. A straightfor-
ward discretization of the Liouville / Kolmogorov equa-
tions is therefore impractical. An ensemble approach to
solving these equations can be taken; however, quite of-
ten, 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 real-
ization may be prohibitive.
The traditional approach to making these prob-
lems computationally tractable is to replace the Liou-
ville/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 phys-
ical insight. For example, in the superposition approxi-
mation and its extensions [7] for dense liquids and plas-
mas, 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 [8]
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 ensem-
ble 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
http://arxiv.org/abs/0704.0804v1
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 in-
formation (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 turbu-
lence [10], phase-ordering dynamics [11] 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 [12].
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 com-
puter 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 up-
dated, 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 clo-
sure equations analytically. This is especially likely when
PDF’s are not Gaussian, which is the case for most use-
ful closures. Monte Carlo or other numerical methods
may be needed in order to evaluate integrals for the mo-
ments [13]. 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 numer-
ical code with an inaccessible source. Clearly, a need ex-
ists 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 [14] [15]. 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 clo-
sure equations that we may not explicitly have. The
equation-free approach extends the applicability of sta-
tistical 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 clo-
sures. 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 sys-
tem, the stochastic Ginzburg-Landau (GL) equations us-
ing a particular PDF-based closure scheme, the entropy
method of Eyink and Levermore [16]. We conclude with
a discussion of closure quality, computational issues, and
the application of our approach to large-scale systems.
PDF-BASED MOMENT CLOSURES
We consider the very general class of dynamical sys-
tems, including maps, formally represented by
Ẋ = U(X(t),N(t), t) (1)
Xt+1 = Ut(Xt,Nt) (2)
where N(t) is a stochastic process with prescribed statis-
tics. The stochastic component arises from unknown pa-
rameters, 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
N(t) is a Markov process (Brownian motion, Poisson pro-
cess, etc.) and—more specifically still—Itô stochastic dif-
ferential equations of the form:
dX = U(X, t)dt +
2S(X, t)dW(t). (3)
The deterministic component of the state, X, is gov-
erned by the continuously differentiable vector field, U :
N × R → RN . For many problems of interest (e.g., cli-
mate) U is a highly nonlinear function. The noise com-
ponent is modeled by the standard mean 0, covariance
matrix I Wiener process, W ∈ RN , possibly modulated
by a state-dependent matrix S : RN×R → RN×N . Equa-
tion (3) encompasses a wide class of systems including
deterministic (S = 0) ones.
In many cases one is interested in knowing the low
order statistics of equation (3), for example an instanta-
neous mean value or possibly multi-point covariance of
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 (X, t):
∂tP = L∗(t)P, (4)
where P satisfies the conditions: P (X, t) ≥ 0, and
P (X, t) dX = 1, and where L∗ is the generator of the
Markov process. In the case of equation (3) this operator
takes the form
L∗(t)ψ(X) = −∇X·(U(X, t)ψ(X))+∇2X : (D(X, t)ψ(X)).
The forward Kolmogorov equation then becomes a
Fokker-Planck equation
∂tP +∇X · (UP ) = ∇2X : (DP ) (6)
where D(X, t) = S(X, t)S(X, t)T is the nonnegative-
definite diffusion matrix arising from the noise term. Un-
like the original dynamical equation (3), the forward Kol-
morogov equation (FKE) is both linear and determinis-
tic. Dealing with it, therefore, has apparent advantages
over the original ensemble of stochastic systems simu-
lations. 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 reduc-
tion should simplify the computation as much as pos-
sible while retaining fidelity to the original dynamical
processes. The reduction proceeds by taking moments of
the FKE with respect to a vector-valued function ξ(X, t)
from RN × R+ → RM . The ξ selected should include
the relevant variables in the system (slow modes, con-
served quantities, etc.). The moments µ(t) of ξ(X, t) are
defined by
µ(t) =
ξ(X, t)P (X, t)dX (7)
and give rise to
µ̇(t) =
ξ̇(X, t)P (X, t)dX, (8)
where
ξ̇(X, t) = ∂tξ(X, t) + L(t)ξ(X, t) (9)
and L is the adjoint of L∗ or the backward Kolmogorov
operator. The result (8) can be obtained by averaging
over an ensemble of realizations of the stochastic dynam-
ics (3). In general, however, (8) is not a closed equa-
tion for the moments, µ. One can close this equation by
choosing a PDF, P (X, t,µ), which itself is a function of
the moments µ.
µ̇(t) = V(µ, t) ≡
ξ̇(X, t)P (X, t,µ)dX. (10)
Alternatively, one can select a family of probability den-
sities P (X, t,α), specified by parameters α = α(µ, t)
rather than directly by the moments µ. 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 param-
eters 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 (X, t,α(µ)(t)) will be close to the exact solution of
the Liouville/Kolmogorov equation (4). The mapping
closure approach of Chen et al [10] and the Gaussian
mapping method of Yeung et al. [11] are based on this
type of parametric PDF closure [19]. In fact, perhaps
the most familiar application of the parametric approach
is the use of the Rayleigh-Ritz method in quantum me-
chanical calculations. This is the essential approach of
our paper.
EQUATION-FREE COMPUTATION
Although we now have obtained a closed moment equa-
tion (equation 10), we still need to determine the dynam-
ical vector field 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
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 collabora-
tors [14] 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 for-
mulae for the closure equations. Instead, one must only
be able to sample an ensemble of system states X dis-
tributed according to the closure ansatz P (X, t;α) and
then evolve each of these via equation (3) for short inter-
vals 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 avail-
able which, given a specific state variable X(t) as input,
returns the value of the state X(t + δt) after a short
time δt. The ensemble of systems, each of which satis-
fies equation (3), is evolved over a time interval δt. The
moments/parameters µ or α are determined at the be-
ginning and end of this interval and the time derivative
µ̇ 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 equa-
tions, namely V(µ, t).
Coarse projective forward Euler (arguably the simplest
of equation-free algorithms) which we will use below il-
lustrates the approach succinctly: Starting from a set
of coarse-grained initial conditions specified by moments
µ(t) we first (a) lift to a consistent fine scale descrip-
tion, that is, sample the PDF ansatz P (X, t;α(t)) to
generate ensembles of initial conditions X for equation
(3) consistent with the set µ(t); (b) starting with these
consistent initial conditions we evolve the fine scale de-
scription for a (relatively short) time δt; we subsequently
restrict back to coarse observables by evaluating the mo-
ments µ(t+ δt) as ensemble-averages and (d) use the re-
sults to estimate locally the time derivative dµ/dt. This
is precisely the right hand-side of the explicitly unavail-
able closure, obtained not through a closed form formula,
but rather through short, judicious computational exper-
iments 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:
µ(t+∆t) = µ(t) + ∆t
µ(t+ δt)− µ(t)
. (11)
The procedure then repeats itself: lifting, fine scale evo-
lution, restriction, estimation, and then (connecting with
continuum traditional numerical analysis) a new for-
ward Euler step. Beyond coarse projective forward Eu-
ler, many other coarse initial-value solvers (e.g. coarse
projective Adams-Bashforth, and even implicit coarse
solvers) have been implemented; the stability and accu-
racy study of such algorithms is progressing [14]. These
developments allow us to construct a nonintrusive imple-
mentation of PDF moment closures, nonintrusive in the
sense that we compute with the closures without explic-
itly 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 µ(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 X from
any chosen distribution P0(X), evolve each of these real-
izations 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 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 appli-
cations 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 µ(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 t0 to a final time t0+T
would require O(T/δt) integration steps, where the time-
step δt is required to be very small by the intrinsic stiff-
ness of the micro-dynamics. In the closure approach, the
evolution of the moment equations (10) from time t0 to
time t0+T requires only O(T/∆t) integration steps, with
(hopefully) ∆t ≫ δt. Each of these closure integration
steps by an increment ∆t requires in the equation-free
approach just one (or just a few) fine-scale integration
step by an increment δt. Thus, there is an over-all savings
by a (hopefully) large factor O(∆t/δt). This crude esti-
mate 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µ/dt from full direct
simulation to the time that we can (linearly or even poly-
nomially) extrapolate µ(t) in time. It is also worth not-
ing that a variety of additional computational tasks, be-
yond 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.
A NUMERICAL EXAMPLE
We illustrate here the equation-free implementation
of moment-closures for a canonical equation of phase-
ordering kinetics [17], the stochastic time-dependent
Ginzburg-Landau (TDGL) equation in one spatial di-
mension. This is written as
∂φ(x, t)
= D∆φ(x, t) − V ′(φ(x, t)) + η(x, t) (12)
where φ(x, t) represents a local order parameter, e.g. a
magnetization. The noise has mean zero and covariance
〈η(x, t)η(x′, t′)〉 = 2kT δ(x−x′)δ(t− t′). The potential V
shall be chosen as
V (φ) =
to represent a single quartic/quadratic well. This
stochastic dynamics has an invariant measure which is
formally of Hamiltonian form P∗[φ] ∝ exp(−H [φ]/kT )
where
H [φ] =
D|∇φ(x)|2 + V (φ(x))] dx. (13)
The Gibbsian measure P∗[φ] is approached at long times
for any random distribution P0[φ] of initial states.
One of the simplest dynamical quantities of interest is
the bulk magnetization φ(t) = (1/V )
φ(x, t)dx, where
V is the total volume. If the initial statistics are space-
homogeneous, then the ensemble average µ(t) = 〈φ(t)〉 is
also given by µ(t) = 〈φ(x, t)〉 for any space point x. Equa-
tion (12) leads to a hierarchy of equations for statistical
moments of φ(x, t). For example, the first moment satis-
fies the equation
∂〈φ(x, t)〉
= ∆〈φ(x, t)〉 − 〈φ(x, t)〉 − 〈φ3(x, t)〉. (14)
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 〈φ3〉, but it
would involve a higher order term 〈φ5〉, and so on. Each
equation contains higher moments and therefore the hi-
erarchy does not close.
To close the equation for µ(t) we assume a parametric
PDF of the form P [φ;α] ∝ exp(−H [φ;α]/kT ) where
H [φ;α] = H [φ] + α
φ(x) dx
is a perturbation of the Hamiltonian (13) by a term pro-
portional to the moment variable ξ[φ] = (1/V )
φ(x) dx.
This is a special case of a general “entropy-based” clo-
sure prescription proposed by Eyink and Levermore [16].
This closure scheme guarantees that α(t) → 0 at long
times and therefore the PDF ansatz P [φ;α(t)] relaxes to
the correct stationary distribution P∗[φ] of the stochastic
process. The determination of the parameter α given the
moment µ is here accomplished by Legendre transform
α = argmaxα[αµ− F (α)], (15)
where the “moment-generating function” F (α) =
log〈exp[α
φ(x) dx]〉∗ and 〈·〉∗ denotes average with re-
spect to the invariant measure P∗[φ]. The numerical op-
timization required for the Legendre transform is well-
suited to gradient descent algorithms such as the conju-
gate gradient method, since
(∂/∂α)[αµ− F (α)] = µ− µ(α),
where µ(α) = 〈ξ〉α is the average of the moment-function
in the PDF ansatz P [φ;α]. In simple cases, F (α) and
µ(α) = F ′(α) may be given by closed analytical expres-
sions. If not, then both of these averages may be deter-
mined together by Monte Carlo sampling techniques.
In the numerical calculations below, we discretize equa-
tion (12) using a forward Euler-Maruyama stochastic in-
tegrator and 3-point stencil for the Laplacian (other dis-
cretizations are possible).
φ(x, t + δt) = φ(x, t)− δt[φ(x, t) + φ3(x, t)] + (16)
(δx)2
[φ(x + δx, t)− 2φ(x, t) + φ(x− δx, t)] +
2kT (δt/δx)N(x, t)
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 ∝ exp(−Hδ/kT ) with discrete Hamiltonian
〈x,x′〉
(φ(x) − φ(x′))2 (17)
φ2(x) +
φ4(x)]
where 〈x, x′〉 are nearest-neighbor pairs. The closure
ansatz can be adopted in the consistently discretized form
Pδ[φ;α] ∝ exp(−Hδ[φ;α]/kT ) where
Hδ[φ;α] = Hδ[φ] + α
δxφ(x).
In this numerical experiment, we integrate an N =
1000 member ensemble of solutions of equation (17),
and measure the ensemble-averaged, global magnetiza-
tion µ(t) = 〈φ(t)〉 = (1/V )
〈φ(x, t)〉 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 inte-
gration scheme works as follows: Suppose we are given
the parameter α(t) at time t. The mean µ(t) is first cal-
culated from the parametric ensemble at time t by Monte
Carlo sampling. Next all N samples are integrated over
a short time-step δt to create a time-advanced ensemble.
From this ensemble µ(t + δt) is calculated, which yields
an estimate of the local time derivative.
µ̇app(t) = [µ(t+ δt)− µ(t)]/δt.
A large, projective Euler time-step of the moment aver-
age is then taken via
µ(t+∆t) = µ(t) + ∆t µ̇app(t).
The parameter is finally updated by using the Legendre
transform inversion to obtain α(t+∆t) from the known
value µ(t + ∆t). The cycle may now be repeated to in-
tegrate the closure equations by successive time-steps of
length ∆t.
A critical issue in general application of projective inte-
gration is the criterion to determine the projective time-
step ∆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/slow-
est “fast group” eigenvalue). Variants of the approach
have been developed for problems with several gaps in
their spectrum [18]. 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 poste-
riori 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 hypothe-
sis testing (are the data consistent with a local linear
model?) must be considered. These are important re-
search issues that are currently explored by several re-
search groups. Nevertheless, the main factor in computa-
tional 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 δt and
the “coarse-scale” projective integration step ∆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 δt and for ∆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 δt = 0.0004. On
the other hand, for projective integration of the closure
equation a time step ∆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 equa-
tions 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 Equa-
tion (14) as long as the statistics are homogeneous and
the Laplacian term vanishes.
CONCLUSIONS
In this paper, we have described how one can combine
recently developed equation-free methods with statisti-
cal moment closures to model nonlinear problems. With
this method we can numerically integrate complex non-
linear 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-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
FIG. 1: Mean total field as a function of time. Line (sym-
bols): traditional (coarse projective) integration, respectively.
See the text for a description of the stepsize selection.
FIG. 2: Comparison of the time dependent PDF’s of the local
field φ(x, t) for the exact solution (blue) and for the projective
integration / closure solution (red).
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 statis-
tical moment closure approach. Since little or no an-
alytic work is required, the sophistication of statistical
moment closures can greatly enhanced beyond Gaussian
PDF ansätze. The “practical usefulness” criterion for
parametric PDF models that they permit analytical cal-
culations 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 ap-
proaches. One will likely have to use nonlocal, acceler-
ated sampling methods. One will also likely employ the
latest in adaptive time and adaptive mesh methods to
optimize performance for large-scale problems.
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 un-
der 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.
[1] T. Schlick, R. D. Skeel, A. T. Brunger, L. V. Kale, J. A.
Board, Jr. , J. Hermans, and K. Schulten J. Comp. Phys.
151, 9, (1999).
[2] M. Karplus and J. A. McCammon, Nature, Structural
and Molecular Biology, 9 , 646, (2002).
[3] P. G. Debenedetti and F. H. Stillinger, Nature, 410, 259,
(2001).
[4] D. T. Gillespie, J. Phys. Chem., 81 , 2340, (1977).
[5] D. J. Wilkinson, Stochastic Modeling for Systems Biol-
ogy, Chapman & Hall / CRC Press, Boca Raton, (2006).
[6] A. J. Majda and X. Wang, Nonlinear Dynamics and Sta-
tistical Theories for Basic Geophysical Flows, Cambridge
Univeristy Press, Cambridge UK, 2006
[7] J. P. Hansen and I. R. MacDonald, Theory of Simple
Liquids, Academic, New York, (1986).
[8] S. B. Pope, Turbulent Flows, Cambridge University
Press, Cambridge, UK, (2000).
[9] A. J. Chorin, O. H. Hald, and R. Kupferman, Proceedings
of the National Academy of Sciences of the United States
of America, 97, 2968, (2000).
[10] H. Chen, S. Chen, and R. H. Kraichnan, Phys. Rev. Lett.,
63, 2657–2660, 1989.
[11] C. Yeung, Y. Oono, and A. Shinozaki Phys. Rev. E, 49,
2693 (1994)
[12] G. L. Eyink, Phys. Rev. E 54 (1996) 3419–3435.
[13] C.D. Levermore, J. Stat. Phys., 86 (1996), 1021–1065.
[14] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G.
Kevrekidis, O. Runborg and K. Theodoropoulos, Comm.
Math. Sciences 1(4) pp.715-762 (2003); S. L. Lee and C.
W. Gear, J. Comp. App. Math., 201, 258, (2007).
[15] I. G. Kevrekidis, C. William Gear and G. Hummer,
A.I.Ch.E Journal, 50(7) pp.1346-1354 (2004)
[16] G. L. Eyink and C. D. Levermore (preprint) Entropy-
Based Closures of Nonlinear Stochastic Dynamics. sub-
mitted to ”Communications in Mathematical Sciences”
(2006).
[17] A. J. Bray, Adv. in Phys., 43, 357, (1994)
[18] C. W. Gear and I. G. Kevrekidis, J. Comp. Phys., 187,
95, (2003)
[19] In the case of [10] the dynamics is an advection-reaction-
diffusion equation for a scalar concentration field X(t) =
{θ(x, t) : x ∈ Rd}. The moment functions are the “fine-
grained PDF” ξϑ,x[X, t] = δ(θ(x, t) − ϑ), labelled by
space point x and scalar value ϑ. The moment average
µϑ,x(t) = 〈δ(θ(x, t) − ϑ)〉 is the 1-point PDF p(ϑ;x, t)
which gives the distribution of scalar values ϑ at space-
time point (x, t). The parametric model P [X;α, t] is
the distribution over scalar fields obtained by the ansatz
θ(x, t) = X(θ0(x, t),x, t) where θ0(x, t) is a reference ran-
dom field of known (Gaussian) statistics and X(·,x, t) :
R → R is a “mapping function”. The latter function is
the “parameter” αϑ0,x(t) = X(ϑ0,x, t) which determines
(and is determined by) the “moment” µϑ,x(t) from the re-
lation p(X(ϑ0,x, t);x, t)|∂X/∂ϑ0| = p0(ϑ0,x, t). Here p0
is the 1-point PDF of the reference Gaussian field θ0(x, t).
The approach of [11] is similar. The problem is phase-
ordering dynamics as given, for example, by our equation
(12) and X(t) = {φ(x, t) : x ∈ Rd}. The moment func-
tions are the quadratic products ξr[X, t] = φ(r, t)φ(0, t),
labelled by the displacement r ∈ Rd and the moment av-
erages µr(t) are the spatial correlation function C(r, t).
The parametric model P [X;α, t] is the distribution ob-
tained by the ansatz φ(x, t) = f(u(x, t)) where u(x, t) is a
homogeneous Gaussian random field with mean zero and
covariance G(r, t) = 〈u(r, t)u(0, t)〉 and f(z) is the sta-
tionary planar interface solution of the TDGL equation
(12). In this case, it is the auxiliary correlation function
G(r, t) which plays the role of the “parameter” αr(t). It
is shown in [11] for various cases how this function may
be uniquely related to the “moment” µr(t) = C(r, t).
