remarkRemark hypothesisHypothesis claimClaim On MLMC for hyperbolic
equationsJunpeng Hu, Shi Jin, Jinglai Li and Lei Zhang

# On multilevel Monte Carlo methods for deterministic and uncertain hyperbolic
systems††thanks: Submitted to the editors DATE. SJ was supported by National
Key R&D Program of China Project No. 2020YFA0712000, the Strategic Priority
Research Program of Chinese Academy of Sciences, XDA25010404 and NSFC grant
No. 12031013. LZ’s work was partially supported by NSFC grants No. 11871339
and 11861131004.

Junpeng Hu School of Mathematical Sciences, Institute of Natural Sciences, and
MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, China (, , ).
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Shi
Jin22footnotemark: 2 Jinglai Li School of Mathematics, University of
Birmingham, Birmingham, B15 2TT, United Kingdom<EMAIL_ADDRESS>Lei
Zhang22footnotemark: 2

###### Abstract

In this paper, we evaluate the performance of the multilevel Monte Carlo
method (MLMC) for deterministic and uncertain hyperbolic systems, where
randomness is introduced either in the modeling parameters or in the
approximation algorithms. MLMC is a well known variance reduction method
widely used to accelerate Monte Carlo (MC) sampling. However, we demonstrate
in this paper that for hyperbolic systems, whether MLMC can achieve a real
boost turns out to be delicate. The computational costs of MLMC and MC depend
on the interplay among the accuracy (bias) and the computational cost of the
numerical method for a single sample, as well as the variances of the sampled
MLMC corrections or MC solutions. We characterize three regimes for the MLMC
and MC performances using those parameters, and show that MLMC may not
accelerate MC and can even have a higher cost when the variances of MC
solutions and MLMC corrections are of the same order. Our studies are carried
out by a few prototype hyperbolic systems: a linear scalar equation, the Euler
and shallow water equations, and a linear relaxation model, the above
statements are proved analytically in some cases, and demonstrated numerically
for the cases of the stochastic hyperbolic equations driven by white noise
parameters and Glimm’s random choice method for deterministic hyperbolic
equations.

###### keywords:

uncertainty quantification, multilevel Monte Carlo method, hyperbolic
equations, variance reduction

35R60, 35Q20, 65C05

##  1 Introduction

The Monte Carlo (MC) method is a widely used computational tool from
multidimensional integration, to statistical physics and finance. It is a
sampling method which generates a large number of samples to retrieve
statistical information that can be a good approximation of a given random or
deterministic problem, without curse-of-dimensionality. The accuracy of MC is
governed by the variance of those samples and the bias induced by the
discretization method. To further reduce the variance and to accelerate the
sampling algorithm, Multilevel Monte Carlo (MLMC) method [12] can be adopted
by decomposing the calculation of the quantity of interest (QoI) to
hierarchical levels. MLMC has achieved great success for the efficient
simulation of stochastic differential equations (SDE) [1, 9, 12], elliptic
PDEs with random coefficients [5, 8], and particle methods for transport
equations [24].

In this paper, we consider linear and nonlinear hyperbolic systems where MC
and MLMC can be applied in the following circumstances: i), for hyperbolic
equations with random parameters; ii), for randomized algorithms of
deterministic hyperbolic equations. These two cases have interesting
applications in random algorithms for deterministic problems [6, 13], or
problems in uncertainty quantification [2, 3, 16, 27]. Although MLMC has be
shown to accelerate MC for many stochastic PDEs and randomized algorithms, we
will demonstrate that, at least in the context of hyperbolic systems, the
situation is more subtle.

Let us start by considering the stochastic differential equation ${\rm
d}u(t)=a(u,t){\rm d}t+b(u,t){\rm d}W$. To compute certain QoI, for example,
$\mathrm{E}\left[u(T)\right]$, one can employ a discretization method such as
the explicit Euler method, $u_{n+1}=u_{n}+a(u_{n},t_{n})\Delta
t+b(u_{n},t_{n})\Delta W$, with time step $\Delta t$. Let $N$ be the number of
samples, and $\sigma$ be the standard deviation of each sample. It is well
known that, the mean square error ($\mathrm{MSE}$) is bounded by
$c_{1}N^{-1}\sigma+c_{2}\Delta t^{2}$, where $c_{1}$, $c_{2}$ are positive
constants. The first term comes from the statistical error of the MC
estimation, and the second term comes from the bias, which is in turn bounded
by the discretization error of the numerical method. For SDE, $\sigma$ is a
constant independent of $\Delta t$, and we can take $\Delta t\approx
O(\delta)$, $N\approx O(\delta^{-2})$, such that the $\mathrm{MSE}$ is
$O(\delta^{2})$ and the computational cost is $O(\delta^{-3})$.

The multilevel Monte Carlo (MLMC) method proposed by Giles [12] et. al. is one
of the state-of-the-art variance reduction techniques to further reduce the
computational cost. MLMC samples the solution $P_{0}^{(0)}$ on the coarsest
level and corrections $P^{(l)}_{l}-P^{(l-1)}_{l}$ on the finer levels, then
assembles an estimator for the QoI on the finest mesh by a telescoping sum. A
sample correction $P^{(l)}_{l}-P^{(l-1)}_{l}$ can be obtained using solutions
$P^{(l)}_{l}$, $P^{(l-1)}_{l}$ with the same set of random variables (random
trajectory) on the $l$th level, but with different time steps $\Delta t_{l}$
and $\Delta t_{l-1}$. The key to MLMC is the reduction of
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ without changing the
expectation $\mathrm{E}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$. It can be
proved that for SDE, the computational cost of MLMC can be reduced to
$O(\delta^{-2}(\log\delta)^{2})$ [12, Theorem 3.1], under certain conditions
for the variances of corrections.

For hyperbolic equations, as one needs to consider both spatial and temporal
variables, the situation is more complex. It turns out that, whether or not
MLMC can enhance the performance of a MC sampling algorithm for hyperbolic
problems, depends on the interplay between the bias induced by the numerical
method for a single sample, variances of the sampled MC solutions or MLMC
corrections, and the computational cost of the numerical method for a single
sample, which can be characterized by problem dependent parameters.

To be more specific, let $P_{l}^{(l)}$ be the MC solution at the $l$th level,
and $Y_{l}:=P^{(l)}_{l}-P^{(l-1)}_{l}$ be the MLMC correction at the $l$th
level ($l\geq 1$). Let

  * •

$\alpha$ denote the order of accuracy for the numerical method of a single
sample (bias). In this paper, we always use first order method with $\alpha=1$
unless otherwise stated.

  * •

$\beta_{0}$ and $\beta$ denote the orders of the variances of MC solutions and
MLMC corrections, respectively. Namely
$\mathsf{V}_{l}:=\mathrm{Var}\left[P_{l}^{(l)}\right]\leq c_{2}(\Delta
t_{l})^{\beta_{0}}$, and $V_{l}:=\mathrm{Var}\left[Y_{l}\right]\leq
c_{2}(\Delta t_{l})^{\beta}$.

  * •

$\gamma$ denote the order of the computational cost for a single sample such
that the cost is $O(\Delta t^{-\gamma})$. For explicit numerical methods,
$\gamma$ corresponds to the physical dimension of the problem, for example,
$\gamma=1$ for SDE, and $\gamma=2$ for 1D SPDE.

We can estimate the cost of MLMC and MC algorithms through parameters
$\alpha$, $\beta$ or $\beta_{0}$, $\gamma$, and the desired accuracy
$\delta>0$, using Eq. 17 and Eq. 19 in § 2. In general, when
$\beta_{0}\leq\beta\leq\gamma$, up to a logarithmic factor, the cost of MC is
$O(\delta^{-2-(\gamma-\beta_{0})/\alpha})$, and the cost of MLMC is
$O(\delta^{-2-(\gamma-\beta)/\alpha})$. Therefore, when $\beta_{0}=\beta$, the
computational costs of MC and MLMC are of the same order. Namely, MLMC does
not accelerate and may even cost more than MC. In fact, we have identified a
few such examples, such as hyperbolic equations with white noise parameters,
as well as the random choice algorithm for the (deterministic) Jin-Xin model
[18] which is a hyperbolic system with relaxation.

We list the results in Table 1 for the problems from literature [1, 5, 8, 9,
12, 24] and the hyperbolic problems that we study in this paper.

  * •

Case 1: stochastic differential equation (SDE) [1, 9, 12];

  * •

Case 2: random elliptic equation [5, 8];

  * •

Case 3: particle method for transport equations [24];

  * •

Case 4: scalar advection equation with random time dependent velocity
$a(t,\omega)$, see § 3,

    * –

4.1, $a$ is piecewise constant in finite number of time intervals,

    * –

4.2, $a$ is uniform white noise in time;

  * •

Case 5: Euler equations with random time dependent adiabatic constant
$\lambda(t,\omega)$, see § 4.1,

    * –

5.1, $\lambda$ is piecewise constant in finite number of time intervals,

    * –

5.2, $\lambda$ is uniform white noise in time;

  * •

Case 6: shallow water equation with random bottom topography $B(x,\omega)$,
see § 4.2,

    * –

6.1, $B$ has a finite number of random parameters,

    * –

6.2, $B$ has uniform white noise parameters;

  * •

Case 7: random choice method for deterministic Jin-Xin model, see § 5,

    * –

7.1, semi-random case, only use random choice for the relaxation step,

    * –

7.2, fully-random case, use random choice for both the convection and the
relaxation steps.

Cases | Dimension | Parameters | cost | Regime
---|---|---|---|---
physical | random ($K$) | $\alpha$ | $\beta_{0}$ | $\beta$ | $\gamma$ | MC | MLMC
1 | $t$ |  $T/\Delta t$, in $t$ | 1 | 0 | 1 | 1 | $\delta^{-3}$ | $\delta^{-2}(\log\delta^{-1})^{2}$ | I
2 | $x$ |  $O(1)$, in $x$ | $1.75^{*}$ | $0^{*}$ | $2^{*}$ | 1 | $\delta^{-2.57}$ | $\delta^{-2}$ | I
3 | $x,t$ |  $T/\Delta t$, in $t$ | 1 | $0^{*}$ | 1 | 1 | $\delta^{-3}$ | $\delta^{-2}(\log\delta^{-1})^{2}$ | I
4 | 4.1 | $x,t$ |  $O(1)$, in $t$ | 1 | 0 | 2 | 2 | $\delta^{-4}$ | $\delta^{-2}(\log\delta^{-1})^{2}$ | I
4.2 | $x,t$ |  $T/\Delta t$, in $t$ | 1 | 1 | 1 | 2 | $\delta^{-3}$ | $\delta^{-3}$ | II
5 | 5.1 | $x,t$ |  $O(1)$, in $t$ | 1 | $0^{*}$ | $1.5^{*}$ | 2 | $\delta^{-4}$ | $\delta^{-2.5}$ | I
5.2 | $x,t$ |  $T/\Delta t$, in $t$ | 1 | $1^{*}$ | $1^{*}$ | 2 | $\delta^{-3}$ | $\delta^{-3}$ | II
6 | 6.1 | $x,t$ |  $O(1)$, in $x$ | 1 | $0^{*}$ | $1^{*}$ | 2 | $\delta^{-4}$ | $\delta^{-3}$ | I
6.2 | $x,t$ |  $T/\Delta t$, in $x$ | 1 | $0^{*}$ | $0^{*}$ | 2 | $\delta^{-4}$ | $\delta^{-4}$ | III
7 | 7.1 | $x,t$ |  $\left(T/\Delta t\right)^{2}$, in $x,t$ | 1 | $0.5^{*}$ | $0.5^{*}$ | 2 | $\delta^{-3.5}$ | $\delta^{-3.5}$ | II
7.2 | $x,t$ |  $\left(T/\Delta t\right)^{2}$, in $x,t$ | 1 | $0^{*}$ | $0^{*}$ | 2 | $\delta^{-4}$ | $\delta^{-4}$ | III
Table 1: Comparison of MLML and MC. Parameters $\alpha$, $\beta_{0}$, $\beta$,
and $\gamma$ are either analytically derived or obtained by numerical
experiments which are marked with a $*$. The complexities are computed using
Eq. 17 and Eq. 19 in § 2, and numerically validated in literature or in the
numerical experiments in § 3, § 4 and § 5. The expressions ‘in $x$’, ‘in $t$’,
‘in $x,t$’, in the ‘random ($K$)’ column indicate the dependence on the random
parameters, $K$ is the number of random parameters (random dimension) in the
problem, $T$ is the computational time, and $\Delta t$ is the time step of the
numerical method.

From above results in Table 1, we can identify three regimes:

  * •

Regime I, $\beta>\beta_{0}=0$, for example, SDE and SPDEs with finite number
of random parameters, where MLMC outperforms MC by orders of magnitude.

  * •

Regime II, $\beta_{0}=\beta>0$, where we can observe the decay of variances of
MC solutions and MLMC corrections. The costs of both MLMC and MC are of the
order $O(\delta^{-4+\beta})$, which is the case for the scalar advection
equation and the Euler equations with uniform white noise parameters in time,
and the semi-random choice method for the Jin-Xin model.

  * •

Regime III, $\beta_{0}=\beta\simeq 0$, namely, the variances of MC solutions
and MLMC corrections do not decay with respect to discretization parameters
(time step and mesh size). The costs of both MLMC and MC are of the worst
possible order $\delta^{-4}$, which is the case for the shallow water equation
with uniform white noise parametrized topography and the fully-random choice
method for the Jin-Xin model.

#### Outline

The rest of this paper is organized as follows. In § 2, we describe the
principle and algorithm of MLMC, and present the main theorems – the cost
analysis of MLMC and standard MC methods. In § 3, we first carry out an in-
depth analysis for the scalar advection equation with a random velocity, then
justify it with numerical experiments. In § 4 and § 5, we provide
comprehensive numerical experiments for other test examples, including
hyperbolic equations with random parameters such as the Euler equations and
the shallow water equations, and the random choice method for the multiscale
Jin-Xin model. We conclude the paper in Section § 6.

#### Notation

The symbol $C$ (or $c$) denotes generic positive constants that may change
from one line of an estimate to the next. When estimating rates of decay or
convergence, $C$ will always remain independent of approximation parameters.
The dependence of $C$ will be clear from the context or stated explicitly. The
expectation and variance of a random variable $X$ is abbreviated as
$\mathrm{E}\left[X\right]$ and $\mathrm{Var}\left[X\right]$ respectively. The
covariance of two random variables $X$, $Y$ is abbreviated as
$\mathrm{Cov}\left[X,Y\right]$. To further simplify notation we will often
write $\approx$ to mean that both sides of $\approx$ are equivalent
infinitesimal.

##  2 The Multilevel Monte Carlo method

In this section, we first introduce the Multilevel Monte Carlo (MLMC) method
and formulate the MC and MLMC algorithms in § 2.1. Then we present the
theorems for their computational costs in § 2.2, which provide the crucial
theoretical foundation of the paper.

###  2.1 Formulation and algorithm

####  2.1.1 MLMC formulation

Given the time dependent hyperbolic problem with randomness, we denote $\Delta
t_{l}$ as the time step on hierarchical levels $l=0,1,\dots,L$, with $\Delta
t_{l-1}=\Gamma\Delta t_{l}$ with the integer refinement ratio $\Gamma>1$, and
$\Delta x_{l}$ as the corresponding grid size which satisfies the CFL
stability condition.

For simplicity, we denote $P^{(l)}$ as the random numerical solution at final
time $T$ computed using time step $\Delta t_{l}$. The QoI
$\mathrm{E}\left[P^{(L)}\right]$ for MLMC can be expanded as a telescoping sum

(1)
$\mathrm{E}\left[P^{(L)}\right]=\mathrm{E}\left[P^{(0)}\right]+\sum_{l=1}^{L}\left(\mathrm{E}\left[P^{(l)}\right]-\mathrm{E}\left[P^{(l-1)}\right]\right).$

Instead of approximating $\mathrm{E}\left[P^{(L)}\right]$ on the finest level
$L$, we can approximate the terms on the right hand side of (1) by independent
samples on hierarchical levels $l=0,1,\dots,L$, and sum them together. The
name of the game is now how to reduce the sample variances.

To be more precise, let $P_{l}^{(l^{\prime})}(s)$ be the $s$th sampled
solution at final time $T$ computed using random variables on the $l$th level,
as well as the time step $\Delta t_{l^{\prime}}$ and the grid size $\Delta
x_{l^{\prime}}$ on the level $l^{\prime}$, with $l^{\prime}=l-1,l$. Let
$N_{l}$ be the number of samples on the $l$th level. That is to say, we have
samples $P_{0}^{(0)}(s)$ on the $0$th level, with $s=1,\dots,N_{0}$, and
$P_{l}^{(l)}(s)$, $P_{l}^{(l-1)}(s)$ on the $l$th level, with
$s=1,\dots,N_{l}$, $l=1,\dots,L$.

We call $P_{l}^{(l)}-P_{l}^{(l-1)}$ as the MLMC correction. Since
$P_{l}^{(l)}(s)$ and $P_{l}^{(l-1)}(s)$ are solutions on successive levels
using the same random variables on the $l$th level, we have
$\mathrm{E}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=\mathrm{E}\left[P^{(l)}\right]-\mathrm{E}\left[P^{(l-1)}\right]$.
At the same time, we hope that
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ is small.

Hence, we have the following MLMC estimate

(2) $\displaystyle\hat{P}^{(L)}$
$\displaystyle:=N_{0}^{-1}\sum_{s=1}^{N_{0}}P_{0}^{(0)}(s)+\sum_{l=1}^{L}N_{l}^{-1}\sum_{s=1}^{N_{l}}\left(P_{l}^{(l)}(s)-P_{l}^{(l-1)}(s)\right)$
$\displaystyle\triangleq\hat{P}_{0}^{(0)}+\sum_{l=1}^{L}\left(\hat{P}_{l}^{(l)}-\hat{P}_{l}^{(l-1)}\right).$

####  2.1.2 Optimal parameters for MLMC

In the MLMC estimate (2), one needs to choose the number of samples $N_{l}$ on
hierarchical levels. The optimal value of $N_{l}$ can be determined in terms
of the cost and variance of MLMC corrections by the following argument.

Let $C_{0}$, $V_{0}$ be the computational cost and variance of $P_{0}^{(0)}$,
and $C_{l}$, $V_{l}$ be the computational cost and variance of
$P_{l}^{(l)}-P_{l}^{(l-1)}$, then the computational cost $\mathfrak{C}$ and
the variance of $\hat{P}^{(L)}$ are

(3)
$\mathfrak{C}=\sum_{l=0}^{L}N_{l}C_{l},\quad\mathrm{Var}\left[\hat{P}^{(L)}\right]=\sum_{l=0}^{L}N_{l}^{-1}V_{l}.$

Let $R=\mathrm{Var}\left[\hat{P}^{(L)}\right]$ be the variance of
$\hat{P}^{(L)}$, and our aim is to minimize the cost with respect to $N_{l}$.
Using the Lagrange function $L$ with the multiplier $\mu^{2}$, we have

(4)
$L(N_{l},\mu^{2})=\sum_{l=0}^{L}N_{l}C_{l}+\mu^{2}\left(\sum_{l=0}^{L}N_{l}^{-1}V_{l}-R\right).$

The Kuhn-Tucker conditions (KKT) condition writes

(5) $\frac{\partial L(N_{l},\mu^{2})}{\partial
N_{l}}=C_{l}-\mu^{2}V_{l}N_{l}^{-2}=0,\quad\Rightarrow
N_{l}=\mu\sqrt{V_{l}/C_{l}},$

which means that the optimal $N_{l}$ is proportional to $\sqrt{V_{l}/C_{l}}$.
By substituting (5) into $R=\mathrm{Var}\left[\hat{P}^{(L)}\right]$, one gets

(6) $\mu=R^{-1}\sum_{l=0}^{L}\sqrt{V_{l}C_{l}},\quad
N_{l}=R^{-1}\sqrt{V_{l}/C_{l}}\left(\sum_{l=0}^{L}\sqrt{V_{l}C_{l}}\right).$

######  Remark 2.1.

The variance $V_{l}$ is not known a priori, and one needs to estimate these
values through samples in the implementation.

####  2.1.3 An adaptive algorithm

In this paper, we adopt the MLMC algorithm (as well as MC) with adaptively
refined discretization levels as formulated in Giles’ paper [12]. We decompose
the $\mathrm{MSE}$ into bias and variance, and require that both the squared
bias and the variance of the estimator are bounded by $\delta^{2}/2$. This
requirement leads to the adaptive criteria for the finest level $L$. The
difference of MLMC and MC is that MLMC combines information from all the
levels, while MC only uses the samples on the finest level and the coarser
levels are only used to determine the finest level.

For simplicity, we denote $P$ as the exact solution at the final time $T$,
$\hat{P}_{l}^{(l)}$, $\hat{P}_{l}^{(l-1)}$ as the estimated MLMC solutions on
level $l$ with time step $\Delta t_{l}$, $\Delta t_{l-1}=\Gamma\Delta t_{l}$
respectively, and

(7) $\hat{Y}_{l}=\left\\{\begin{aligned} &\hat{P}_{0}^{(0)},&l=0,\\\
&\hat{P}_{l}^{(l)}-\hat{P}_{l}^{(l-1)},&l>0,\end{aligned}\right.$

as estimated MLMC corrections.

Given the accuracy threshold $\delta>0$ for the MLMC algorithm, we start from
the initial $0$th level. Using an initial $N_{0}=N_{I}$ number of samples, we
obtain the Monte Carlo estimates $\hat{V}_{0}$, $\hat{Y}_{0}$ of $V_{0}$,
$Y_{0}$ respectively. Then the updated number of samples on $0$th level
$N_{0}^{\prime}$ can be determined by $N_{0}^{\prime}=\left\lceil
2\delta^{-2}\hat{V}_{0}\right\rceil$. If $N_{0}^{\prime}>N_{0}$, we evaluate
extra samples and update $\hat{V}_{0}$, $\hat{Y}_{0}$, $N_{0}=N_{0}^{\prime}$.

For the iteration $l\rightarrow l+1$, $V_{l+1}$, $Y_{l+1}$ can be estimated by
$\hat{V}_{l+1}$, $\hat{Y}_{l+1}$ at the new level $l+1$ using $N_{l+1}=N_{I}$
samples, then the optimal number of samples on the $k$th level
$N_{k}^{\prime}$, with $k=0,\dots,l+1$, can be determined by

(8) $N_{k}^{\prime}=\left\lceil
2\delta^{-2}\sqrt{\hat{V}_{k}/C_{k}}\left(\sum_{j=0}^{l+1}\sqrt{\hat{V}_{j}C_{j}}\right)\right\rceil,$

where $C_{j}$ is the computational cost on level $j$. If
$N_{k}^{\prime}>N_{k}$, we let $N_{k}=N_{k}^{\prime}$, evaluate extra samples
and update $\hat{V}_{k}$, $\hat{Y}_{k}$.

We can identify the finest level $L$ by an a posterior error estimate through
the requirement that the squared bias of the estimator is bounded by
$\delta^{2}/2$. For a numerical scheme with (weak) convergence rate $\alpha$,
we have

(9)
$\mathrm{E}\left[\hat{P}_{l}^{(l)}-P\right]=\mathrm{E}\left[P_{l}^{(l)}-P\right]\approx
c_{1}(\Delta t_{l})^{\alpha}.$

The following estimate of the correction
$\mathrm{E}\left[\hat{P}_{l}^{(l)}-\hat{P}_{l}^{(l-1)}\right]$

(10)
$\mathrm{E}\left[\hat{P}_{l}^{(l)}-\hat{P}_{l}^{(l-1)}\right]\approx(1-\Gamma^{\alpha})c_{1}(\Delta
t_{l})^{\alpha}\approx(1-\Gamma^{\alpha})\mathrm{E}\left[\hat{P}_{l}^{(l)}-P\right],$

can be used to obtain the stopping criteria

(11)
$\left|\hat{Y}_{l}\right|=\left|\hat{P}_{l}^{(l)}-\hat{P}_{l}^{(l-1)}\right|\leq\frac{\Gamma^{\alpha}-1}{\sqrt{2}}\delta,$

which leads to
$\left|\mathrm{E}\left[\hat{P}_{l}^{(l)}-P\right]\right|^{2}\leq\delta^{2}/2$.

We summarize the adaptive MLMC algorithm [12] as follows.

Algorithm 1 Multilevel Monte Carlo method

Given $N_{I}\in\mathbb{N}$.

1: for $l=0,1,\dots$ do

2: Draw $N_{l}=N_{I}$ samples to obtain $\hat{V}_{l}$, $\hat{Y}_{l}$;

3: Define optimal $N_{k}^{\prime}$, $k=0,\dots,l$ using (8);

4: If $N_{k}^{\prime}>N_{k}$, generate extra samples and update $\hat{V}_{k}$,
$\hat{Y}_{k}$, $N_{k}=N_{k}^{\prime}$ for $k=0,\dots,l$;

5: if $l\geq 1$ & (11) is true then

5: Stop with $L=l$;

6: end if

7: $l=l+1$;

8: end for

For MC, we only consider the solution $\hat{P}^{(L)}=\hat{P}_{L}^{(L)}$ on the
finest level. Denoting $\mathsf{V}_{l}=\mathrm{Var}\left[P_{l}^{(l)}\right]$,
we require

(12)
$\mathrm{Var}\left[\hat{P}_{L}^{(L)}\right]=N_{L}^{-1}\mathsf{V}_{L}\leq\delta^{2}/2,$

in MC. For comparison, we require
$\mathrm{Var}\left[\hat{P}^{(L)}\right]=\sum_{l=0}^{L}N_{l}^{-1}V_{l}\leq\delta^{2}/2$
in MLMC.

By (12), the optimal number of samples $N_{l}^{\prime}$ can be taken as

(13) $N_{l}^{\prime}=\left\lceil
2\delta^{-2}\hat{\mathsf{V}}_{l}\right\rceil,$

where $\hat{\mathsf{V}}_{l}$ is the Monte Carlo estimate of $\mathsf{V}_{l}$
using initial $N_{I}$ samples. We use a similar a posterior estimate as the
stopping criteria

(14)
$\left|\hat{P}_{l}^{(l)}-\hat{P}_{l-1}^{(l-1)}\right|<\frac{\Gamma^{\alpha}-1}{\sqrt{2}}\delta.$

The adaptive MC algorithm can be given in the spirit of adaptive MLMC
Algorithm 1.

Algorithm 2 Monte Carlo method

Given $N_{I}\in\mathbb{N}$.

1: for $l=0,1,\dots$ do

2: Draw $N_{l}=N_{I}$ samples to obtain $\hat{\mathsf{V}}_{l}$,
$\hat{P}_{l}^{(l)}$;

3: Define $N_{l}^{\prime}$ using (13);

4: If $N_{l}^{\prime}>N_{l}$, generate extra samples at level $l$ and update
$\hat{\mathsf{V}}_{l}$, $\hat{P}_{l}^{(l)}$, $N_{l}=N_{l}^{\prime}$;

5: if $l\geq 1$ & (14) is true then

5: Stop with $L=l$;

6: end if

7: $l=l+1$;

8: end for

###  2.2 Cost analysis

In this section, we present cost theorems for MLMC and MC methods, as
prescribed in Algorithm 1 and Algorithm 2, which are crucial for the
hyperbolic equation applications. The proofs are given in Appendix B.1 and
B.2, respectively.

######  Theorem 2.2 (MLMC Computational Cost).

Assume that the MLMC estimators $\hat{Y}_{l}$ ($l\in\mathbb{N}$) defined in
(7) has the following properties:

  1. A1.

$\mathrm{E}\left[P^{(l)}-P\right]\leq c_{1}(\Delta t_{l})^{\alpha}$,

  2. A2.

$\mathrm{Var}\left[\hat{Y}_{l}\right]=N_{l}^{-1}V_{l}\leq
c_{2}N_{l}^{-1}(\Delta t_{l})^{\beta}$,

  3. A3.

$\mathfrak{C}\left[\hat{Y}_{l}\right]=N_{l}C_{l}\leq c_{3}N_{l}(\Delta
t_{l})^{-\gamma}$,

where $2\alpha\geq\gamma$, $\beta$, $\gamma$, $c_{1}$, $c_{2}$, $c_{3}$ are
positive constants. There exist $L\in\mathbb{N}^{+}$,
$N_{l}\in\mathbb{N}^{+}$, $l=0,...,L$, and a constant $c_{4}>0$ such that for
any $\delta<e^{-1}$, the multilevel estimator

(15) $\hat{P}^{(L)}=\sum_{l=0}^{L}\hat{Y}_{l},$

has the following $\mathrm{MSE}$ bound

(16)
$\mathrm{MSE}\equiv\mathrm{E}\left[\left(\hat{P}^{(L)}-\mathrm{E}\left[P\right]\right)^{2}\right]<\delta^{2},$

and the computational cost $\mathfrak{C}$ is bounded by

(17) $\displaystyle\mathfrak{C}$ $\displaystyle\leq\left\\{\begin{aligned}
&c_{4}\delta^{-2},&\beta>\gamma,\\\
&c_{4}\delta^{-2}(\log\delta^{-1})^{2},&\beta=\gamma,\\\
&c_{4}\delta^{-2-(\gamma-\beta)/\alpha},&0<\beta<\gamma.\end{aligned}\right.$

######  Theorem 2.3 (MC Computational Cost).

Assume that the MC estimators $\hat{P}_{l}^{(l)}$ ($l\in\mathbb{N}$) has the
following properties:

  1. A4.

$\mathrm{E}\left[P^{(l)}-P\right]\leq c_{1}(\Delta t_{l})^{\alpha}$,

  2. A5.

$\mathrm{Var}\left[\hat{P}_{l}^{(l)}\right]=N_{l}^{-1}\mathsf{V}_{l}\leq
c_{2}N_{l}^{-1}(\Delta t_{l})^{\beta_{0}}$,

  3. A6.

$\mathfrak{C}\left[\hat{P}_{l}^{(l)}\right]=N_{l}\mathfrak{C}\left[P_{l}^{(l)}\right]\leq
c_{3}N_{l}(\Delta t_{l})^{-\gamma}$,

where $2\alpha\geq\gamma$, $\beta_{0}$, $\gamma$, $c_{1}$, $c_{2}$, $c_{3}$
are positive constants. There exist $L\in\mathbb{N}^{+}$,
$N_{l}\in\mathbb{N}^{+}$, $l=0,...,L$, and a constant $c_{4}>0$ such that for
any $\delta<e^{-1}$, the MC estimator $\hat{P}_{L}^{(L)}$ has the following
$\mathrm{MSE}$ bound

(18)
$\mathrm{MSE}\equiv\mathrm{E}\left[\left(\hat{P}_{L}^{(L)}-\mathrm{E}\left[P\right]\right)^{2}\right]<\delta^{2},$

and the computational cost $\mathfrak{C}^{*}$ is bounded by

(19) $\displaystyle\mathfrak{C}^{*}$ $\displaystyle\leq\left\\{\begin{aligned}
&c_{4}\delta^{-2},&\beta_{0}>\gamma,\\\
&c_{4}\delta^{-2}\log\delta^{-1},&\beta_{0}=\gamma,\\\
&c_{4}\delta^{-2-(\gamma-\beta_{0})/\alpha},&0<\beta_{0}<\gamma.\end{aligned}\right.$

Comparing the results of Eq. 17 and Eq. 19, we can see that when
$0<\beta_{0}<\beta<\gamma$, the MLMC method outperforms the MC method.

##  3  Comparing MLMC and MC for the scalar advection equation with random
velocity

In sections § 3, § 4 and § 5, we use four examples to compare MLMC and MC
methods for both linear and nonlinear problems. In this section, we consider
the first example, the scalar advection equation with a time dependent random
velocity. For this example, we analyze the variances of the MC solutions and
MLMC corrections, and provide the complexity estimates of MLMC and MC methods
in Theorem 3.3 for two cases: 1) the number of random parameters is finite, 2)
the random parameter is modeled as white noise. In the first case, MLMC does
outperform MC, while in the second case, MLMC has the same order of
computational cost $\delta^{-3}$ as MC. We then use numerical experiments to
justify those theoretical results.

We first consider the scalar advection equation

(20) $u_{t}+a(t,\omega)u_{x}=0,$

with $x\in[-1,1]$, $t\in[0,T]$, initial condition $u(0,x)=u_{0}(x)$, and
periodic boundary condition $u(-1,t)=u(1,t)$.

We assume that the random variable is of the form
$\omega=(\omega_{1},\dots,\omega_{K})$, such that $K$ is the random dimension.
The random time dependent velocity $a=a(t,\omega)$ can be modeled as piecewise
constants over $K$ equal length time intervals

(21) $a(t,\omega)=\bar{a}+\omega_{k}\quad\text{ for }\quad\frac{k-1}{K}T\leq
t<\frac{k}{K}T,\quad k=1,\dots,K,$

where $\bar{a}=1$, $\omega_{k}\sim\mathcal{U}(-1,1)$. The corresponding
solution is denoted as $u(x,t,\omega)$.

Let $a(t,\omega(s))$ be a sample of $a(t,\omega)$. We denote $v^{n}_{i}(s)$ as
the numerical solution at $(x_{i},t_{n})$, with the random variable
$\omega(s)$, $x_{i}=-1+i\Delta x_{l}$, $t_{n}=n\Delta t_{l}$, where $\Delta
t_{l}$ is the time step and $\Delta x_{l}$ is the mesh size, with
$\frac{\Delta t_{l}}{\Delta x_{l}}=\kappa\leq 1$ by the CFL condition. In the
following, when no confusion arises, we also use notations $v^{n}_{i}$ or
$v^{n}$ for simplicity. We denote $X_{l}:=\frac{2}{\Delta x_{l}}$ and
$M_{l}:=\frac{T}{\Delta t_{l}}$ as the spatial and temporal degrees of
freedom.

Equation (20) can be discretized by the first order upwind scheme

(22) $\frac{v_{i}^{n+1}(s)-v_{i}^{n}(s)}{\Delta
t}+a(t_{n},\omega{(s)})\frac{v_{i}^{n}(s)-v_{i-1}^{n}(s)}{\Delta x}=0,\quad
a>0.$

###  3.1 Analysis

We now carry out an analysis for this simple model. Denoting
$u(\cdot,t,\omega(s))$ as the exact solution of (20) using the random variable
$\omega(s)$, which can be obtained by method of characteristics. At the final
time $T$, we have the following expression

(23)
$u\left(x,T,\omega{(s)}\right)=u\left(x-\left(\bar{a}+\omega_{K}{(s)}\right)\frac{T}{K},\frac{K-1}{K}T\right)=\dots=u_{0}\left(x-\bar{a}T-\frac{\sum_{k=1}^{K}\omega_{k}{(s)}}{K}T\right).$

######  Proposition 3.1.

Denote $P=u\left(\cdot,T,\omega(s)\right)$ as the exact solution and $P^{(l)}$
as the numerical solution using time step $\Delta t_{l}$ and mesh size $\Delta
x_{l}$ at final time $T$. The variances of $P$, $P^{(l)}$ satisfy the
following properties

  1. P1.

$\mathrm{Var}\left[P\right]:=\sum_{i}\mathrm{Var}\left[P(x_{i})\right]\Delta
x_{l}=O(K^{-1})$,

  2. P2.

$\mathrm{Var}\left[P^{(l)}-P\right]:=\sum_{i}\mathrm{Var}\left[(P^{(l)}-P)(x_{i})\right]\Delta
x_{l}=O((\Delta t_{l})^{2})$,

  3. P3.

$\mathrm{Var}\left[P^{(L)}\right]:=\sum_{i}\mathrm{Var}\left[P^{(l)}(x_{i})\right]\Delta
x_{l}=O(K^{-1})$.

######  Proof 3.2.

  1. (a)

Proof of P1. Defining $\bar{u}=u_{0}(x-\bar{a}T)$, we calculate $P-\bar{u}$ as
follows

(24)
$\left|(P-\bar{u})(x_{i})\right|^{2}=\left|u(x_{i},T,\omega(s))-\bar{u}(x_{i})\right|^{2}\leq\left\|u_{0}^{\prime}(\cdot)\right\|_{\infty}^{2}T^{2}\left|\frac{\sum_{k=1}^{K}\omega_{k}^{(s)}}{K}\right|^{2},$

and

(25)
$\mathrm{E}\left[\left|\frac{\sum_{k=1}^{K}\omega_{k}^{(s)}}{K}\right|^{2}\right]=\mathrm{Var}\left[\frac{\sum_{k=1}^{K}\omega_{k}^{(s)}}{K}\right]=\frac{\mathrm{Var}\left[\omega_{1}\right]}{K}.$

Property P1 then follows by

(26)
$\mathrm{Var}\left[P\right]=\sum_{i}\mathrm{Var}\left[\left(P-\bar{u}\right)(x_{i})\right]\Delta
x_{l}\leq\sum_{i}\mathrm{E}\left[\left(P-\bar{u}\right)^{2}(x_{i})\right]\Delta
x_{l}=O(K^{-1}).$

  2. (b)

Proof of P2. We define the following intermediate variable with respect to the
upwind scheme (22),

$\tilde{u}(x_{i},t_{n+1},\omega(s)):=u(x_{i},t_{n},\omega{(s)})-a(t_{n},\omega(s))\frac{\Delta
t_{l}}{\Delta
x_{l}}(u(x_{i},t_{n},\omega{(s)})-u(x_{i-1},t_{n},\omega{(s)})).$

The local truncation error of the upwind scheme is

(27) $\displaystyle\tau_{i}^{n}$
$\displaystyle:=u(x_{i},t_{n+1},\omega{(s)})-\tilde{u}(x_{i},t_{n+1},\omega{(s)})$
$\displaystyle=\frac{a(t_{n},\omega(s))(\Delta
t_{l})^{2}}{2}\left(a(t_{n},\omega(s))u_{xx}(x_{i},t_{\eta})-\frac{\Delta
x_{l}}{\Delta t_{l}}u_{xx}(x_{\xi},t_{n})\right)$
$\displaystyle\leq\frac{\mathcal{C}_{1}a(t_{n},\omega(s))}{2}\left(a(t_{n},\omega(s))+\frac{\Delta
x_{l}}{\Delta t_{l}}\right)(\Delta t_{l})^{2}$
$\displaystyle\leq\left(2+\frac{1}{\kappa}\right)\mathcal{C}_{1}(\Delta
t_{l})^{2},$

where $\mathcal{C}_{1}=\sup|u_{xx}|$. Then

$\left\|\tau^{n}\right\|_{2}\leq\left(2+\frac{1}{\kappa}\right)\mathcal{C}_{1}\sqrt{X_{l}}(\Delta
t_{l})^{2}.$

The stability of the upwind scheme is obvious since

(28)
$\left\|v^{n+1}\right\|_{2}^{2}=\sum_{i}\left|v_{i}^{n+1}\right|^{2}\leq\sum_{i}(1-\kappa
a(t_{n},\omega(s)))\left|v_{i}^{n}\right|^{2}+\kappa
a(t_{n},\omega(s))\left|v_{i-1}^{n}\right|^{2}=\left\|v^{n}\right\|_{2}^{2},$

where the CFL condition $\kappa\leq 1$ is used. The error $P^{(l)}-P$ can be
bounded by

(29) $\displaystyle\left\|P^{(l)}-P\right\|_{2}$
$\displaystyle=\left\|v^{M_{l}}-u\left(\cdot,t_{M_{l}},\omega(s)\right)\right\|_{2}$
$\displaystyle\leq\left\|v^{M_{l}}-\tilde{u}\left(\cdot,t_{M_{l}},\omega(s)\right)\right\|_{2}+\left\|u\left(\cdot,t_{M_{l}},\omega(s)\right)-\tilde{u}\left(\cdot,t_{M_{l}},\omega{(s)}\right)\right\|_{2}$
$\displaystyle(\text{stability})$
$\displaystyle\leq\left\|v^{M_{l}-1}-u\left(\cdot,t_{M_{l}-1},\omega(s)\right)\right\|_{2}+\left\|\tau^{M_{l}-1}\right\|_{2}$
$\displaystyle\leq\cdots$ $\displaystyle\leq
0+\sum_{j=0}^{M_{l}-1}\left\|\tau^{j}\right\|_{2}$
$\displaystyle\leq\left(2+\frac{1}{\kappa}\right)\mathcal{C}_{1}T\sqrt{X_{l}}\Delta
t_{l}.$

Property P2 follows by

(30)
$\mathrm{Var}\left[P^{(l)}-P\right]\leq\sum_{i}\mathrm{E}\left[\left(P^{(l)}-P\right)^{2}(x_{i})\right]\Delta
x_{l}=\mathrm{E}\left[\left\|P^{(l)}-P\right\|_{2}^{2}\right]\Delta
x_{l}=O((\Delta t_{l})^{2}).$

  3. (c)

Proof of P3. Property P3 follows by the assumption that $\Delta
t_{l}=O(K^{-1})$, and the inequality

$\left(\sqrt{\mathrm{Var}\left[P^{(l)}-P\right]}-\sqrt{\mathrm{Var}\left[P\right]}\right)^{2}\leq\mathrm{Var}\left[P^{(l)}\right]\leq\left(\sqrt{\mathrm{Var}\left[P^{(l)}-P\right]}+\sqrt{\mathrm{Var}\left[P\right]}\right)^{2}.$

We compute MLMC corrections, which are the difference of solutions
$P^{(l)}_{l}$ and $P^{(l-1)}_{l}$, the numerical solutions at final time $T$
of (20) on the $l$th level, with time step $\Delta t_{l}$, $\Delta
t_{l-1}=\Gamma\Delta t_{l}$ respectively. We have the following theorem for
the computational costs of MLMC and MC for the scalar random advection
equation.

######  Theorem 3.3.

We consider two cases,

  * •

Case I: When the number of random variables $K$ is fixed, and $M_{0}\geq K$,
we have $\beta_{0}=0$, $\beta=2$, i.e.

(31)
$\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1),\quad\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=O((\Delta
t_{l})^{2}).$

Given the accuracy threshold $\delta>0$, the cost of MLMC is
$O\left(\delta^{-2}\left(\log\delta^{-1}\right)^{2}\right)$ and the cost of MC
is $O\left(\delta^{-4}\right)$.

  * •

Case II: The random velocity is modeled as white noise, namely the number of
random variables depends on $l$, $K_{l}=\frac{T}{\Delta t_{l}}=M_{l}$ for each
level $l$. We have $\beta_{0}=\beta=1$, i.e.

(32) $\mathrm{Var}\left[P_{l}^{(l)}\right]=O(\Delta
t_{l}),\quad\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=O(\Delta
t_{l}).$

Given the accuracy threshold $\delta>0$, the computational costs of MLMC and
MC are both $O\left(\delta^{-3}\right)$.

######  Proof 3.4.

The values of $\beta_{0}$ can be obtained directly from Property P3. In the
following, we will explain how random variables are used in both cases and
evaluate $\beta$.

  * •

Proof of Case I. We first generate $K$ random variables for one sample on
level $l$. Since $K$ is fixed, we use the same random variables
$\left\\{\omega_{k}(s)\right\\}_{k=1}^{K}$ to calculate $P_{l}^{(l)}(s)$ and
$P_{l}^{(l-1)}(s)$. Therefore similar error analysis as in (27), (29), (30)
leads to the following estimate

(33)
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]\leq\mathrm{E}\left[\left\|P_{l}^{(l)}-P_{l}^{(l-1)}\right\|_{2}^{2}\right]=O((\Delta
t_{l})^{2}).$

Thus Property A2 holds with $\beta=2$. The computational costs of MLMC and MC
follow by the application of Eq. 17, Eq. 19, and the fact that $\alpha=1$,
$\gamma=2$, $\beta_{0}=0$, $\beta=2$.

In this case, $\left\\{\omega_{k}(s)\right\\}_{k=1}^{K}$ can be directly used
to calculate $P_{l}^{(l-1)}$, and the numerical implementation of Algorithm 1
is straightforward.

  * •

Proof of Case II. We need to generate $K_{l}=T/\Delta t_{l}$ random variables
for one sample on level $l$, denoted by
$\left\\{\omega_{k}\right\\}_{i=1}^{K_{l}}$. For example, let $a$ be a uniform
white noise

(34) $a(t)=\bar{a}+\omega(t),\quad\omega(t)\sim\mathcal{U}(-1,1),$

with correlation
$B(t)=\mathrm{Var}\left[\omega\right]\delta(t)=\frac{1}{3}\delta(t)$. At each
time $t_{n-1}$, $a(t_{n-1})=\bar{a}+\omega_{n}$, $n=1,\dots,K_{l}$. The
computation of $P^{(l)}_{l}$ is also straightforward, while the computation of
$P^{(l-1)}_{l}$ is a little tricky. To compute $P^{(l-1)}_{l}$, we need to
construct $K_{l-1}$ random variables from the $K_{l}$ random variables, such
that

    1. 1.

They are strongly correlated with $\left\\{\omega_{k}\right\\}_{k=1}^{K_{l}}$,
so that $\mathrm{Var}\left[P^{(l)}_{l}-P^{(l-1)}_{l}\right]$ is kept small;

    2. 2.

The distribution of the coarse random variables is the same as its finer
counterpart, so that
$\mathrm{E}\left[P^{(l-1)}_{l}\right]=\mathrm{E}\left[P^{(l-1)}_{l-1}\right]$
still holds.

To that end, we use the following trick [24] to construct those coarse random
variables

(35) $\frac{\tilde{a}(t_{(n-1)\Gamma})}{2}=\left(\max_{1\leq
j\leq\Gamma}\left\\{\frac{a(t_{n\Gamma-j})}{2}\right\\}\right)^{\Gamma}\sim\mathcal{U}(0,1),\quad
n=1,\dots,K_{l-1}(=K_{l}/\Gamma).$

The inequality (33) does not work since the random velocities
$\left\\{\tilde{a}(t_{n-1})\right\\}_{n=1}^{K_{l-1}}$,
$\left\\{a(t_{n-1})\right\\}_{n=1}^{K_{l}}$ are not the same now, although
strongly correlated. Nevertheless, the Cauchy-Schwarz inequality implies that

(36)
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]\leq\left(\sqrt{\mathrm{Var}\left[P_{l}^{(l)}\right]}+\sqrt{\mathrm{Var}\left[P_{l}^{(l-1)}\right]}\right)^{2}=O(\Delta
t_{l}).$

The Property A2 holds with $\beta=1$. The computational costs of MLMC and MC
follow by the application of Eq. 17, Eq. 19, and the fact that $\alpha=1$,
$\gamma=2$, $\beta_{0}=1$, $\beta=1$.

###  3.2 Numerical results: Case I

In the numerical experiments for Case I, we choose two constants $K=1$ and
$K=32$. We take

(37) $\kappa=\frac{\Delta t_{l}}{\Delta x_{l}}=0.5,\quad\Gamma=2,\quad\Delta
x_{0}=\frac{1}{32},\quad N_{I}=500,$

where $\Gamma$ is the refinement ratio, $N_{I}$ is the initial number of
sampling, the CFL condition is satisfied since $|a|\leq 2$. We also choose the
initial condition

(38) $u_{0}(x)=\frac{\sin(x)+1}{2},$

and the periodic boundary condition.

In view of our theoretical results, we will show the computational costs of MC
and MLMC vs. the accuracy level $\delta$, also their variances vs. level. Let
$N_{l}$, $N_{l}^{\ast}$ represent the number of samples of MLMC and MC
respectively. The computational cost of MLMC can be estimated as the total
number of time steps and gird points on all levels normalized by that on the
level $0$,

(39)
$\mathfrak{C}=N_{0}+\sum_{l=1}^{L}N_{l}\left(\Gamma^{2l}+\Gamma^{2(l-1)}\right).$

Namely, the computational cost of one sample on level $0$ ($C_{0}$) is
regarded as $1$ unit, and $C_{l+1}/C_{l}=\Gamma^{2},l\geq 1$ since we refine
the grid in both spatial and temporal dimensions.

The computational cost of MC can be estimated similiarly as

(40) $\mathfrak{C}^{\ast}=N_{0}^{\ast}+\sum_{l=1}^{L}N_{l}^{\ast}\Gamma^{2l}.$

In Fig. 1 and Fig. 2, the numerical results are shown in solid lines, while
the theoretical order is plotted in dashed lines. We observe that as long as
$K$ is a constant, we have $\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1)$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=O((\Delta t_{l})^{2})$,
i.e. $\beta_{0}=0$, $\beta=2$. The costs of MC and MLMC are of order
$O(\delta^{-4})$ and $O(\delta^{-2})$, respectively. In this case, to achieve
the same accuracy, MLMC outperforms MC by two orders of magnitude in terms of
the computational cost.

  1. 1.

$K=1$

Figure 1: Example 1: Case I, $K=1$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variance on different levels (right).

  2. 2.

$K=32$

Figure 2: Example 1: Case I, $K=32$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variance on different levels (right).

###  3.3 Numerical results: Case II

In this case, we need to generate $K_{l}=T/\Delta t_{l}$ random variables
$\omega_{k}\sim\mathcal{U}(-1,1)$ on the $l$ th level. Using (35), we can
construct $K_{l-1}$ random variables and use them to calculate
$P_{l}^{(l-1)}$. Fig. 3 shows that both $\mathrm{Var}\left[P_{l}^{(l)}\right]$
and $\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ decays with order
$O(\Delta t_{l})$, i.e. $\beta_{0}=\beta=1$, which is consistent with the
derivation in (36). The costs of MC and MLMC are both of the order
$O(\delta^{-3})$, while MLMC appears to have a larger constant than MC.

Figure 3: Example 1: Case II, the computational cost w.r.t accuracy threshold
$\delta$ (left) and variance on different levels (right).

##  4 Comparing MLMC and MC for the Euler equation and the shallow water
equation

In this section, we extend our study to more practical examples such as the
Euler equation with a time dependent random adiabatic constant, and the
shallow water equations with space dependent random tomography. The Euler
equation example is qualitatively similar to the scalar advection example,
albeit more complex. However, the shallow water equation example is
qualitatively distinct, and has computational cost $\delta^{-4}$ if the random
tomography is modeled as white noise. In [25], Mishra et. al. applied MLMC to
Euler equations with random initial data, while we consider Euler equations
with time dependent random parameters in § 4.1. Mishra et. al. also applied
MLMC to shallow water equations with uncertain topography in [26], where the
bottom topography is treated as a fixed number of random variables with
uniform distribution. That is similar to our Case I, but we also consider the
white noise Case II in § 4.2.

###  4.1 Euler equation with random adiabatic constant $\lambda(t,\omega)$

In this subsection, we consider the Euler equations as a prototypical
nonlinear partial differential equations with time dependent random
parameters. When the number of random variables $K$ is fixed, the analysis of
Property A2 is similar as (33); while the white noise case is now hard to
analyze due to nonlinearity, we can numerically estimate
$\mathrm{Var}\left[P_{l}^{(l)}\right]$,
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ and their decay orders
$\beta_{0}$, $\beta$.

Given the Euler equations with gravity

(41) $\left\\{\begin{aligned} \rho_{t}+(\rho v)_{x}&=0,\\\ (\rho v)_{t}+(v\rho
v)_{x}+p_{x}&=-\rho\phi_{x},\\\ E_{t}+[(E+p)v]_{x}&=-\rho
v\phi_{x},\end{aligned}\right.$

where

$E=\rho e+\frac{\rho}{2}v^{2},\quad\phi(x)=gx.$

Here $\rho(t,x)$ is the density, $v(t,x)$ is the velocity, $E(t,x)$ is the
total energy, $p(t,x)$ is the pressure, $e(t,x)$ is the specific internal
energy, and $g$ is the gravitational constant. The adiabatic constant
$\lambda$ can be modeled as a time dependent random parameter
$\lambda(t,\omega)$ in the equation of state

$p=\rho e(\lambda-1)=A\rho^{\lambda}.$

Problems with uncertain $\lambda$ were considered in, for example, [29].

We impose the initial condition

(42) $\rho_{0}(x)=2-\phi(x),\quad
p_{0}(x)=A_{0}\rho(x)^{\lambda_{0}},\quad\lambda_{0}=\frac{4}{3},\quad
A_{0}=1,$

and the Neumann boundary condition on $x\in[0,2]$.

To be more precise, we let

(43) $\lambda=\bar{\lambda}+\omega_{k}\quad\text{ for }\quad\frac{k-1}{K}T\leq
t<\frac{k}{K}T,\quad k=1,\dots,K,$

with $\bar{\lambda}=\frac{4}{3}$,
$\omega_{k}\sim\mathcal{U}(-\frac{1}{3},\frac{1}{3})$.

We use a well-balanced scheme introduced by Käppeli and Mishra [19], where the
hydrostatic steady state of (41) is

(44) $v\equiv 0,\quad p_{x}=-\rho\phi_{x},\quad h+\phi=\text{const},\quad
h=e+\frac{p}{\rho}.$

The numerical parameters are chosen as follows

(45) $\kappa=\frac{\Delta t_{l}}{\Delta x_{l}}=0.1,\quad\Gamma=2,\quad\Delta
x_{0}=0.25,\quad T=2,\quad N_{I}=500.$

Here we show the numerical results for Case I: $\lambda$ with fixed $K$, and
Case II: white noise $\lambda$.

####  4.1.1 Case I

Fig. 4 ($K=1$) and Fig. 5 ($K=8$) show that
$\mathrm{Var}\left[P_{l}^{(l)}\right]\approx O(1)$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]\approx O(\left(\Delta
t\right)^{1.5})$, i.e. $\beta_{0}=0$, $\beta=1.5$. The costs of MC and MLMC
are $O(\delta^{-4})$ and $O(\delta^{-2.5})$ respectively. Both results are
consistent with the analysis.

  1. 1.

$K=1$

Figure 4: Example 2: Case I, $K=1$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variances on different levels (right).

  2. 2.

$K=8$

Figure 5: Example 2: Case I, $K=8$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variance on different levels (right).

####  4.1.2 Case II

In this case, we generate white noise $\lambda$ with $K_{l}=T/\Delta t_{l}$ on
the $l$th level, and use similar trick as (35) to construct coarser random
variables.

From Fig. 6, we observe that $\mathrm{Var}\left[P_{l}^{(l)}\right]$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ are both $O(\Delta
t_{l})$, i.e. $\beta_{0}=\beta=1$. Comparing with the variances in Fig. 4,
Fig. 5, it is evident that the temporal white noise leads to the decay of
$\mathrm{Var}\left[P_{l}^{(l)}\right]$. The costs of MC and MLMC are both
$O(\delta^{-3})$.

Figure 6: Example 2: Case II, the computational cost w.r.t accuracy threshold
$\delta$ (left) and variance on different levels (right).

###  4.2 Shallow water equations with random tomography $B(x,\omega)$

In this example, we consider space dependent random parameters instead of time
dependent random parameters, using the shallow water equations

(46) $\left\\{\begin{aligned} &h_{t}+(hu)_{x}=0,\\\
&(hu)_{t}+\left(hu^{2}+\frac{gh^{2}}{2}\right)_{x}=-ghB_{x}.\end{aligned}\right.$

Here $h(t,x)$ is the height, $u(t,x)$ is the mean velocity, with the
gravitational acceleration constant $g$, and the bottom topography $B(x)$. In
view of the measurement error [22], we can introduce the following random
tomography $B(x,\omega)$ defined via $K$ random variables, such that

(47) $B(x,\omega)=1+\omega_{k}\sin(\pi x)\quad\text{ for
}\quad\frac{k-1}{K}\leq x<\frac{k}{K},\quad k=1,\dots,K,$

with $\omega_{k}\sim\mathcal{U}(0,1)$ and the steady state

(48) $q:=hu\equiv\widehat{q}=\text{const},\quad
E:=\frac{u^{2}}{2}+g(h+B)\equiv\widehat{E}=\text{const}.$

We treat the still water case with $q\equiv 0$ ($v\equiv 0$) using the well-
balanced central-upwind scheme [21]. The initial conditions

(49) $h(x,0)=5+e^{\cos(2\pi x)},\quad u(x,0)=\sin(\cos(2\pi x)),$

and the periodic boundary condition on $x\in[0,1]$ are imposed.

We use the following numerical parameters

(50) $\kappa=\frac{\Delta t}{\Delta x}=0.05,\quad\Gamma=2,\quad\Delta
x_{0}=\frac{1}{64},\quad T=0.1,\quad N_{I}=500.$

Similarly as before, we have two qualitatively different cases: Case I:
$B(x,\omega)$ with fixed $K$, and Case II: spatial white noise $B(x,\omega)$.

####  4.2.1 Case I

We use $K=1$ in Fig. 7 and $K=8$ in Fig. 8. We observe that
$\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1)$ ($\beta_{0}=0$), the numerical
values for the decay orders of
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ in these two experiments
are $\beta=2.5\,(K=1)$ and $\beta=1.0\,(K=8)$ respectively.

The cost of MC is approximately $O(\delta^{-4})$ since $\beta_{0}=0$. In the
case $K=1$, since $\beta=2.5>\gamma$, the cost of MLMC is approximately
$O(\delta^{-2})$. While in the case $K=8$, the cost of MLMC is approximately
$O(\delta^{-3})$. As we expected, MLMC is better than MC method by one or two
orders of magnitude.

  1. 1.

$K=1$

Figure 7: Example 3: Case I, $K=1$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variance on different levels (right).

  2. 2.

$K=8$

Figure 8: Example 3: Case I, $K=8$, the computational cost w.r.t accuracy
threshold $\delta$ (left) and variance on different levels (right).

####  4.2.2 Case II

When $B(x,\omega)$ is modeled by spatially dependent white noise, we have
$K_{l}=X_{l}$ for the $l$th level. When a sample $P_{l}^{(l)}$ is computed on
level $l$, the value of $B(x_{i},\omega)$ is generated by

(51) $B(x_{i-1},\omega)=1+\omega_{i}\sin(\pi
x_{i-1}),\quad\omega_{i}\sim\mathcal{U}(0,1),\quad i=1,\dots,X_{l}.$

To calculate $P_{l}^{(l-1)}$, we first construct
$\left\\{\tilde{\omega}_{i}\right\\}_{i=1}^{X_{l-1}}$ using (35), then let

(52) $\tilde{B}(x_{(i-1)\Gamma},\omega)=1+\tilde{\omega}_{i}\sin(\pi
x_{(i-1)\Gamma}),\quad i=1,\dots,X_{l-1}(=X_{l}/\Gamma).$

Fig. 9 shows that the decays of $\mathrm{Var}\left[P_{l}^{(l)}\right]$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ are very slow, namely, we
have $\beta\simeq 0$ and $\beta_{0}\simeq 0$. In this case, both MC and MLMC
methods are expensive. In fact, MLMC has a higher cost compared to MC.

Figure 9: Example 3: Case II, the computational cost w.r.t accuracy threshold
$\delta$ (left) and variance on different levels (right).

##  5 Comparing MLMC and MC for the linear Jin-Xin relaxation model with
random choice method

While the first three examples in § 3 and § 4 are hyperbolic equations with
random parameters, the last example is of different nature. In this example,
we apply a randomized algorithm to a multiscale deterministic equation–the
Jin-Xin hyperbolic relaxation model [18]. This algorithm is in the spirit of
the Glimm’s random choice method [13], and the random variables are in both
spatial and temporal dimensions. We identify similar behavior of MLMC in this
example as the third example. This simple model can shed some light on the
application of randomized algorithms to more complicated problems.

In particular, we deal with the following one-dimensional linear Jin-Xin
relaxation model [18]

(53) $\left\\{\begin{aligned} \frac{\partial u}{\partial t}+\frac{\partial
v}{\partial x}&=0,\\\ \frac{\partial v}{\partial t}+a\frac{\partial
u}{\partial x}&=-\frac{1}{\varepsilon}(v-bu),\end{aligned}\right.$

where $(x,t)\in[-1,1]\times[0,T]$ with the periodic boundary condition, and
the initial data $u(x,0)=u_{0}(x)$, $\quad v(x,0)=v_{0}(x)$. $a>0$,
$b\in\mathbb{R}$, and $\varepsilon>0$ is the relaxation parameter. Asymptotic
analysis [18, 28] shows that when $\varepsilon$ is small, the stability
condition $a\geq b^{2}$ ensures the well-posedness of model. In the limit
$\varepsilon\rightarrow 0$, the second equation of (53) relaxes to the local
equilibrium $v=bu$, and the first equation of (53) reduces to the scalar
advection equation $u_{t}+bu_{x}=0.$ The Jin-Xin model is a prototypical
multiscale hyperbolic equation, and can be treated by the so-called asymptotic
preserving (AP) method [17].

In the following, we will develop the asymptotic preserving Monte Carlo method
(APMC) for the equation (53), by combining the Glimm’s random choice method
[13] and the asymptotic preserving principle. A random algorithm is
interesting not only for theoretical interest [13] but also of practical
interest [4, 7, 10]. We will compare the APMC method with its multilevel
version (MLAPMC). We will demonstrate that, since the random variables in the
APMC method involves both spatial and temporal dimensions, the cost of MLAPMC
is of the same order as APMC.

###  5.1 Asymptotic preserving Monte Carlo method (APMC)

The equation (53) is stiff in the sense that the time step of any explicit
method has to satisfy $\Delta t=O(\varepsilon)$ by the stability requirement,
which is prohibitively expensive as $\varepsilon\rightarrow 0$. The asymptotic
preserving (AP) method [15] allows to choose $\Delta t$ independent of
$\varepsilon$, while it still attains the correct asymptotic limit. To be
precise, (53) can be split into two steps,

(54) $\text{(convection step)}\quad\left\\{\begin{aligned} \frac{\partial
u}{\partial t}+\frac{\partial v}{\partial x}&=0,\\\ \frac{\partial v}{\partial
t}+a\frac{\partial u}{\partial x}&=0,\end{aligned}\right.$

which is a linear convection equation, and

(55) $\text{(relaxation step)}\quad\left\\{\begin{aligned} \frac{\partial
u}{\partial t}&=0,\\\ \frac{\partial v}{\partial
t}&=-\frac{1}{\varepsilon}(v-bu),\end{aligned}\right.$

which is a relaxation equation with a stiff source only in the equation for
$v$.

We can construct a random choice method for the split system,

  * •

convection step: (54) can be transformed to a diagonal form as

(56) $\left\\{\begin{aligned} \frac{\partial\rho}{\partial
t}+\sqrt{a}\frac{\partial\rho}{\partial x}=0,\\\ \frac{\partial l}{\partial
t}-\sqrt{a}\frac{\partial l}{\partial x}=0,\end{aligned}\right.$

by introducing the characteristic variables

(57) $\rho=\sqrt{a}u+v,\quad l=\sqrt{a}u-v.$

The explicit upwind scheme to (56) writes

(58) $\displaystyle\frac{\rho_{i}^{n+\frac{1}{2}}-\rho_{i}^{n}}{\Delta t}$
$\displaystyle=-\sqrt{a}\frac{\rho_{i}^{n}-\rho_{i-1}^{n}}{\Delta x},$
$\displaystyle\frac{l_{i}^{n+\frac{1}{2}}-l_{i}^{n}}{\Delta t}$
$\displaystyle=\sqrt{a}\frac{l_{i+1}^{n}-l_{i}^{n}}{\Delta x},$

namely,

(59) $\displaystyle\rho_{i}^{n+\frac{1}{2}}$
$\displaystyle=(1-\nu)\rho_{i}^{n}+\nu\rho_{i-1}^{n},$ $\displaystyle
l_{i}^{n+\frac{1}{2}}$ $\displaystyle=(1-\nu)l_{i}^{n}+\nu l_{i+1}^{n},$

where $\displaystyle\nu=\sqrt{a}\frac{\Delta t}{\Delta x}$ is the CFL number.
$u_{i}^{n+\frac{1}{2}}$ and $v_{i}^{n+\frac{1}{2}}$ can be obtained by

(60) $\displaystyle u_{i}^{n+\frac{1}{2}}$
$\displaystyle=\frac{\rho_{i}^{n+\frac{1}{2}}+l_{i}^{n+\frac{1}{2}}}{2\sqrt{a}},$
$\displaystyle v_{i}^{n+\frac{1}{2}}$
$\displaystyle=\frac{\rho_{i}^{n+\frac{1}{2}}-l_{i}^{n+\frac{1}{2}}}{2}.$

The Monte Carlo method can be applied to update $\rho$ and $l$, since (59) are
convex combinations. For example, we generate samples according to the
following distribution

(61) $\displaystyle\rho_{i}^{n+\frac{1}{2}}$
$\displaystyle=\left\\{\begin{array}[]{ll}{\rho_{i}^{n},}&{\text{with
probability }1-\nu},\\\ {\rho_{i-1}^{n},}&{\text{with probability
}\nu},\end{array}\right.$ $\displaystyle l_{i}^{n+\frac{1}{2}}$
$\displaystyle=\left\\{\begin{array}[]{ll}{l_{i}^{n},}&{\text{with probability
}1-\nu},\\\ {l_{i+1}^{n}.}&{\text{with probability }\nu}.\end{array}\right.$

In the implementation, we let $\rho_{i}^{n+\frac{1}{2}}=\rho_{i}^{n}$ and
sample a random variable $\xi_{i}^{n}\sim\mathcal{U}(0,1)$. If
$\xi_{i}^{n}<\nu$, we update $\rho_{i}^{n+\frac{1}{2}}=\rho_{i-1}^{n}$.

  * •

relaxation step: we use the fully implicit Euler’s method to treat the stiff
relaxation equation (55),

(62) $\displaystyle\frac{u_{i}^{n+1}-u_{i}^{n+\frac{1}{2}}}{\Delta t}$
$\displaystyle=0,$
$\displaystyle\frac{v_{i}^{n+1}-v_{i}^{n+\frac{1}{2}}}{\Delta t}$
$\displaystyle=-\frac{1}{\varepsilon}\left(v_{i}^{n+1}-bu_{i}^{n+1}\right).$

Thus

(63) $\displaystyle u_{i}^{n+1}$ $\displaystyle=u_{i}^{n+\frac{1}{2}},$
$\displaystyle v_{i}^{n+1}$
$\displaystyle=\frac{\varepsilon}{\varepsilon+\Delta
t}v_{i}^{n+\frac{1}{2}}+\frac{\Delta t}{\varepsilon+\Delta
t}bu_{i}^{n+\frac{1}{2}}.$

$v_{i}^{n+1}$ is a convex combination of $v_{i}^{n+\frac{1}{2}}$ and
$bu_{i}^{n+\frac{1}{2}}$, which has the following Monte Carlo approximation

(64)
$v_{i}^{n+1}=\left\\{\begin{array}[]{ll}{v_{i}^{n+\frac{1}{2}},}&{\text{with
probability }\displaystyle p,}\\\ {bu_{i}^{n+\frac{1}{2}}.}&{\text{with
probability }1-\displaystyle p,}\end{array}\right.$

with $\displaystyle p=\frac{\varepsilon}{\varepsilon+\Delta t}$.

The APMC algorithm can be summarized as follows

Algorithm 3 Asymptotic Preserving Monte Carlo (APMC) - fully-random version

Given the number of samples $N\in\mathbb{N}$, $T>0$, time step $\Delta t$,
mesh size $\Delta x$ (or equivalently, $\nu=\sqrt{a}\Delta t/\Delta x$).
Generating $N$ samples $\left\\{u(k),v(k)\right\\}_{k=1}^{N}$.

1: Initialize $u^{0}_{i}(k)=u_{0}(x_{i}),v^{0}_{i}(k)=v_{0}(x_{i})$,
$k=1,\dots,N$;

2: for $n=0,\dots,T/\Delta t-1$ do

3: for $k=1,\dots,N$ do

4: Compute $\rho^{n}_{i}(k),l^{n}_{i}(k)$ as in (57);

5: Sample $\rho^{n+\frac{1}{2}}_{i}(k)$ and $l^{n+\frac{1}{2}}_{i}(k)$ as in
(61);

6: Compute $u^{n+\frac{1}{2}}_{i}(k),v^{n+\frac{1}{2}}_{i}(k)$ as in (60);

7: $u^{n+1}_{i}(k)=u^{n+\frac{1}{2}}_{i}(k)$;

8: Sample $v^{n+1}_{i}(k)=v^{n+\frac{1}{2}}_{i}(k)$ as in (64);

9: end for

10: $\hat{u}_{i}^{n}=\frac{1}{N}\sum_{k=1}^{N}u_{i}^{n}(k)$,
$\hat{v}_{i}^{n}=\frac{1}{N}\sum_{k=1}^{N}v_{i}^{n}(k)$;

11: end for

We can replace the convection Step 5 in Algorithm 3 by the deterministic
scheme (59). The resulting algorithm is called the _semi-random_ version of
APMC. Also, we can obtain the deterministic method if we replace the Step 5 by
(59) and the Step 8 by (62), respectively.

###  5.2 Variance analysis

The deterministic method is asymptotic preserving since it has the correct
$\varepsilon\to 0$ limit, as proved in Appendix C. To apply the Monte Carlo
method as $\varepsilon\rightarrow 0$, the variance should be uniformly
bounded, i.e.

(65)
$\mathrm{Var}\left[u_{i}^{n}\right],\mathrm{Var}\left[v_{i}^{n}\right]\leq C,$

where $C$ is a constant independent of $\varepsilon$.

We present the following lemmas for the expectations and variances. The proofs
are given in Appendix D.

######  Lemma 5.1.

Under the condition $a\geq b^{2}$, the energy

(66)
$\displaystyle\mathcal{E}_{n}:=\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\right)^{2}\right]\Delta
x,$

is monotonically decreasing,

(67) $\mathcal{E}_{n}\leq\mathcal{E}_{n-1}\leq\cdots\leq\mathcal{E}_{0}.$

Moreover, the $l_{2}$ norms of
$\mathrm{E}\left[u^{n}\right],\mathrm{E}\left[v^{n}\right]$ remain bounded in
the sense that

(68)
$\left\|\mathrm{E}\left[u^{n}\right]\right\|_{2}=\sqrt{\sum_{i}\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}\Delta
x}\leq\sqrt{\frac{\mathcal{E}_{0}}{a-b^{2}}},$ (69)
$\left\|\mathrm{E}\left[v^{n}\right]\right\|_{2}=\sqrt{\sum_{i}\left(\mathrm{E}\left[v_{i}^{n}\right]\right)^{2}\Delta
x}\leq\sqrt{\frac{2a\mathcal{E}_{0}}{a-b^{2}}}.$

######  Remark 5.2.

We have the following estimate for
$\mathrm{E}\left[u_{i+1}^{n}-u_{i}^{n}\right]$ and
$\mathrm{E}\left[v_{i+1}^{n}-v_{i}^{n}\right]$,

(70)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i+1}^{n+1}-u_{i}^{n+1}\right]\right)^{2}+\left(\mathrm{E}\left[b(u_{i+1}^{n+1}-u_{i}^{n+1})-(v_{i+1}^{n+1}-v_{i}^{n+1})\right]\right)^{2}\right]\Delta
x$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\left(u_{i+1}^{0}-u_{i}^{0}\right)^{2}+\left(b\left(u_{i+1}^{0}-u_{i}^{0}\right)-\left(v_{i+1}^{0}-v_{i}^{0}\right)\right)^{2}\right]\Delta
x:=\mathcal{C}_{0}(\Delta x)^{2}.$

Similar as (68), (69), we have

(71)
$\sqrt{\sum_{i}\left(\mathrm{E}\left[\frac{u_{i+1}^{n+1}-u_{i}^{n+1}}{\Delta
x}\right]\right)^{2}\Delta
x}\leq\sqrt{\frac{\mathcal{C}_{0}}{a-b^{2}}},\quad\sqrt{\sum_{i}\left(\mathrm{E}\left[\frac{v_{i+1}^{n+1}-v_{i}^{n+1}}{\Delta
x}\right]\right)^{2}\Delta x}\leq\sqrt{\frac{2a\mathcal{C}_{0}}{a-b^{2}}}.$

######  Lemma 5.3.

As $\varepsilon\rightarrow 0^{+}$,
$\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\to 0$, and for any fixed
$\varepsilon>0$,
$\sum_{j}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j}-v_{i}^{j}\right]\right)^{2}\Delta
x$ is uniformly bounded in the sense that

(72)
$\sum_{j=1}^{\infty}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j}-v_{i}^{j}\right]\right)^{2}\Delta
x\leq\frac{\displaystyle p^{2}}{1-\displaystyle p^{2}}\mathcal{E}_{0}.$

In particular,
$\sum_{i}\left(\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\right)^{2}\Delta
x\to 0$ as $n\rightarrow\infty$.

To better understand the variances introduced in the convection step and in
the relaxation step, respectively, we first consider the semi-random APMC with
the random relaxation step and the deterministic convection step. We have the
following lemma for variances.

######  Lemma 5.4.

The variance of the semi-random APMC solution is bounded in the sense that

(73)
$\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n}-\tilde{v}_{i}^{n}\right]\right]\Delta
x\leq\frac{\displaystyle p}{1+\displaystyle p}\mathcal{E}_{0},\quad\forall n,$

where
$\mathcal{E}_{0}=\sum_{i}\left[(a-b^{2})\left(u_{i}^{0}\right)^{2}+\left(bu_{i}^{0}-v_{i}^{0}\right)^{2}\right]\Delta
x$. Namely, for any $\varepsilon\ll 1$, the variance is uniformly bounded and
tends to zero as $\varepsilon\rightarrow 0$.

######  Remark 5.5.

We can use the same trick as in the proof of Lemma 5.3 to get an estimate of
$\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]$,

(74)
$\displaystyle\sum_{j=1}^{n+1}\sum_{i}\mathrm{Var}\left[bu_{i}^{j}-v_{i}^{j}\right]\Delta
x\leq\frac{1}{\displaystyle
p}\sum_{j=1}^{n+1}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j}-v_{i}^{j}\right]\right)^{2}\Delta
x\leq\frac{\displaystyle p}{1-\displaystyle p^{2}}\mathcal{E}_{0},\quad\forall
n,$

in particular,

(75) $\sum_{i}\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]\Delta x\rightarrow
0,\quad\text{as}\quad n\rightarrow\infty.$

The main difference for the variance analysis of the fully-random APMC in the
following Lemma 5.6, is to use the conditional variance formula for
characteristic variables in the random convection step.

######  Lemma 5.6.

The variance of the fully-random APMC solutions is bounded in the sense that

(76)
$\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n}\right]+\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]\right]\Delta
x\leq\frac{\displaystyle p}{1+\displaystyle
p}\mathcal{E}_{0}+\frac{a(1-\nu)}{\nu}T\mathcal{C}_{0}\Delta t,\quad\forall
n,$

where

$\displaystyle\mathcal{E}_{0}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(u_{i}^{0}\right)^{2}+\left(bu_{i}^{0}-v_{i}^{0}\right)^{2}\right]\Delta
x,$ $\displaystyle\mathcal{C}_{0}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(\frac{u_{i+1}^{0}-u_{i}^{0}}{\Delta
x}\right)^{2}+\left[b\left(\frac{u_{i+1}^{0}-u_{i}^{0}}{\Delta
x}\right)-\left(\frac{v_{i+1}^{0}-v_{i}^{0}}{\Delta
x}\right)\right]^{2}\right]\Delta x.$

Then for any $\varepsilon\ll 1$, the variance is uniformly bounded and tends
to zero as $\varepsilon\rightarrow 0$, $\Delta t\rightarrow 0$. Therefore, the
APMC method still works in the regime when $\varepsilon\ll 1$ and $\Delta t$
independent of $\varepsilon$, i.e. has the asymptotic preserving property.

###  5.3 Multilevel APMC

We introduce the multilevel version of APMC (MLAPMC) algorithm here. Since the
random variables are distinct at the finest temporal and spatial
discretization level, we encounter similar situation as the Case II in Example
3, Example 4.1 and Example 4.2. Therefore, we need to use the trick introduced
in [24] and in (35) to construct coarser random variables and compute MLMC
corrections.

The fully-random MLAPMC algorithm of one sample on level $l$ can be stated as
follows.

Algorithm 4 Multilevel APMC (MLAPMC) - fully-random version

Given $T>0$, time step $\Delta t_{l}$, CFL number $\nu=\sqrt{a}\Delta
t_{l}/\Delta x_{l}$.

1: Initialize $u^{0}_{i}=u_{0}(x_{i}),v^{0}_{i}=v_{0}(x_{i})$,
$\bar{u}^{0}_{i_{1}}=u_{0}(x_{i_{1}\Gamma}),\bar{v}^{0}_{i_{1}}=v_{0}(x_{i_{1}\Gamma})$,
for $i_{1}=0,\dots,X_{l-1}-1$;

2: for $n=0,\dots,M_{l-1}-1$ do

3: for $k=0,\dots,\Gamma-1$ do

4: compute $\rho^{n^{\prime}}_{i},l^{n^{\prime}}_{i}$ as in (57), where
$n^{\prime}=n\Gamma+k$;

5: sample $\rho^{n^{\prime}+\frac{1}{2}}_{i}$ with
$\xi_{i,k}\sim\mathcal{U}(0,1)$ and $l^{n^{\prime}+\frac{1}{2}}_{i}$ with
$\eta_{i,k}\sim\mathcal{U}(0,1)$ as in (61);

6: compute $u^{n^{\prime}+\frac{1}{2}}_{i}$, $v^{n^{\prime}+\frac{1}{2}}_{i}$
as in (60);

7: $u^{n^{\prime}+1}_{i}=u^{n^{\prime}+\frac{1}{2}}_{i}$;

8: sample $v^{n^{\prime}+1}_{i}$ with $\zeta_{i,k}\sim\mathcal{U}(0,1)$ as in
(64);

9: end for

10: compute $\bar{\rho}^{n}_{i_{1}},\bar{l}^{n}_{i_{1}}$ as in (57);

11: $\xi_{i_{1}}=\left(\max_{i_{1}\leq i<i_{1}+\Gamma,0\leq
k<\Gamma}\\{\xi_{i,k}\\}\right)^{\Gamma^{2}}$,
$\eta_{i_{1}}=\left(\max_{i_{1}\leq i<i_{1}+\Gamma,0\leq
k<\Gamma}\\{\eta_{i,k}\\}\right)^{\Gamma^{2}}$;

12: sample $\bar{\rho}^{n+\frac{1}{2}}_{i_{1}}$ with $\xi_{i_{1}}$ and
$\bar{l}^{n+\frac{1}{2}}_{i_{1}}$ with $\eta_{i_{1}}$ as in (61);

13: compute $\bar{u}^{n+\frac{1}{2}}_{i_{1}}$,
$\bar{v}^{n+\frac{1}{2}}_{i_{1}}$ as in (60);

14: $\bar{u}^{n+1}_{i_{1}}=\bar{u}^{n+\frac{1}{2}}_{i_{1}}$;

15: $\zeta_{i_{1}}=\left(\max_{i_{1}\leq i<i_{1}+\Gamma,0\leq
k<\Gamma}\\{\zeta_{i,k}\\}\right)^{\Gamma^{2}}$;

16: sample $\bar{v}^{n+1}_{i_{1}}$ with $\zeta_{i_{1}}$ as in (64);

17: end for

As discussed before, the variance in the asymptotic regime ($\varepsilon\ll
1$) is small, and APMC performs well enough. For comparison, we apply MLAPMC
in the non-stiff regime ($\varepsilon=1$) with semi-random and fully-random
schemes respectively. We consider the initial value

(77) $u_{0}(x)=\frac{\sin(x)+1}{2},\quad v_{0}(x)=0,$

with $x\in[-1,1]$, $t\in[0,T]$, and periodic boundary condtion.

We take the numerical parameters

(78) $\nu=\sqrt{a}\frac{\Delta t_{l}}{\Delta
x_{l}}=\frac{1}{2},\quad\Gamma=2,\quad\Delta x_{0}=\frac{1}{32},\quad
T=1,\quad N_{I}=500,\quad a=1,\quad b=2.$

We implement both MLAPMC and APMC algorithms with accuracy thresholds
$\delta=0.02$, $0.01$, $0.005$, $0.002$, $0.001$. The results are shown in
Fig. 10 and Fig. 11. The costs of two algorithms with respect to $\delta$ are
plotted in the left subfigure, and the right subfigure shows the variances of
$P_{l}^{(l)}$ and $P_{l}^{(l)}-P_{l}^{(l-1)}$. We can see that both
$\mathrm{Var}\left[P_{l}^{(l)}\right]$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ decay as $\Delta t_{l}$
decreases in semi-random case, and the decay orders are about $0.5$, i.e.
$\beta_{0}=\beta=0.5$. The costs of APMC and MLAPMC are both
$O(\delta^{-3.5})$. While in the fully-random case,
$\mathrm{Var}\left[P_{l}^{(l)}\right]$ and
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$ decay very slowly, i.e.
$\beta\simeq 0$ and $\beta_{0}\simeq 0$. The costs of APMC and MLAPMC in
fully-random case are both $O(\delta^{-4})$.

Figure 10: The computational cost w.r.t accuracy threshold $\delta$ (left) and
variance on different levels (right) using semi-random APMC and MLAPMC
methods.

Figure 11: The computational cost w.r.t accuracy threshold $\delta$ (left) and
variance on different levels (right) using fully-random APMC and MLAPMC
methods.

##  6 Conclusion

In this paper, we study the Multilevel Monte Carlo method, and compare it with
Monte Carlo method for applications in linear and nonlinear hyperbolic systems
with randomness originated from either random parameters or randomized
algorithms. In many cases, MLMC has lower cost than MC. Nevertheless, if the
variance of MC solutions decays at the same rate as MLMC corrections, the
standard MC performs well and a multilevel method may not provide improvement.

We conclude through both analytic studies and a few numerical experiments with
different levels of randomness, that the possible cause is white noise like
random variables in the SPDE or random algorithm, which may influence the
decay rates for MLMC corrections and lead to the deterioration of the MLMC
complexity. If the white noise random variables depend only on time, the
variances have first order decay with respect to the discretization level or
time step, and both MLMC and MC have $O(\delta^{-3})$ computational cost.
While if the white noise random variables depend further on space, the
variances do not decay and both MLMC and MC have $O(\delta^{-4})$ cost. A
rigorous analysis of the scalar advection equation with a random velocity is
given and is consistent with our numerical results.

In summary, we have identified three regimes for the performances of MLMC and
MC for hyperbolic equations. More research is needed to investigate SPDEs
driven by Gaussian and Levy white noise [14, 30], and to develop acceleration
methods for SPDEs and random algorithms for interesting physical applications.

## References

  * [1] A. Abdulle and A. Blumenthal, Stabilized multilevel monte carlo method for stiff stochastic differential equations, Journal of Computational Physics, 251 (2013), pp. 445–460.
  * [2] R. Abgrall and S. Mishra, Uncertainty quantification for hyperbolic systems of conservation laws, in Handbook of Numerical Analysis, vol. 18, Elsevier, 2017, pp. 507–544.
  * [3] N. Ayi and E. Faou, Analysis of an asymptotic preserving scheme for stochastic linear kinetic equations in the diffusion limit, SIAM/ASA Journal on Uncertainty Quantification, 7 (2019), pp. 760–785.
  * [4] W. Bao and S. Jin, The random projection method for hyperbolic conservation laws with stiff reaction terms, Journal of Computational Physics, 163 (2000), pp. 216–248.
  * [5] A. Barth, C. Schwab, and N. Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients, Numerische Mathematik, 119 (2011), pp. 123–161.
  * [6] A. J. Chorin, Random choice solution of hyperbolic systems, Journal of computational physics, 22 (1976), pp. 517–533.
  * [7] A. J. Chorin, Random choice methods with applications to reacting gas flow, Journal of computational physics, 25 (1977), pp. 253–272.
  * [8] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients, Computing and Visualization in Science, 14 (2011), pp. 3–15.
  * [9] S. Dereich and F. Heidenreich, A multilevel monte carlo algorithm for lévy-driven stochastic differential equations, Stochastic Processes and their Applications, 121 (2011), pp. 1565–1587.
  * [10] G. Dimarco and L. Pareschi, Hybrid multiscale methods for hyperbolic problems i. hyperbolic relaxation problems, Communications in Mathematical Sciences, 4 (2006), pp. 155–177.
  * [11] G. Dimarco, L. Pareschi, and G. Samaey, Asymptotic-Preserving Monte Carlo methods for transport equations in the diffusive limit, SIAM Journal on Scientific Computing, 40 (2018), pp. A504–A528.
  * [12] M. B. Giles, Multilevel Monte Carlo path simulation, Operations research, 56 (2008), pp. 607–617.
  * [13] J. Glimm, Solutions in the large for nonlinear hyperbolic systems of equations, Communications on Pure and Applied Mathematics, 18 (1965), pp. 697–715.
  * [14] H. Holden, B. Øksendal, J. Ubøe, and T. Zhang, Stochastic partial differential equations, Birkhäuser Boston, Boston, MA, 1996, pp. 141–191.
  * [15] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM Journal on Scientific Computing, 21 (1999), pp. 441–454.
  * [16] S. Jin and L. Pareschi, Uncertainty quantification for hyperbolic and kinetic equations, vol. 14, Springer, 2018.
  * [17] S. Jin, L. Pareschi, and G. Toscani, Uniformly accurate diffusive relaxation schemes for multiscale transport equations, SIAM Journal on Numerical Analysis, 38 (2000), pp. 913–936.
  * [18] S. Jin and Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Communications on pure and applied mathematics, 48 (1995), pp. 235–276.
  * [19] R. Käppeli and S. Mishra, Well-balanced schemes for the euler equations with gravitation, Journal of Computational Physics, 259 (2014), pp. 199–219.
  * [20] K. Karhunen, Über lineare Methoden in der Wahrscheinlichkeitsrechnung, Ann. Acad. Sci. Fennicae Ser. A. I. Math.-Phys., 1947 (1947), p. 79.
  * [21] A. Kurganov, Finite-volume schemes for shallow-water equations, Acta Numerica, 27 (2018), pp. 289–351.
  * [22] D. Liu, Uncertainty quantification with shallow water equations, University of Florence, (2010).
  * [23] M. Loève, Fonctions aléatoires de second ordre, Revue Sci., 84 (1946), pp. 195–206.
  * [24] E. Løvbak, G. Samaey, and S. Vandewalle, A multilevel Monte Carlo method for asymptotic-preserving particle schemes in the diffusive limit, Numerische Mathematik, 148 (2021), pp. 141–186.
  * [25] S. Mishra, C. Schwab, and J. Šukys, Multi-level monte carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions, Journal of Computational Physics, 231 (2012), pp. 3365–3388.
  * [26] S. Mishra, C. Schwab, and J. Šukys, Multilevel monte carlo finite volume methods for shallow water equations with uncertain topography in multi-dimensions, SIAM Journal on Scientific Computing, 34 (2012), pp. B761–B784.
  * [27] H. N. Najm, Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics, Annual review of fluid mechanics, 41 (2009), pp. 35–52.
  * [28] R. Natalini, Recent mathematical results on hyperbolic relaxation problems, in Analysis of Systems of Conservation Laws, 1999, pp. 128–198.
  * [29] J. Tryoen, O. Le Maitre, M. Ndjinga, and A. Ern, Intrusive galerkin methods with upwinding for uncertain nonlinear hyperbolic systems, Journal of Computational Physics, 229 (2010), pp. 6485–6511.
  * [30] J. B. Walsh, An introduction to stochastic partial differential equations, in École d’Été de Probabilités de Saint Flour XIV \- 1984, P. L. Hennequin, ed., Berlin, Heidelberg, 1986, Springer Berlin Heidelberg, pp. 265–439.

##  Appendix A Some known results of MLMC for differential equations

We review some applications of MLMC in the literature. For consistency of
notations, we always denote $P_{l}^{(l)}$, $P_{l}^{(l-1)}$ as the numerical
solutions (at final time $T$ if time-dependent) using time step $\Delta t_{l}$
or $\Delta t_{l-1}$ and mesh size $\Delta x_{l}$ or $\Delta x_{l-1}$ on level
$l$ respectively.

  * •

Stochastic differential equations [1, 9, 12].

For example, we consider the forward Euler method with time step $\Delta t$ of
the SDE

$\mathrm{d}u(t)=a(u,t)\mathrm{d}t+b(u,t)\mathrm{d}W,$ (79)
$u_{n+1}=u_{n}+a(u_{n},t_{n})\Delta t+b(u_{n},t_{n})\Delta W_{n}.$

In this case, the random variables are Gaussian white noises. Due to the
additivity of Gaussian distribution, it is natural to obtain a random variable
on the coarser level by adding up the random variables on the finer level,
leading to the reduction of
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]$. As proved in [12], we
have $\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1)$,
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=O(\Delta t_{l})$, i.e.
$\beta_{0}=0$, $\beta=1$.

  * •

PDE with random coefficients [5, 8].

In subsurface flow applications, the hydraulic conductivity can be modelled as
a random field $k=k({\bf{x}},\omega)$ on $D\times\Omega$. We have the
following Darcy’s equation with random coefficients in the form

(80) $-\nabla\cdot(k({\bf{x}},\omega)\nabla
p({\bf{x}},\omega))=f({\bf{x}}),\quad\text{in }D.$

To apply MLMC, we need to sample from the input random field
$k({\bf{x}},\omega)$, e.g. by a truncation of the Karhunen-Loève (KL)
expansion [20, 23]. Then for each sample with fixed $\omega$, we solve the PDE
with multigrid methods. In this case, the dimension of random variables is
fixed by the truncation order of the KL expansion and same random variables
are used to calculate $P_{l}^{(l-1)}$, $P_{l}^{(l)}$. The numerical results in
[8] show that $\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1)$,
$\mathrm{Var}\left[P_{l}^{(l)}-p_{l}^{(l-1)}\right]=O((\Delta t)^{2})$, i.e.
$\beta_{0}=0$, $\beta=2$.

  * •

Particle method for transport equations [24]. For the BGK model

(81)
$\partial_{t}f(x,v,t)+\frac{v}{\varepsilon}\partial_{x}f(x,v,t)=\frac{1}{\varepsilon^{2}}\left(\mathcal{M}(v)\rho(x,t)-f(x,v,t)\right),$

where $f(x,v,t)$ is the particle density distribution at time $t$, position
$x$ with particle velocity $v$ and $\mathcal{M}(v)$ is the absolute
Maxwellian. $\varepsilon$ is the Knudsen number, the dimensionaless mean free
path. An asymptotic preserving particle scheme was proposed in [11]. Given the
current state of a particle in the position-velocity phase space, its position
is updated based on its velocity in the transport step. Then, the particle
changes its velocity with the probability $p_{\Delta t_{l}}=\frac{\Delta
t_{l}}{\varepsilon^{2}+\Delta t_{l}}$ in the collision step, such that the
collision is performed if $\xi_{l}\leq p_{\Delta t_{l}}$ with
$\xi_{l}\sim\mathcal{U}(0,1)$.

To apply MLMC, Løvbak, Samaey, and Vandewalle [24] proposed to use random
variables $\left\\{\xi_{l,k}\right\\}_{k=1}^{\Gamma}$ in the finer level to
construct coarser level random varibales,

(82) $\xi_{l-1}=\left(\max_{1\leq
k\leq\Gamma}\xi_{l,k}\right)^{\Gamma}\sim\mathcal{U}(0,1).$

If $\xi_{l-1}\leq p_{\Delta t_{l-1}}$, collision is performed on the coarser
level. Therefore the trajectories of particles in the finer level and the
coarser level are correlated. In this case, particles are independent of each
other, similar to the SDE case.
$\mathrm{Var}\left[P_{l}^{(l)}-P_{l}^{(l-1)}\right]=O(\Delta t)$, i.e.
$\beta=1$ is proved in [24, Lemma 9]. The numerical results show that
$\mathrm{Var}\left[P_{l}^{(l)}\right]=O(1)$, i.e. $\beta_{0}=0$.

##  Appendix B Analysis of MLMC and MC in § 2.2

###  B.1 Proof of Eq. 17, cost of MLMC

######  Proof B.1.

Let
$L=\left\lceil\frac{\log(\sqrt{2}c_{1}T^{\alpha}\delta^{-1})}{\alpha\log\Gamma}\right\rceil$,
the bias term $c_{1}(\Delta t_{L})^{\alpha}$ is controlled from above and
below by

(83) $\frac{1}{\sqrt{2}}\Gamma^{-\alpha}\delta<c_{1}(\Delta
t_{L})^{\alpha}\leq\frac{1}{\sqrt{2}}\delta,$

where $\Gamma=\Delta t_{l-1}/\Delta t_{l}$. Equivalently, we have

(84)
$\left|\mathrm{E}\left[\hat{P}^{(L)}-P\right]\right|^{2}=\left|\mathrm{E}\left[P^{(L)}-P\right]\right|^{2}\leq\frac{1}{2}\delta^{2}.$

According to (5), the variance of the MLMC estimate $\hat{P}^{(L)}$ is

(85)
$\mathrm{Var}\left[\hat{P}^{(L)}\right]=\sum_{l=0}^{L}N_{l}^{-1}V_{l}=\mu^{-1}\sum_{l=0}^{L}\sqrt{C_{l}V_{l}}=\frac{1}{2}\delta^{2},$

which implies $\mu=2\delta^{-2}\sum_{l=0}^{L}\sqrt{C_{l}V_{l}}$, and
$N_{l}=\left\lceil\mu\sqrt{V_{l}/C_{l}}\right\rceil$, $l=0,\dots,L$. The total
computational cost is

(86)
$\mathfrak{C}=\sum_{l=0}^{L}N_{l}C_{l}\leq\sum_{l=0}^{L}\left(\mu\sqrt{C_{l}V_{l}}+C_{l}\right)=2\delta^{-2}\left(\sum_{l=0}^{L}\sqrt{C_{l}V_{l}}\right)^{2}+\sum_{l=0}^{L}C_{l}.$

By properties A2, A3, we have

(87)
$\sum_{l=0}^{L}\sqrt{C_{l}V_{l}}\leq\sqrt{c_{2}c_{3}}\sum_{l=0}^{L}(\Delta
t_{l})^{\frac{\beta-\gamma}{2}}=\left\\{\begin{aligned}
&c_{5},&\beta>\gamma,\\\ &L+1\leq c_{6}\log\delta^{-1},&\beta=\gamma,\\\
&c_{7}(\Delta t_{L})^{\frac{\beta-\gamma}{2}}\leq
c_{8}\delta^{-(\gamma-\beta)/(2\alpha)},&0<\beta<\gamma,\end{aligned}\right.$

and

(88) $\sum_{l=0}^{L}C_{l}\leq c_{3}\sum_{l=0}^{L}(\Delta t_{l})^{-\gamma}\leq
c_{9}\delta^{-\gamma/\alpha}\leq c_{9}\delta^{-2}.$

In summary, the computational cost $\mathfrak{C}$ is bounded by

(89) $\displaystyle\mathfrak{C}$ $\displaystyle\leq\left\\{\begin{aligned}
&c_{4}\delta^{-2},&\beta>\gamma,\\\
&c_{4}\delta^{-2}\left(\log\delta^{-1}\right)^{2},&\beta=\gamma,\\\
&c_{4}\delta^{-2-(\gamma-\beta)/\alpha},&0<\beta<\gamma.\end{aligned}\right.$

###  B.2 Proof of Eq. 19, cost of MC

######  Proof B.2.

We can perform similar analysis for the adaptive MC given in Algorithm 2. The
difference is that, for MLMC the bound for the total variance in (85) is
needed, while for MC, we only need the bound for the variance on each level,
as $\mathrm{Var}\left[\hat{P}_{l}^{(l)}\right]$ in (92).

We have similar estimate for the bias by taking
$L=\left\lceil\frac{\log(\sqrt{2}c_{1}T^{\alpha}\delta^{-1})}{\alpha\log\Gamma}\right\rceil$,

(90) $\frac{1}{\sqrt{2}}\Gamma^{-\alpha}\delta<c_{1}(\Delta
t_{L})^{\alpha}\leq\frac{1}{\sqrt{2}}\delta,$

where $\Gamma=\Delta t_{l-1}/\Delta t_{l}$, and equivalently,

(91)
$\left|\mathrm{E}\left[\hat{P}_{L}^{(L)}-P\right]\right|^{2}=\left|\mathrm{E}\left[P^{(L)}-P\right]\right|^{2}\leq\frac{1}{2}\delta^{2}.$

The variances of each level should be bounded by

(92)
$\mathrm{Var}\left[\hat{P}_{l}^{(l)}\right]=N_{l}^{-1}\mathsf{V}_{l}\leq\frac{1}{2}\delta^{2},\quad\Rightarrow\quad
N_{l}=\left\lceil 2\delta^{-2}\mathsf{V}_{l}\right\rceil,\quad l=0,\dots,L.$

to make sure that $\mathrm{MSE}<\delta^{2}$. The total computational cost
$\mathfrak{C}^{*}$ is

(93)
$\mathfrak{C}^{*}=\sum_{l=0}^{L}N_{l}\mathfrak{C}\left[P_{l}^{(l)}\right]\leq
2\delta^{-2}\sum_{l=0}^{L}\mathfrak{C}\left[P_{l}^{(l)}\right]\mathsf{V}_{l}+\sum_{l=0}^{L}\mathfrak{C}\left[P_{l}^{(l)}\right].$

With properties A5, A6, we have

(94) $\sum_{l=0}^{L}\mathfrak{C}\left[P_{l}^{(l)}\right]\mathsf{V}_{l}\leq
c_{2}c_{3}\sum_{l=0}^{L}(\Delta
t_{l})^{\beta_{0}-\gamma}=\left\\{\begin{aligned} &c_{5},&\beta_{0}>\gamma,\\\
&L+1\leq c_{6}\log\delta^{-1},&\beta_{0}=\gamma,\\\ &c_{7}(\Delta
t_{L})^{\beta_{0}-\gamma}\leq
c_{8}\delta^{-(\gamma-\beta_{0})/\alpha},&0<\beta_{0}<\gamma,\end{aligned}\right.$

and

(95) $\sum_{l=0}^{L}\mathfrak{C}\left[P_{l}^{(l)}\right]\leq
c_{3}\sum_{l=0}^{L}(\Delta t_{l})^{-\gamma}\leq
c_{9}\delta^{-\gamma/\alpha}\leq c_{9}\delta^{-2}.$

In summary, the computational cost $\mathfrak{C}^{*}$ is bounded by

(96) $\displaystyle\mathfrak{C}^{*}$ $\displaystyle\leq\left\\{\begin{aligned}
&c_{4}\delta^{-2},&\beta_{0}>\gamma,\\\
&c_{4}\delta^{-2}\log\delta^{-1},&\beta_{0}=\gamma,\\\
&c_{4}\delta^{-2-(\gamma-\beta_{0})/\alpha},&0<\beta_{0}<\gamma.\end{aligned}\right.$

##  Appendix C AP property for the deterministic method of the Jin-Xin model

For completeness, we formulate the AP property of the deterministic numerical
method in § 5.1 in the following proposition.

######  Proposition C.1.

The deterministic numerical method in § 5.1 is asymptotic preserving since it
has the correct numerical $\varepsilon\to 0$ limit, namely, in the limit
$\varepsilon\rightarrow 0$, the numerical solutions $u_{i}^{n}$, $v_{i}^{n}$
tends to the numerical solutions of the limit equations

(97) $v=bu,\quad\frac{\partial u}{\partial t}+b\frac{\partial u}{\partial
x}=0.$

We provide a heuristic argument in the following, for complete theory, please
refer to [15]. As $\varepsilon\rightarrow 0$, the implicit Euler’s method used
in relaxation step collapses to

(98) $v_{i}^{n+1}=bu_{i}^{n+1},\quad n=0,\dots,M_{l}-1.$

Applying (98) into the convection step gives

(99) $\rho_{i}^{n}=(\sqrt{a}+b)u_{i}^{n},\quad
l_{i}^{n}=(\sqrt{a}-b)u_{i}^{n},$

and

(100)
$u_{i}^{n+1}=(1-\nu)u_{i}^{n}+\frac{(\sqrt{a}+b)\nu}{2\sqrt{a}}u_{i-1}^{n}+\frac{(\sqrt{a}-b)\nu}{2\sqrt{a}}u_{i+1}^{n},$

which is equivalent to

(101) $\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta
t}+b\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2\Delta x}=\frac{\sqrt{a}\Delta
x}{2}\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{(\Delta x)^{2}}.$

(98), (101) is a first-order consistent and stable discretization of (97),
with leading truncation error

(102) $\frac{\sqrt{a}\Delta x}{2}\frac{\partial^{2}u}{\partial x^{2}}.$

Thus this scheme has the correct asymptotic limit and is expected to capture
the correct limiting behavior in the regime when $\varepsilon\ll 1$ and
$\Delta t\gg\varepsilon$, with a first-order accuracy.

##  Appendix D Analysis for the APMC method in § 5.2

###  D.1 Proof of Lemma 5.1

######  Proof D.1.

Combining (59) and (60), we obtain the formula of
$\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]$,
$\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]$,

(103) $\displaystyle\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]$
$\displaystyle=(1-\nu)\mathrm{E}\left[v_{i}^{n}\right]+\frac{\nu}{2}\left(\mathrm{E}\left[v_{i+1}^{n}\right]+\mathrm{E}\left[v_{i-1}^{n}\right]\right)-\frac{\nu\sqrt{a}}{2}\left(\mathrm{E}\left[u_{i+1}^{n}\right]-\mathrm{E}\left[u_{i-1}^{n}\right]\right),$
$\displaystyle\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]$
$\displaystyle=(1-\nu)\mathrm{E}\left[u_{i}^{n}\right]+\frac{\nu}{2}\left(\mathrm{E}\left[u_{i+1}^{n}\right]+\mathrm{E}\left[u_{i-1}^{n}\right]\right)-\frac{\nu}{2\sqrt{a}}\left(\mathrm{E}\left[v_{i+1}^{n}\right]-\mathrm{E}\left[v_{i-1}^{n}\right]\right).$

We compute the following quadratic sums and simplify them as follows

(104)
$\displaystyle\quad\sum_{i}\left[a\left(\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\right)^{2}+\left(\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]\right)^{2}\right]=\frac{1}{2}\sum_{i}\left[\left(\mathrm{E}\left[\rho_{i}^{n+\frac{1}{2}}\right]\right)^{2}+\left(\mathrm{E}\left[l_{i}^{n+\frac{1}{2}}\right]\right)^{2}\right]$
$\displaystyle=\frac{(1-\nu)^{2}+\nu^{2}}{2}\sum_{i}\left[\left(\mathrm{E}\left[\rho_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[l_{i}^{n}\right]\right)^{2}\right]+\nu(1-\nu)\sum_{i}\left[\mathrm{E}\left[\rho_{i}^{n}\right]\mathrm{E}\left[\rho_{i+1}^{n}\right]+\mathrm{E}\left[l_{i}^{n}\right]\mathrm{E}\left[l_{i+1}^{n}\right]\right]$
$\displaystyle=\left[(1-\nu)^{2}+\nu^{2}\right]\sum_{i}\left[a\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[v_{i}^{n}\right]\right)^{2}\right]+2\nu(1-\nu)\sum_{i}\left[a\mathrm{E}\left[u_{i}^{n}\right]\mathrm{E}\left[u_{i+1}^{n}\right]+\mathrm{E}\left[v_{i}^{n}\right]\mathrm{E}\left[v_{i+1}^{n}\right]\right],$
(105)
$\displaystyle\quad\sum_{i}\left[\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]\right]=\frac{1}{4\sqrt{a}}\sum_{i}\left[\left(\mathrm{E}\left[\rho_{i}^{n+\frac{1}{2}}\right]\right)^{2}-\left(\mathrm{E}\left[l_{i}^{n+\frac{1}{2}}\right]\right)^{2}\right]$
$\displaystyle=\frac{(1-\nu)^{2}+\nu^{2}}{4\sqrt{a}}\sum_{i}\left[\left(\mathrm{E}\left[\rho_{i}^{n}\right]\right)^{2}-\left(\mathrm{E}\left[l_{i}^{n}\right]\right)^{2}\right]+\frac{\nu(1-\nu)}{2\sqrt{a}}\sum_{i}\left[\mathrm{E}\left[\rho_{i}^{n}\right]\mathrm{E}\left[\rho_{i+1}^{n}\right]-\mathrm{E}\left[l_{i}^{n}\right]\mathrm{E}\left[l_{i+1}^{n}\right]\right]$
$\displaystyle=\left[(1-\nu)^{2}+\nu^{2}\right]\sum_{i}\left[\mathrm{E}\left[u_{i}^{n}\right]\mathrm{E}\left[v_{i}^{n}\right]\right]+\nu(1-\nu)\sum_{i}\left[\mathrm{E}\left[u_{i}^{n}\right]\mathrm{E}\left[v_{i+1}^{n}\right]+\mathrm{E}\left[v_{i}^{n}\right]\mathrm{E}\left[u_{i+1}^{n}\right]\right],$

which lead to the following estimate

(106) $\displaystyle\mathcal{E}_{n+\frac{1}{2}}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\right)^{2}+\left(\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right)^{2}\right]\Delta
x$
$\displaystyle=\sum_{i}\left[a\left(\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\right)^{2}+\left(\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]\right)^{2}-2b\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\mathrm{E}\left[v_{i}^{n+\frac{1}{2}}\right]\right]\Delta
x$
$\displaystyle=[(1-\nu)^{2}+\nu^{2}]\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\right)^{2}\right]\Delta
x$
$\displaystyle\quad+2\nu(1-\nu)\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n}\right]\mathrm{E}\left[u_{i+1}^{n}\right]\right)+\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\mathrm{E}\left[bu_{i+1}^{n}-v_{i+1}^{n}\right]\right]\Delta
x$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\right)^{2}\right]\Delta
x=\mathcal{E}_{n}.$

In the relaxation step, by (63), we have

(107)
$\mathrm{E}\left[u_{i}^{n+1}\right]=\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right],\quad\mathrm{E}\left[bu_{i}^{n+1}-v_{i}^{n+1}\right]=\displaystyle
p\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right],$

thus

(108) $\displaystyle\mathcal{E}_{n+1}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i}^{n+\frac{1}{2}}\right]\right)^{2}+p^{2}\left(\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right)^{2}\right]\Delta
x\leq\mathcal{E}_{n+\frac{1}{2}}\leq\mathcal{E}_{n}.$

Therefore the energy decreases monotonically. And (68) can be easily derived
from (67). The boundedness of
$\left\|\mathrm{E}\left[v^{n}\right]\right\|_{2}$ (69) follows by the fact
that

(109) $\sum_{i}\left(\mathrm{E}\left[v_{i}^{n}\right]\right)^{2}\Delta x\leq
2\sum_{i}\left[\left(\mathrm{E}\left[bu_{i}^{n}-v_{i}^{n}\right]\right)^{2}+b^{2}\left(\mathrm{E}\left[u_{i}^{n}\right]\right)^{2}\right]\Delta
x\leq\frac{2a}{a-b^{2}}\mathcal{E}_{n}.$

###  D.2 Proof of Lemma 5.3

######  Proof D.2.

(108) implies that

(110)
$\displaystyle\mathcal{E}_{n+1}+\sum_{i}(1-p^{2})\left(\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right)^{2}\Delta
x=\mathcal{E}_{n+\frac{1}{2}}\leq\mathcal{E}_{n}.$

Therefore

(111)
$\displaystyle\quad\mathcal{E}_{n+1}+\sum_{j=0}^{n}\sum_{i}(1-p^{2})\left(\mathrm{E}\left[bu_{i}^{j+\frac{1}{2}}-v_{i}^{j+\frac{1}{2}}\right]\right)^{2}\Delta
x$
$\displaystyle\leq\mathcal{E}_{n}+\sum_{j=0}^{n-1}\sum_{i}(1-p^{2})\left(\mathrm{E}\left[bu_{i}^{j+\frac{1}{2}}-v_{i}^{j+\frac{1}{2}}\right]\right)^{2}\Delta
x$ $\displaystyle\leq\cdots$
$\displaystyle\leq\mathcal{E}_{1}+\sum_{i}(1-p^{2})\left(\mathrm{E}\left[bu_{i}^{\frac{1}{2}}-v_{i}^{\frac{1}{2}}\right]\right)^{2}\Delta
x$ $\displaystyle\leq\mathcal{E}_{0}.$

Noting that $\mathrm{E}\left[bu_{i}^{n+1}-v_{i}^{n+1}\right]=\displaystyle
p\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]$, we have

(112)
$\sum_{j=1}^{n+1}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j}-v_{i}^{j}\right]\right)^{2}\Delta
x=\displaystyle
p^{2}\sum_{j=0}^{n}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j+\frac{1}{2}}-v_{i}^{j+\frac{1}{2}}\right]\right)^{2}\Delta
x\leq\frac{\displaystyle p^{2}}{1-\displaystyle
p^{2}}\mathcal{E}_{0},\quad\forall n,$

namely, (72) holds. Therefore, as $\varepsilon\rightarrow 0$, $\displaystyle
p\rightarrow 0$, $\mathrm{E}\left[bu^{n}-v^{n}\right]$ converges to zero in
$\ell^{2}$.

###  D.3 Proof of Lemma 5.4, for the semi-random APMC

######  Proof D.3.

We have the following identities, for the change of variances in the
convection step,

(113)
$\displaystyle\quad\sum_{i}\left[a\mathrm{Var}\left[\tilde{u}_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[\tilde{v}_{i}^{n+\frac{1}{2}}\right]\right]=\frac{1}{2}\sum_{i}\left[\mathrm{Var}\left[\tilde{\rho}_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[\tilde{l}_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\left[(1-\nu)^{2}+\nu^{2}\right]\sum_{i}\left[a\mathrm{Var}\left[\tilde{u}_{i}^{n}\right]+\mathrm{Var}\left[\tilde{v}_{i}^{n}\right]\right]+2\nu(1-\nu)\sum_{i}a\mathrm{Cov}\left[\tilde{u}_{i}^{n},\tilde{u}_{i+1}^{n}\right]+\mathrm{Cov}\left[\tilde{v}_{i}^{n},\tilde{v}_{i+1}^{n}\right],$
$\displaystyle\quad\sum_{i}\mathrm{Cov}\left[\tilde{u}_{i}^{n+\frac{1}{2}},\tilde{v}_{i}^{n+\frac{1}{2}}\right]=\frac{1}{4\sqrt{a}}\sum_{i}\left[\mathrm{Var}\left[\tilde{\rho}_{i}^{n+\frac{1}{2}}\right]-\mathrm{Var}\left[\tilde{l}_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\left[(1-\nu)^{2}+\nu^{2}\right]\sum_{i}\mathrm{Cov}\left[\tilde{u}_{i}^{n},\tilde{v}_{i}^{n}\right]+\nu(1-\nu)\sum_{i}\left[\mathrm{Cov}\left[\tilde{u}_{i}^{n},\tilde{v}_{i+1}^{n}\right]+\mathrm{Cov}\left[\tilde{v}_{i}^{n},\tilde{u}_{i+1}^{n}\right]\right].$

Combining above identities, we have

(114)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n+\frac{1}{2}}-\tilde{v}_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=[(1-\nu)^{2}+\nu^{2}]\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n}-\tilde{v}_{i}^{n}\right]\right]$
$\displaystyle\quad+2\nu(1-\nu)\sum_{i}\left[(a-b^{2})\mathrm{Cov}\left[\tilde{u}_{i}^{n},\tilde{u}_{i+1}^{n}\right]+\mathrm{Cov}\left[b\tilde{u}_{i}^{n}-\tilde{v}_{i}^{n},b\tilde{u}_{i+1}^{n}-\tilde{v}_{i+1}^{n}\right]\right]$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n}-\tilde{v}_{i}^{n}\right]\right].$

We then treat the relaxation step, using the conditional variance formula

(115) $\displaystyle\mathrm{Var}\left[\tilde{u}_{i}^{n+1}\right]$
$\displaystyle=\mathrm{Var}\left[\tilde{u}_{i}^{n+\frac{1}{2}}\right],$
$\displaystyle\mathrm{Var}\left[b\tilde{u}_{i}^{n+1}-\tilde{v}_{i}^{n+1}\right]$
$\displaystyle=p\mathrm{Var}\left[b\tilde{u}_{i}^{n+\frac{1}{2}}-\tilde{v}_{i}^{n+\frac{1}{2}}\right]+p(1-p)\left(\mathrm{E}\left[b\tilde{u}_{i}^{n+\frac{1}{2}}-\tilde{v}_{i}^{n+\frac{1}{2}}\right]\right)^{2}.$

By substituting (114) into (115), we have

(116)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n+1}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n+1}-\tilde{v}_{i}^{n+1}\right]\right]$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[\tilde{u}_{i}^{n}\right]+\mathrm{Var}\left[b\tilde{u}_{i}^{n}-\tilde{v}_{i}^{n}\right]\right]+p(1-p)\sum_{i}\left(\mathrm{E}\left[b\tilde{u}_{i}^{n+\frac{1}{2}}-\tilde{v}_{i}^{n+\frac{1}{2}}\right]\right)^{2}$
$\displaystyle\leq\cdots$ $\displaystyle\leq
0+p(1-p)\sum_{j=0}^{n}\sum_{i}\left(\mathrm{E}\left[b\tilde{u}_{i}^{j+\frac{1}{2}}-\tilde{v}_{i}^{j+\frac{1}{2}}\right]\right)^{2}$
$\displaystyle=\frac{1-\displaystyle p}{\displaystyle
p}\sum_{j=0}^{n}\sum_{i}\left(\mathrm{E}\left[b\tilde{u}_{i}^{j+1}-\tilde{v}_{i}^{j+1}\right]\right)^{2}.$

(73) follows by a combination of (116) and Lemma 5.3.

###  D.4 Proof of Lemma 5.6, for the fully-random APMC

######  Proof D.4.

We have the similar identities for the change of the variances, when we
compute $u_{i}^{n+\frac{1}{2}},v_{i}^{n+\frac{1}{2}}$ from the characteristic
variables $\rho_{i}^{n+\frac{1}{2}}$, $l_{i}^{n+\frac{1}{2}}$,

(117)
$\displaystyle\sum_{i}\left[a\mathrm{Var}\left[u_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[v_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\sum_{i}\left[\mathrm{Var}\left[\frac{\rho_{i}^{n+\frac{1}{2}}+l_{i}^{n+\frac{1}{2}}}{2}\right]+\mathrm{Var}\left[\frac{\rho_{i}^{n+\frac{1}{2}}-l_{i}^{n+\frac{1}{2}}}{2}\right]\right]$
$\displaystyle=\frac{1}{2}\sum_{i}\left[\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]\right],$
$\displaystyle\sum_{i}\mathrm{Cov}\left[u_{i}^{n+\frac{1}{2}},v_{i}^{n+\frac{1}{2}}\right]$
$\displaystyle=\sum_{i}\mathrm{Cov}\left[\frac{\rho_{i}^{n+\frac{1}{2}}+l_{i}^{n+\frac{1}{2}}}{2\sqrt{a}},\frac{\rho_{i}^{n+\frac{1}{2}}-l_{i}^{n+\frac{1}{2}}}{2}\right]$
$\displaystyle=\frac{1}{4\sqrt{a}}\sum_{i}\left[\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]-\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]\right].$

Combining the identities in (117), we have

(118)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\frac{1}{2}\sum_{i}\left[\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]\right]-\frac{b}{2\sqrt{a}}\sum_{i}\left[\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]-\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\frac{1}{2}\sum_{i}\left[\left(1-\frac{b}{\sqrt{a}}\right)\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]+\left(1+\frac{b}{\sqrt{a}}\right)\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]\right].$

We apply the conditional variance formula

(119)
$\displaystyle\mathrm{Var}\left[\rho_{i}^{n+\frac{1}{2}}\right]=(1-\nu)\mathrm{Var}\left[\rho_{i}^{n}\right]+\nu\mathrm{Var}\left[\rho_{i-1}^{n}\right]+\nu(1-\nu)\left(\mathrm{E}\left[\rho_{i}^{n}-\rho_{i-1}^{n}\right]\right)^{2},$
$\displaystyle\mathrm{Var}\left[l_{i}^{n+\frac{1}{2}}\right]=(1-\nu)\mathrm{Var}\left[l_{i}^{n}\right]+\nu\mathrm{Var}\left[l_{i+1}^{n}\right]+\nu(1-\nu)\left(\mathrm{E}\left[l_{i+1}^{n}-l_{i}^{n}\right]\right)^{2},$

and write

(120)
$\displaystyle\mathrm{Var}\left[\rho_{i}^{n}\right]=\mathrm{Var}\left[\sqrt{a}u_{i}^{n}\right]+\mathrm{Var}\left[v_{i}^{n}\right]+2\sqrt{a}\mathrm{Cov}\left[u_{i}^{n},v_{i}^{n}\right],$
$\displaystyle\mathrm{Var}\left[l_{i}^{n}\right]=\mathrm{Var}\left[\sqrt{a}u_{i}^{n}\right]+\mathrm{Var}\left[v_{i}^{n}\right]-2\sqrt{a}\mathrm{Cov}\left[u_{i}^{n},v_{i}^{n}\right].$

By substituting (119), (120) into (118), it leads to the following identity

(121)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\sum_{i}\left[a\mathrm{Var}\left[u_{i}^{n}\right]+\mathrm{Var}\left[v_{i}^{n}\right]-2b\mathrm{Cov}\left[u_{i}^{n},v_{i}^{n}\right]\right]$
$\displaystyle\quad+\frac{\nu(1-\nu)}{2}\left[\left(1-\frac{b}{\sqrt{a}}\right)\sum_{i}\left(\mathrm{E}\left[\rho_{i}^{n}-\rho_{i-1}^{n}\right]\right)^{2}+\left(1+\frac{b}{\sqrt{a}}\right)\sum_{i}\left(\mathrm{E}\left[l_{i+1}^{n}-l_{i}^{n}\right]\right)^{2}\right].$

We then derive

(122)
$\displaystyle\quad\frac{1}{2}\left[\left(1-\frac{b}{\sqrt{a}}\right)\sum_{i}\left(\mathrm{E}\left[\rho_{i}^{n}-\rho_{i-1}^{n}\right]\right)^{2}+\left(1+\frac{b}{\sqrt{a}}\right)\sum_{i}\left(\mathrm{E}\left[l_{i+1}^{n}-l_{i}^{n}\right]\right)^{2}\right]$
$\displaystyle=\sum_{i}\left[a\left(\mathrm{E}\left[u_{i+1}^{n}-u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[v_{i+1}^{n}-v_{i}^{n}\right]\right)^{2}-2b\mathrm{E}\left[u_{i+1}^{n}-u_{i}^{n}\right]\mathrm{E}\left[v_{i+1}^{n}-v_{i}^{n}\right]\right]$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(\mathrm{E}\left[u_{i+1}^{n}-u_{i}^{n}\right]\right)^{2}+\left(\mathrm{E}\left[b\left(u_{i+1}^{n}-u_{i}^{n+1}\right)-\left(v_{i+1}^{n}-v_{i}^{n+1}\right)\right]\right)^{2}\right]:=V_{n},$

by noticing that

(123)
$\displaystyle\left(\mathrm{E}\left[\rho_{i}^{n}-\rho_{i-1}^{n}\right]\right)^{2}$
$\displaystyle=\left(\mathrm{E}\left[\sqrt{a}u_{i}^{n}-\sqrt{a}u_{i-1}^{n}+v_{i}^{n}-v_{i-1}^{n}\right]\right)^{2},$
$\displaystyle\left(\mathrm{E}\left[l_{i+1}^{n}-l_{i}^{n}\right]\right)^{2}$
$\displaystyle=\left(\mathrm{E}\left[\sqrt{a}u_{i+1}^{n}-\sqrt{a}u_{i}^{n}-v_{i+1}^{n}-v_{i}^{n}\right]\right)^{2}.$

In the convective part, we have

(124)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right]$
$\displaystyle=\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n}\right]+\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]\right]+\nu(1-\nu)V_{n},$

and also $V_{n}\Delta x\leq\mathcal{C}_{0}(\Delta x)^{2}$ by (70).

The analysis of the relaxation part is similar to that of (115), namely,

(125)
$\displaystyle\quad\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n+1}\right]+\mathrm{Var}\left[bu_{i}^{n+1}-v_{i}^{n+1}\right]\right]$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n+\frac{1}{2}}\right]+\mathrm{Var}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right]+p(1-p)\sum_{i}\left(\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right)^{2}$
$\displaystyle\leq\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n}\right]+\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]\right]+p(1-p)\sum_{i}\left(\mathrm{E}\left[bu_{i}^{n+\frac{1}{2}}-v_{i}^{n+\frac{1}{2}}\right]\right)^{2}+\nu(1-\nu)V_{n}$
$\displaystyle\leq\cdots$ $\displaystyle\leq\frac{(1-\displaystyle
p)}{\displaystyle
p}\sum_{j=1}^{n+1}\sum_{i}\left(\mathrm{E}\left[bu_{i}^{j}-v_{i}^{j}\right]\right)^{2}+\nu(1-\nu)\sum_{j=0}^{n}V_{n}.$

Thus

(126)
$\displaystyle\sum_{i}\left[(a-b^{2})\mathrm{Var}\left[u_{i}^{n}\right]+\mathrm{Var}\left[bu_{i}^{n}-v_{i}^{n}\right]\right]\Delta
x$ $\displaystyle\leq\frac{\displaystyle p}{1+\displaystyle
p}\mathcal{E}_{0}+\nu(1-\nu)n\mathcal{C}_{0}(\Delta x)^{2}$
$\displaystyle\leq\frac{\displaystyle p}{1+\displaystyle
p}\mathcal{E}_{0}+\frac{a(1-\nu)}{\nu}T\mathcal{C}_{0}\Delta t,\quad\forall
n,$

where

$\displaystyle\mathcal{E}_{0}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(u_{i}^{0}\right)^{2}+\left(bu_{i}^{0}-v_{i}^{0}\right)^{2}\right]\Delta
x,$ $\displaystyle\mathcal{C}_{0}$
$\displaystyle=\sum_{i}\left[(a-b^{2})\left(\frac{u_{i+1}^{0}-u_{i}^{0}}{\Delta
x}\right)^{2}+\left[b\left(\frac{u_{i+1}^{0}-u_{i}^{0}}{\Delta
x}\right)-\left(\frac{v_{i+1}^{0}-v_{i}^{0}}{\Delta
x}\right)\right]^{2}\right]\Delta x.$