scispace - formally typeset

Book ChapterDOI

A Multilevel Monte Carlo Asymptotic-Preserving Particle Method for Kinetic Equations in the Diffusion Limit

AbstractWe propose a multilevel Monte Carlo method for a particle-based asymptotic-preserving scheme for kinetic equations. Kinetic equations model transport and collision of particles in a position-velocity phase-space. With a diffusive scaling, the kinetic equation converges to an advection-diffusion equation in the limit of zero mean free path. Classical particle-based techniques suffer from a strict time-step restriction to maintain stability in this limit. Asymptotic-preserving schemes provide a solution to this time step restriction, but introduce a first-order error in the time step size. We demonstrate how the multilevel Monte Carlo method can be used as a bias reduction technique to perform accurate simulations in the diffusive regime, while leveraging the reduced simulation cost given by the asymptotic-preserving scheme. We describe how to achieve the necessary correlation between simulation paths at different levels and demonstrate the potential of the approach via numerical experiments.

Topics: Monte Carlo method (62%), Limit (mathematics) (51%)

Summary (2 min read)

1 Introduction

  • Kinetic equations, modeling particle behavior in a position-velocity phase space, occur in many domains.
  • Deterministic methods become prohibitively expensive for higher dimensional applications.
  • In the diffusive scaling, there are only of two works [9, 16] so far, to the best of their knowledge.
  • In Section 3, the authors cover the multilevel Monte Carlo method that is the core contribution of this paper.

2.1 Model equation in the diffusive limit

  • (6) Page:3 job:MLMCforAP macro:svmult.cls date/time:15-Jan-2019/18:31 Equation (6) is also known as the Goldstein-Taylor model, and can be solved using a particle scheme.
  • Equation (6) is then solved via operator splitting as 1.
  • Transport step. (9) To maintain stability, this approximation requires a time step restriction ∆ t = O(ε2), leading to unacceptably high computational costs.

2.2 Asymptotic-preserving Monte Carlo scheme

  • Discretizing this equation, using operator splitting as above, again leads to a Monte Carlo scheme.
  • For each particle Xp and for each time step n, one time step now consists of a transport-diffusion and a collision step: 1. Transport-diffusion step. (13).
  • For more details, the authors refer the reader to [16].

3.1 Method and notation

  • That is a choice the authors make for notational convenience, and is not essential for the method they present.
  • On the one hand, a small time step reduces the bias of the simulation of each sampled trajectory, and thus of the estimated quantity of interest.
  • The key idea behind the Multilevel Monte Carlo method [17] is to generate a sequence of estimates with varying discretization accuracy and a varying number of realizations.
  • The estimators (17) estimate the bias induced by sampling with a simulation time step size ∆ t`−1 by comparing the sample results with a simulation using a time step size ∆ t`.

3.2.1 Coupled trajectories and notation

  • The differences in (17) will only have low variance if the simulated paths Xn,m∆ t`,p and Xn∆ t`−1,p are correlated.
  • In each time step using the asymptotic-preserving particle scheme (11)–(13), there are two sources of stochastic behavior.
  • Particle trajectories can be coupled by separately correlating the random numbers used for each individual particle in the transport-diffusion and collision phase of each time step.
  • The authors present the complete algorithm in Section 3.2.4.

3.2.2 Coupling the transport-diffusion phase

  • The authors observe that the paths have similar behavior, i.e., if the fine simulation tends towards negative values, so does the coarse simulation and vice versa.
  • Still, there is an observable difference between them.
  • This is due to the bias caused by the paths having different diffusion coefficients.

3.2.3 Coupling the collision phase

  • While correlating the Brownian paths is relatively straightforward, the coupling of the velocities in the collision phase is more involved.
  • Afterwards, when the authors decide to perform a collision both at level ` and `−1, they will correlate the new velocities generated in both simulations.
  • The maximum of a set of uniformly distributed random number is not uniformly distributed.
  • This approach to selecting the sign of V n+1p,∆ t`−1 means that the velocities going into the next time step will have the same sign.
  • This is one source of the bias that the authors want to estimate using the multilevel Monte Carlo method.

4 Experimental Results

  • The authors will now demonstrate the viability of the suggested approach through some numerical experiments.
  • The authors will simulate the model given by (10), using the multilevel Monte Carlo method to estimate a selected quantity of interest, which is the expected Page:10 job:.
  • The ensemble of particles is initialized at the origin with equal probability of having a left and right velocity.
  • When discussing results the authors will replace the full expression for a sample of the quantity of interest, based on an arbitrary particle p, F(XN,0∆ t`,p), with the symbol F̀ to simplify notation.

4.1 Model correlation behavior

  • The authors fix the number of samples per difference estimator at 100 000.
  • This matches the weak convergence order of the Euler-Maruyama scheme, used to simulate the model (11)–(13), as well as the expected behavior from the time step dependent bias in the asymptotic-preserving model.
  • This means that the existing general theory for multilevel Monte Carlo methods [18] concerning, e.g. samples per level, convergence criteria and conditions for adding levels, can be applied in this regime.
  • This explains the decrease in the means of the differences and in the variance of the differences for both small and large ∆ t which is clearly visible in Figure 6.

5 Conclusion

  • The authors have derived a new multilevel scheme for asymptotic-preserving particle schemes of the form given in (10).
  • The authors have demonstrated that this scheme has interesting convergence behavior as the time step is refined, which is apparent in the expected value and variance of sampled differences of the quantity of interest.
  • The authors have proposed a strategy for selecting the coarsest level in the multilevel Monte Carlo method in this context, and shown that this simulation approach gives a significant speedup over a classical Monte Carlo simulation.
  • More complicated models including, for example, absorption terms can also be studied.

Did you find this useful? Give us your feedback

...read more

Content maybe subject to copyright    Report

A Multilevel Monte Carlo
Asymptotic-Preserving Particle Method for
Kinetic Equations in the Diffusion Limit
Emil Løvbak, Giovanni Samaey, and Stefan Vandewalle
Abstract
We propose a multilevel Monte Carlo method for a particle-based asymptotic-
preserving scheme for kinetic equations. Kinetic equations model transport and
collisions of particles in a position-velocity phase-space. With a diffusive scaling,
the kinetic equation converges to an advection-diffusion equation in the limit of zero
mean free path. Classical particle-based techniques suffer from a strict time-step
restriction to maintain stability in this limit. Asymptotic-preserving schemes provide
a solution to this time step restriction, but introduce a first-order error in the time step
size. We demonstrate how the multilevel Monte Carlo method can be used as a bias
reduction technique to perform accurate simulations in the diffusive regime, while
leveraging the reduced simulation cost given by the asymptotic-preserving scheme.
We describe how to achieve the necessary correlation between simulation paths at
different levels and demonstrate the potential of the approach via computational
experiments.
1 Introduction
Kinetic equations, modeling particle behavior in a position-velocity phase space,
occur in many domains. Examples are plasma physics [
4
], bacterial chemotaxis [
33
]
and computational fluid dynamics [
32
]. Many of these applications exhibit a strong
time-scale separation, leading to an unacceptably high simulation cost [
7
]. However,
one typically is only interested in computing the evolution of some macroscopic
quantities of interest. These are usually some moments of the particle distribution,
which can be computed as averages over velocity space. The time-scale at which
these quantities of interest change is often much slower than the time-scale governing
Emil Løvbak · Giovanni Samaey · Stefan Vandewalle
KU Leuven, Department of Computer Science, NUMA Section, Celestijnenlaan 200A
box 2402, 3001 Leuven, Belgium e-mail: emil.loevbak@cs.kuleuven.be, e-mail: gio-
vanni.samaey@cs.kuleuven.be, e-mail: stefan.vandewalle@cs.kuleuven.be
1

2 Emil Løvbak, Giovanni Samaey, and Stefan Vandewalle
the particle dynamics. The nature of the macroscopic dynamics depends on the
scaling of the problem, which can be either hyperbolic or diffusive [15].
The model problem considered in this work is a one-dimensional kinetic equation
of the form
t
f (x,v,t) + v
x
f (x,v,t) = Q ( f (x,v,t)), (1)
where
f (x,v,t)
represents the distribution
f
of particles as a function of position
x R
and velocity
v R
as it evolves in time
t
. The left-hand side of equation
(1)
represents
transport, while
Q( f (x,v,t))
is a collision operator that results in discontinuous
velocity changes. As the collision operator, we take the BGK model [
3
], which
represents linear relaxation to an equilibrium distribution that only depends on the
position density
ρ(x,t) =
Z
f (x,v,t)dv. (2)
We introduce a parameter
ε
that represents the mean free path. When decreasing
ε
,
the average time between collisions decreases. In this paper, we consider the diffusive
scaling. In that case, we simultaneously increase the time scale at which we observe
evolution of the particle distribution, arriving at
ε
t
f (x,v,t) + v
x
f (x,v,t) =
1
ε
(ρ(x,t) f (x, v,t)). (3)
It can be shown that when taking the limit
ε 0
, the behavior of equations of the
form (5) is fully described by the diffusion equation [25]
t
ρ(x,t) =
xx
ρ(x,t). (4)
The simulation of kinetic equations can be done with deterministic methods,
solving the partial differential equation (PDE) that describes evolution of the particle
distribution in position-velocity phase space. Alternatively, one can use stochastic
methods that simulate a large number of particle trajectories. Deterministic methods
become prohibitively expensive for higher dimensional applications. Particle-based
methods do not suffer from this curse of dimensionality, at the expense of introducing
a statistical error in the computed solution. The issue of time-scale separation is
present in both deterministic and stochastic methods.
One way to avoid th issue of time-scale separation, is through the use of
asymptotic-preserving methods, which aim at reproducing a scheme for the lim-
iting macroscopic equation in the limit of infinite time-scale separation. For de-
terministic discretisation methods, there is a long line of such methods. We refer
to [
2
,
5
,
6
,
10
,
14
,
19
,
20
,
21
,
22
,
23
,
24
,
26
,
27
,
28
] as a representative sample
of such methods in the diffusive scaling. The recent review paper [
15
] contains
an overview of the state of the art on asymptotic-preserving methods for kinetic
equations, and ample additional references. In the particle-based setting, only a few
asymptotic-preserving methods have been developed, mostly in the hyperbolic scal-
ing [
11
,
12
,
13
,
29
,
30
,
31
]. In the diffusive scaling, there are only of two works [
9
,
16
]
so far, to the best of our knowledge. Both methods avoid the time step restrictions
Page: 2 job: MLMCforAP macro: svmult.cls date/time: 15-Jan-2019/18:31

MLMC AP Particle Method for Kinetic Equations in the Diffusion Limit 3
caused by fast problem time-scales, at the expense of introducing a bias, which is of
order one in the time step size.
The goal of the present paper is to combine the asymptotic-preserving scheme
in [
16
] with the multilevel Monte Carlo method. Given a fixed computational budget,
a trade-off typically has to be made between a small bias and a low variance. The
former can be obtained by reducing the time step, the latter by simulating many
trajectories with large time steps. The core idea behind the multilevel Monte Carlo
method [
17
] is to reduce computational cost, by combining estimates computed
with different time step sizes. The multilevel Monte Carlo method was originally
developed in the context of stochastic processes, and has been applied to problems
across many fields, for example, finance [
17
] and data science [
1
]. The method has
also successfully been applied to simulating large PDE’s with random coefficients [
8
].
Recent work has also used multilevel Monte Carlo methods in an optimization
context [34].
The remainder of this paper is organized as follows. In Section 2, we describe the
model kinetic equation on which we will demonstrate our approach, as well as the
asymptotic-preserving Monte Carlo scheme that was introduced in [
16
]. In Section 3,
we cover the multilevel Monte Carlo method that is the core contribution of this paper.
In Section 4, we present some preliminary experimental results, demonstrating the
properties of the new scheme as well as its computational gain. Finally, in Section 5
we will summarize our main results and mention some possible future extensions.
2 Model problem and asymptotic-preserving scheme
2.1 Model equation in the diffusive limit
The model problem considered in this work is a one-dimensional kinetic equation in
the diffusive scaling of the form (5), which we rewrite as
t
f (x,v,t) +
v
ε
x
f (x,v,t) =
1
ε
2
(ρ(x,t) f (x, v,t)). (5)
For ease of exposition, we restrict ourselves to the case of two discrete velocities,
v =
±1
. Then, we can write
f
+
(x,t)
and
f
(x,t)
to represent the distribution of particles
with respective positive and negative velocities, and
ρ(x,t) = f
+
(x,t) + f
(x,t)
represents the total density of particles. In this case, equation (5) simplifies to
t
f
+
(x,t) +
1
ε
x
f
+
(x,t) =
1
ε
2
ρ(x,t)
2
f
+
(x,t)
t
f
(x,t)
1
ε
x
f
(x,t) =
1
ε
2
ρ(x,t)
2
f
(x,t)
. (6)
Page: 3 job: MLMCforAP macro: svmult.cls date/time: 15-Jan-2019/18:31

4 Emil Løvbak, Giovanni Samaey, and Stefan Vandewalle
Equation
(6)
is also known as the Goldstein-Taylor model, and can be solved using a
particle scheme. We then introduce a time step t and an ensemble of P particles

X
n
p,t
,V
n
p,t

P
p=1
. (7)
The particle state (position and velocity) is represented as
(X,V )
,
p
is the particle
index (
1 p P
), and
n
represents the time index, i.e.,
X
n
p,t
X
p
(nt)
. Equation
(6)
is then solved via operator splitting as
1. Transport step. The position of each particle is updated based on its velocity
X
n+1
p,t
= X
n
p,t
±V
n
p,t
t. (8)
2. Collision step. During collisions, each particle’s velocity is updated as:
V
n+1
p,t
=
(
±1/ε, with probability p
c
t
= t/ε
2
and equal probability in the sign,
V
n
p,t
, otherwise.
(9)
To maintain stability, this approximation requires a time step restriction
t = O (ε
2
)
,
leading to unacceptably high computational costs.
2.2 Asymptotic-preserving Monte Carlo scheme
Recently, an asymptotic-preserving Monte Carlo scheme was proposed [
16
], based
on the simulation of a modified equation
t
f
+
+
ε
ε
2
+ t
x
f
+
=
t
ε
2
+ t
xx
f
+
+
1
ε
2
+ t
ρ
2
f
+
t
f
ε
ε
2
+ t
x
f
=
t
ε
2
+ t
xx
f
+
1
ε
2
+ t
ρ
2
f
. (10)
In
(10)
we have dropped the space and time dependency of
f
±
and
ρ
, for conciseness.
The model given by
(10)
reduces to
(6)
in the limit when
t
tends to zero, and has
an O(t) bias. In the limit when ε tends to zero, the equations reduce to (4).
Discretizing this equation, using operator splitting as above, again leads to a
Monte Carlo scheme. For each particle
X
p
and for each time step
n
, one time step
now consists of a transport-diffusion and a collision step:
1. Transport-diffusion step.
The position of the particle is updated based on its
velocity and a Brownian increment
X
n+1
p,t
= X
n
p,t
±
ε
ε
2
+ t
t +
2t
r
t
ε
2
+ t
ξ
n
p
= X
n
p,t
+V
n
p,t
t +
2t
p
D
t
ξ
n
p
, (11)
Page: 4 job: MLMCforAP macro: svmult.cls date/time: 15-Jan-2019/18:31

MLMC AP Particle Method for Kinetic Equations in the Diffusion Limit 5
in which we have taken
ξ
n
p
N (0,1)
and introduced a
t
-dependent velocity
V
n
p,t
and diffusion coefficient D
t
:
V
n
p,t
=
ε
ε
2
+ t
, D
t
=
t
ε
2
+ t
. (12)
2. Collision step. During collisions, each particle’s velocity is updated as:
V
n+1
p,t
=
±
ε
ε
2
+ t
, with probability p
c
t
=
t
ε
2
+ t
and equal probability in the sign,
V
n
p,t
, otherwise.
(13)
For more details, we refer the reader to [16].
3 Multilevel Monte Carlo method
3.1 Method and notation
We want to estimate some quantity of interest
Y
that is a function of the particle
distribution
f (x,v,t)
at some specific moment
t = t
in time, i.e., we are interested in
Y (t
) = E[F(X(t
))] =
Z Z
F(x) f (x,v,t
)dxdv. (14)
Note that, in equation
(14)
, the function
F
only depends on the position
x
and not on
velocity. That is a choice we make for notational convenience, and is not essential for
the method we present. The classical Monte Carlo estimator
ˆ
Y (t
)
for
(14)
is given
by
ˆ
Y (t
) =
1
P
P
p=1
F(X
N
p,t
), t
= Nt. (15)
Here,
P
denotes the number of simulated trajectories,
N
the number of simulated time
steps,
t
the time step size and
X
N
p,t
is generated by the time-discretised process
(11)
(13)
. Given a constrained computational budget, a trade-off has to be made
when selecting the time step size
t
. On the one hand, a small time step reduces the
bias of the simulation of each sampled trajectory, and thus of the estimated quantity
of interest. On the other hand, a large time step reduces the cost per trajectory,
which increases the number of trajectories that can be simulated and thus reduces
the resulting variance on the estimate. The key idea behind the Multilevel Monte
Carlo method [
17
] is to generate a sequence of estimates with varying discretization
accuracy and a varying number of realizations. The method achieves the bias of the
finest discretization, with the variance of the coarsest discretization.
Page: 5 job: MLMCforAP macro: svmult.cls date/time: 15-Jan-2019/18:31

Citations
More filters

Journal ArticleDOI
Abstract: Kinetic equations model distributions of particles in position-velocity phase space. Often, one is interested in studying the long-time behavior of particles in high-collisional regimes in which an approximate (advection)-diffusion model holds. In this paper we consider the diffusive scaling. Classical particle-based techniques suffer from a strict time-step restriction in this limit, to maintain stability. Asymptotic-preserving schemes avoid this problem, but introduce an additional time discretization error, possibly resulting in an unacceptably large bias for larger time steps. Here, we present and analyze a multilevel Monte Carlo scheme that reduces this bias by combining estimates using a hierarchy of different time step sizes. We demonstrate how to correlate trajectories from this scheme, using different time steps. We also present a strategy for selecting the levels in the multilevel scheme. Our approach significantly reduces the computation required to perform accurate simulations of the considered kinetic equations, compared to classical Monte Carlo approaches.

References
More filters

Journal ArticleDOI
Abstract: A kinetic theory approach to collision processes in ionized and neutral gases is presented. This approach is adequate for the unified treatment of the dynamic properties of gases over a continuous range of pressures from the Knudsen limit to the high-pressure limit where the aerodynamic equations are valid. It is also possible to satisfy the correct microscopic boundary conditions. The method consists in altering the collision terms in the Boltzmann equation. The modified collision terms are constructed so that each collision conserves particle number, momentum, and energy; other characteristics such as persistence of velocities and angular dependence may be included. The present article illustrates the technique for a simple model involving the assumption of a collision time independent of velocity; this model is applied to the study of small amplitude oscillations of one-component ionized and neutral gases. The initial value problem for unbounded space is solved by performing a Fourier transformation on the space variables and a Laplace transformation on the time variable. For uncharged gases there results the correct adiabatic limiting law for sound-wave propagation at high pressures and, in addition, one obtains a theory of absorption and dispersion of sound for arbitrary pressures. For ionized gases the difference in the nature of the organization in the low-pressure plasma oscillations and in high-pressure sound-type oscillations is studied. Two important cases are distinguished. If the wavelengths of the oscillations are long compared to either the Debye length or the mean free path, a small change in frequency is obtained as the collision frequency varies from zero to infinity. The accompanying absorption is small; it reaches its maximum value when the collision frequency equals the plasma frequency. The second case refers to waves shorter than both the Debye length and the mean free path; these waves are characterized by a very heavy absorption.

5,976 citations


BookDOI
01 Jan 1991
Abstract: PART 1: PRIMER Why attempting to do plasma physics via computer simulation using particles makes good sense Overall view of a one dimensional electrostatic program A one dimensional electrostatic program ES1 Introduction to the numerical methods used Projects for ES1 A 1d electromagnetic program EM1 Projects for EM1 PART 2: THEORY Effects of the spatial grid Effects of the finitw time ste Energy-conserving simulation models Multipole models Kinetic theory for fluctuations and noise collisions Kinetic properties: theory, experience and heuristic estimates PART 3: PRACTICE Electrostatic programs in two and three dimensions Electromagnetic programs in two and three dimensions Particle loading, injection boudary conditions and external circuit PART 4: APPENDICES

4,558 citations


Book
01 Jan 1988
Abstract: I. Basic Principles of The Kinetic Theory of Gases.- 1. Introduction.- 2. Probability.- 3. Phase space and Liouville's theorem.- 4. Hard spheres and rigid walls. Mean free path.- 5. Scattering of a volume element in phase space.- 6. Time averages, ergodic hypothesis and equilibrium states.- References.- II. The Boltzmann Equation.- 1. The problem of nonequilibrium states.- 2. Equations for the many particle distribution functions for a gas of rigid spheres.- 3. The Boltzmann equation for rigid spheres.- 4. Generalizations.- 5. Details of the collision term.- 6. Elementary properties of the collision operator. Collision invariants.- 7. Solution of the equation Q(f,f) = 0.- 8. Connection between the microscopic description and the macroscopic description of gas dynamics.- 9. Non-cutoff potentials and grazing collisions. Fokker-Planck equation.- 10. Model equations.- References.- III. Gas-Surface Interaction and the H-Theorem.- 1. Boundary conditions and the gas-surface interaction.- 2. Computation of scattering kernels.- 3. Reciprocity.- 4. A remarkable inequality.- 5. Maxwell's boundary conditions. Accommodation coefficients.- 6. Mathematical models for gas-surface interaction.- 7. Physical models for gas-surface interaction.- 8. Scattering of molecular beams.- 9. The H-theorem. Irreversibility.- 10. Equilibrium states and Maxwellian distributions.- References.- IV, Linear Transport.- 1. The linearized collision operator.- 2. The linearized Boltzmann equation.- 3. The linear Boltzmann equation. Neutron transport and radiative transfer.- 4. Uniqueness of the solution for initial and boundary value problems.- 5. Further investigation of the linearized collision term.- 6. The decay to equilibrium and the spectrum of the collision operator.- 7. Steady one-dimensional problems. Transport coefficients.- 8. The general case.- 9. Linearized kinetic models.- 10. The variational principle.- 11. Green's function.- 12. The integral equation approach.- References.- V. Small and Large Mean Free Paths.- 1. The Knudsen number.- 2. The Hilbert expansion.- 3. The Chapman-Enskog expansion.- 4. Criticism of the Chapman-Enskog method.- 5. Initial, boundary and shock layers.- 6. Further remarks on the Chapman-Enskog method and the computation of transport coefficients.- 7. Free molecule flow past a convex body.- 8. Free molecule flow in presence of nonconvex boundaries.- 9. Nearly free-molecule flows.- References.- VI. Analytical Solutions of Models.- 1. The method of elementary solutions.- 2. Splitting of a one-dimensional model equation.- 3. Elementary solutions of the simplest transport equation.- 4. Application of the general method to the Kramers and Milne problems.- 5. Application to the flow between parallel plates and the critical problem of a slab.- 6. Unsteady solutions of kinetic models with constant collision frequency.- 7. Analytical solutions of specific problems.- 8. More general models.- 9. Some special cases.- 10. Unsteady solutions of kinetic models with velocity dependent collision frequency.- 11. Analytic continuation.- 12. Sound propagation in monatomic gases.- 13. Two-dimensional and three-dimensional problems. Flow past solid bodies.- 14. Fluctuations and light scattering.- References.- VII. The Transition Regime.- 1. Introduction.- 2. Moment and discrete ordinate methods.- 3. The variational method.- 4. Monte Carlo methods.- 5. Problems of flow and heat transfer in regions bounded by planes or cylinders.- 6. Shock-wave structure.- 7. External flows.- 8. Expansion of a gas into a vacuum.- References.- VIII. Theorems on the Solutions of the Boltzmann Equation.- 1. Introduction.- 2. The space homogeneous case.- 3. Mollified and other modified versions of the Boltzmann equation.- 4. Nonstandard analysis approach to the Boltzmann equation.- 5. Local existence and validity of the Boltzmann equation.- 6. Global existence near equilibrium.- 7. Perturbations of vacuum.- 8. Homoenergetic solutions.- 9. Boundary value problems. The linearized and weakly nonlinear cases.- 10. Nonlinear boundary value problems.- 11. Concluding remarks.- References.- References.- Author Index.

2,797 citations


Journal ArticleDOI
TL;DR: It is shown that multigrid ideas can be used to reduce the computational complexity of estimating an expected value arising from a stochastic differential equation using Monte Carlo path simulations.
Abstract: We show that multigrid ideas can be used to reduce the computational complexity of estimating an expected value arising from a stochastic differential equation using Monte Carlo path simulations. In the simplest case of a Lipschitz payoff and a Euler discretisation, the computational cost to achieve an accuracy of O(e) is reduced from O(e-3) to O(e-2 (log e)2). The analysis is supported by numerical results showing significant computational savings.

1,393 citations


Journal ArticleDOI
TL;DR: A novel variance reduction technique for the standard Monte Carlo method, called the multilevel Monte Carlo Method, is described, and numerically its superiority is demonstrated.
Abstract: We consider the numerical solution of elliptic partial differential equations with random coefficients Such problems arise, for example, in uncertainty quantification for groundwater flow We describe a novel variance reduction technique for the standard Monte Carlo method, called the multilevel Monte Carlo method, and demonstrate numerically its superiority The asymptotic cost of solving the stochastic problem with the multilevel method is always significantly lower than that of the standard method and grows only proportionally to the cost of solving the deterministic problem in certain circumstances Numerical calculations demonstrating the effectiveness of the method for one- and two-dimensional model problems arising in groundwater flow are presented

496 citations