# Monte Carlo Finite Volume Element Methods for the Convection-Diffusion Equation with a Random Diffusion Coefficient

Authors: Qian Zhang; Zhiyue Zhang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The paper presents a framework for the construction of Monte Carlo finite volume element method (MCFVEM) for the convection-diffusion equation with a random diffusion coefficient, which is described as a random field. We first approximate the continuous stochastic field by a finite number of random variables via the Karhunen-Loève expansion and transform the initial stochastic problem into a deterministic one with a parameter in high dimensions. Then we generate independent identically distributed approximations of the solution by sampling the coefficient of the equation and employing finite volume element variational formulation. Finally the Monte Carlo (MC) method is used to compute corresponding sample averages. Statistic error is estimated analytically and experimentally. A quasi-Monte Carlo (QMC) technique with Sobol sequences is also used to accelerate convergence, and experiments indicate that it can improve the efficiency of the Monte Carlo method.


---

## Body

## 1. Introduction

Mathematical models and computer simulations are widely used in engineering and science. To date, most attention is paid to deterministic mathematical models, where all input data are assumed to be perfectly known. However, in many cases, the parameters in the model are effected by uncertainty, either because they are not perfectly known or because they are intrinsically variable. The uncertainties are usually described as random variables or random fields. There is an increasing interest in including uncertainty in these models and quantifying its effect on the predict quantities of interest for applications, and there have been many numerical methods for partial differential equation whose input data (coefficients, forcing term, boundary conditions, and domain boundary) are uncertain.When the size of noise is relatively small, a perturbation method may be the most popular nonstatistical method, where the random field is expanded via Taylor series around its mean and truncated at certain order. This approach has been used extensively in various engineering fields [1–4]. Another approach is the Neumann expansion, which is based on the inverse of the stochastic operator in a Neumann series [5, 6]. Their applicability is often strongly dependent on the underlying operator and is typically limited to static problems.Much attention has recently been devoted to the development of Stochastic Galerkin (SG) method. Ghanem and Spanos pioneered a polynomial chaos expansion method [7], and its generalization, generalized polynomial chaos (gPC), has become one of the most widely used methods. With gPC, stochastic solutions are expressed as orthogonal polynomials of the input random parameters [8, 9]. The terms of gPC expansion depend on both the dimensionality of the random variable and the order of orthogonal polynomial interpolation. It performs exponential convergence. However gPC Galerkin method is more cumbersome to implement because the resulting equations for the expansion coefficients are almost coupled, especially for the problem that takes highly complex and nonlinear form. Another disadvantage of gPC Galerkin method appears when the dimension of random space is relatively large and the convergence rate heavily deteriorates, which is called curse of dimensionality.One of the most commonly used methods is Monte Carlo (MC) method [10] or one of its variants. MC methods are both general and simple to code, and they are naturally suited for parallelization. They generate a set of independent identically distributed approximations of the solution by sampling the random coefficients of the equation, using a spatial discretization of the partial differential equation, for example, by a Galerkin finite element method or finite difference method. Then the corresponding sample averages of desired statistics can be computed through these approximations. If the noise is described by a small number of random parameters or if the accuracy requirement is sufficiently strict, then the stochastic Galerkin method is preferred; otherwise, a MC method still seems to be the best choice.In this paper, we focus our study on the convection-diffusion equation with homogeneous Dirichlet boundary conditions. LetD be a convex bounded polygonal domain in R
d, d
=
1, 2, 3, and let (Ω
,
F
,
P) be a complete probability space. Here Ω is the set of outcomes, F
∈
2
Ω is the σ-algebra of events, and P
:
F
→
[
0,1
] is a probability measure. Consider the following stochastic linear convection-diffusion problem: find a stochastic function, u
:
D
¯
×
Ω
→
R, such that the following equation holds:
(1)
u
t
(
x
,
t
,
ω
)
+
c
(
x
)
∇
u
(
x
,
t
,
ω
)
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
t
,
ω
)
)
=
f
(
x
,
t
)
,
D
×
(
0
,
T
]
×
Ω
,
(2)
u
(
x
,
0
,
ω
)
=
u
0
(
x
,
ω
)
,
D
×
Ω
,
(3)
u
(
x
,
t
,
ω
)
=
0
,
∂
D
×
(
0
,
T
]
×
Ω
.
Here x
=
(
x
1, x
2, …, x
d
)
∈
D
⊂
R
d, d
=
1, 2, 3, are spatial coordinates, a: D
×
Ω
→
R is a stochastic function with continuous and bounded covariance function, and c is the velocity of flow. Our goal is to compute the quantity of interest (QoI):
(4)
Q
(
u
(
x
,
T
,
ω
)
)
=
∫
D
u
(
x
,
T
,
ω
)
d
x
.We consider the case when the random coefficient is in high random dimensions. We will employ finite volume element method (FVEM) [11] for discretization in physic space and MC method in random space. FVEM is used here because this method has its own advantages: reasonable accuracy, conservation of physical quantity locally, and high efficiency. The convergence rate of the MC is interpreted in the probability sense, and the estimate of its error needs a posterior estimate of variance of the sampled random variable, which in turn is a prior bound on higher statistical moments. The convergence rate of MC for the approximation of the expected value can be improved. In our paper, QMC offers a way to get a better convergence rate than MC.The paper is organized as follows. Our numerical method and corresponding theoretical aspects are described in Section2. In Section 3 we introduce the QMC and compare it with MC for numerical integral. In Section 4, we give the results of numerical examples. The conclusions are summarized in the last section.


## 2. Numerical Method and Theoretical Aspects

### 2.1. Karhunen-Loève Expansion and Finite Dimensional Noise

The Karhunen-Loève (KL) expansion [7] is one of the most widely used techniques for dimension reduction in representing a random process. It is based on the expansion of the correlation function of the process. Consider that the stochastic process a
(
x
,
ω
) in (1) is a stationary random field with continuous covariance function Cov
[
a
]
:
D
¯
×
D
¯
→
R. For well-posedness we require a
(
x
,
ω
)
≥
a
min
⁡
>
0. We expand a
(
x
,
ω
) using truncated KL expansion with N terms
(5)
a
(
x
,
ω
)
=
E
[
a
(
x
,
ω
)
]
+
∑
n
=
1
N
λ
n
z
n
(
ω
)
b
n
(
x
)
,
where E
[
a
(
x
,
ω
)
] is the mean of the input process a
(
x
,
ω
), {
z
n
} is a set of independent and identically distributed random variables, and λ
n, b
n denote the eigenvalues and eigenfunctions of the correlation function C with respect to the eigenvalue problem:
(6)
∫
D
C
(
x
,
y
)
e
n
(
y
)
d
y
=
λ
n
e
n
(
x
)
.
The summation in (5) is truncated at finite number N, and the number of terms N is determined by the decay of eigenvalues to ensure the approximation error is sufficiently small.The solution corresponding to the stochastic partial differential equation (1) can be described by just a finite number of random variables; that is, u
(
x
,
t
,
ω
)
=
u
(
x
,
t
,
z
1
(
ω
)
,
z
2
(
ω
)
,
…
,
z
N
(
ω
)
). We denote by Γ
n the bounded image set of the random variables z
n
(
Ω
), Γ
=
Γ
1
×
Γ
2
×
⋯
Γ
N, and we assume that the random vector z has a joint probability density function ρ
:
Γ
→
R
+ that factorizes as ρ
(
z
)
=
∏
n
=
1
N
‍
ρ
n
(
z
n
), for all z
∈
Γ
⊂
R
N. Then we can rewrite (1)–(3) as follows with an N-dimensional parameter:
(7)
u
t
(
x
,
t
,
z
)
+
c
(
x
)
∇
u
(
x
,
t
,
z
)
-
∇
·
(
a
(
x
,
z
)
∇
u
(
x
,
t
,
z
)
)
=
f
(
x
,
t
)
,
D
×
(
0
,
T
]
×
Γ
,
u
(
x
,
0
,
z
)
=
u
0
(
x
,
z
)
,
D
×
Γ
,
u
(
x
,
t
,
z
)
=
0
,
∂
D
×
(
0
,
T
]
×
Γ
.


### 2.2. Finite Volume Element Approximation in the Physical Space

The physical domainD can be decomposed into a grid T
h with a set of spaced nodes; accordingly we place dual grid T
h
* [11]. Select trial function space U
h
(
D
) as the finite element space of piecewise continuous polynomials on T
h, with dimension I
=
dim
⁡
(
U
h
(
D
)
) and basis {
ϕ
i
}
i
=
1
I. The test function space V
h
(
D
) corresponding to T
h
* is taken as the piecewise constant function space, with the same dimension as U
h
(
D
) and basis {
ψ
i
}
i
=
1
I. We introduce the space L
ρ
2
(
Γ
) of square integrable functions on Γ with respect to the measure ρ
(
z
)
d
z; that is, L
ρ
2
(
Γ
)
=
{
v
:
Γ
→
R
:
∥
v
∥
L
ρ
2
=
(
∫
Γ
‍
v
2
(
z
)
ρ
(
z
)
d
z
)
1
/
2
<
∞
}. Then we introduce the semidiscrete problem: find u
h
∈
U
h
(
D
)
⊗
L
ρ
2
(
Γ
) such that for all v
h
∈
V
h
(
D
)
(8)
∫
D
u
h
t
(
x
,
t
,
z
)
v
h
(
x
)
d
x
+
∫
D
[
c
(
x
)
∇
u
h
(
x
,
t
,
z
)
v
h
(
x
)
hhhhhh
+
a
(
x
,
z
)
∇
u
h
(
x
,
t
,
z
)
·
∇
v
h
(
x
)
]
d
x
=
∫
D
‍
f
(
x
,
t
)
v
h
(
x
)
d
x
.Problem (8) admits a unique solution for every z
∈
Γ. Then the semidiscrete solution can be expressed as u
h
(
x
,
t
,
z
)
=
∑
i
=
1
I
u
h
,
i
(
t
,
z
)
ϕ
i
(
x
). Denoting by U
(
t
,
z
)
=
[
u
h
,
1
(
t
,
z
)
,
u
h
,
2
(
t
,
z
)
,
…
,
u
h
,
I
(
t
,
z
)
] the vector of nodal values as functions of the random variables z and t, problem (8) can be written in algebraic form as
(9)
∂
U
(
t
,
z
)
∂
t
+
A
U
(
t
,
z
)
+
B
(
z
)
U
(
t
,
z
)
=
F
,
where A is deterministic matrix with respect to c
(
x
), matrix B
(
z
) is a stochastic stiffness matrix with respect to the random diffusion coefficient a
(
x
,
z
) which is an affine function of the random variables z
n, and F
i
=
∫
D
‍
f
(
x
,
t
)
ψ
i
(
x
)
d
x is a deterministic right hand side. The discretization in time may be performed by usual methods, for example, explicit or implicit Euler schemes and Runge-Kutta schemes.


### 2.3. Monte Carlo Finite Volume Element Method

In this section we describe the use of the Monte Carlo finite volume element method (MCFVEM) to construct approximations of the expected value of the solution and quantity of interest.Formulation of the MCFVEMStep 1.
Give a number of realizations,M, trial function space U
h
(
D
), and test function space V
h
(
D
), as defined in Section 2.2.Step 2.
For eachi
=
1,2
,
…
,
M, sample independent and identically distributed realizations of the diffusion coefficient a
(
ω
i
,
·
) based on realizations of {
z
n
}
n
=
1
N, and find a corresponding approximation u
h
(
ω
i
,
·
) which satisfies (9).Step 3.
Finally, use the sample average(
1
/
M
)
∑
i
=
1
M
‍
u
h
(
z
(
ω
i
)
) to approximate E
[
u
] and (
1
/
M
)
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
) to approximate E
[
Q
(
u
)
].


### 2.4.  MCFVEM Error Analysis

After the PDE is discretized by finite volume element method, we can compute discrete solutionsu
h
(
z
(
ω
i
)
), i
=
1,2
,
…
,
M. Our goal is to compute E
[
Q
(
u
(
z
)
)
], and then our estimator is
(10)
E
[
Q
(
u
(
z
)
)
]
≈
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
.The computational error naturally separates into two parts:(11)
|
E
[
Q
(
u
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
|
E
[
Q
(
u
(
z
)
)
]
-
E
[
Q
(
u
h
(
z
)
)
]
|
+
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
ε
h
+
ε
s
.
The size of the spatial discretization schemes controls the space discretization error ε
h, while the number of realizations, M of u
h, controls the statistical error ε
s.


#### 2.4.1. Discretization Error

Assumption 1.
The functionalQ
:
V
→
R, with Q
(
0
)
=
0, is globally Lipchitz; that is,
(12)
∃
C
Q
>
0
s
.
t
.
|
Q
(
u
)
-
Q
(
u
′
)
|
≤
C
Q
∥
u
-
u
′
∥
V
,
∀
u
,
u
′
∈
V
.Assumption 2.
There existα
>
0 and C
u
(
z
)
>
0, with ∫
Γ
‍
C
u
(
z
)
p
ρ
(
z
)
d
z
<
∞ for some p
>
1, such that
(13)
∥
u
(
z
)
-
u
h
(
z
)
∥
V
≤
C
u
(
z
)
h
α
,
∀
z
∈
Γ
,
0
<
h
<
h
0
.We apply Assumptions1 and 2 to estimate the discretization error
(14)
ε
h
=
|
E
[
Q
(
u
(
z
)
)
]
-
E
[
Q
(
u
h
(
z
)
)
]
|
=
|
∫
Γ
(
Q
(
u
(
z
)
)
-
Q
(
u
h
(
z
)
)
)
ρ
(
z
)
d
z
|
≤
C
Q
∫
Γ
∥
u
(
z
)
-
u
h
(
z
)
∥
V
ρ
(
z
)
d
z
≤
C
Q
h
α
∫
Γ
C
u
(
z
)
ρ
(
z
)
d
z
≤
C
Q
∥
C
u
∥
L
ρ
p
h
α
,
where ∥
C
u
∥
L
ρ
p
=
(
∫
‍
C
u
p
(
z
)
ρ
(
z
)
d
z
)
1
/
p and α depends on the physical space we choose. Once the trial function space and the test function space are certain, α is a fixed constant, the estimation for α; see [11].


#### 2.4.2. Statistical Error

Theorem 3 (the weak law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are independent, identically distributed (i.i.d) random variables and E
[
Y
i
]
=
μ
<
∞. Then
(15)
∑
j
=
1
M
Y
j
M
→
p
μ
,
a
s
M
⟶
∞
,
where →
p denotes convergence in probability; that is, the convergence (15) means P
(
|
∑
j
=
1
M
‍
(
Y
j
/
M
)
-
E
[
Y
]
|
>
ϵ
)
→
0 for all ϵ
>
0.Theorem 4 (the strong law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are i.i.d random variables with E
[
|
Y
i
|
]
<
∞ and E
[
Y
i
]
=
μ. Then
(16)
∑
j
=
1
M
Y
j
M
→
a
.
s
.
μ
,
a
s
M
⟶
∞
,
where →
a
.
s
. denotes almost sure convergence; that is, the convergence (16) means P
(
∑
j
=
1
M
‍
(
Y
j
/
M
)
→
μ
)
=
1.Theorem 5 (the central limit theorem).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(17)
∑
j
=
1
M
ξ
j
M
⇀
λ
,
where λ is N
(
0,1
) and ⇀ denotes convergence of the distributions, also called weak convergence; that is, the convergence (17) means that the following limit holds as M
→
∞:
(18)
E
[
g
(
1
M
∑
j
=
1
M
ξ
j
)
]
⟶
E
(
g
(
λ
)
)
for all bounded and continuous function g.Proof.
See [9, page 23].Theorem 6 (law of the iterated logarithm).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(19)
sup
⁡
⁡
|
∑
m
=
1
M
‍
ξ
m
|
2
M
log
⁡
⁡
log
⁡
⁡
M
→
a
.
s
.
1
,
a
s
M
⟶
∞
.Let us analyze the statistics error assuming thatQ is a globally Lipschitz function:
(20)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
|
1
M
∑
i
=
1
M
{
E
[
Q
(
u
h
(
z
)
)
]
-
Q
(
u
h
(
z
(
ω
i
)
)
)
}
|
,
E
[
ε
s
(
M
)
]
=
0. Then
(21)
Var
⁡
[
ε
s
(
M
)
]
=
1
M
Var
⁡
[
Q
(
u
h
(
z
)
)
]
and we can estimate
(22)
Var
⁡
[
Q
(
u
h
(
z
)
)
]
=
E
[
Q
(
u
h
(
z
)
)
2
]
≤
C
Q
2
∥
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
∥
u
-
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
)
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
h
2
α
∥
C
u
∥
L
ρ
2
(
Γ
)
2
)
<
∞
,
where ∥
u
∥
L
ρ
2
(
Γ
;
V
)
=
(
∫
Γ
∥
u
(
z
)
∥
V
2
ρ
(
z
)
d
z
)
1
/
2. We conclude that for Lipschitz functionals Q we have
(23)
Var
⁡
[
ε
s
(
M
)
]
≤
C
M
⟶
0
as
M
⟶
∞
with constant C uniformly bounded with respect to h. Since Var
⁡
[
Q
(
u
h
(
z
)
)
]
≤
C, we can apply the law of large numbers, law of iterated logarithms, and the central limit theorem. In particular, the previous result implies, via the Chebyshev inequality, convergence in probalility, that is, for any given ɛ
>
0
(24)
P
(
|
ε
s
(
M
)
|
>
ɛ
)
≤
Var
⁡
[
ε
s
(
M
)
]
ɛ
2
≤
C
M
ɛ
2
⟶
0
hhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
.LetY
=
Q
(
u
h
(
z
)
), Y
j
=
Q
(
u
h
(
z
(
ω
j
)
)
), Y
j, j
=
1
,
2
,
3
,
… be i.i.d random variables with E
[
Y
j
]
=
μ
<
∞ and Var
⁡
[
Y
j
]
=
σ
<
∞. We have
(25)
ε
s
(
M
)
=
|
1
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
.Letξ
j
=
(
Y
j
-
μ
)
/
σ, then ξ
j, j
=
1
,
2
,
… are i.i.d, and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then by Theorem 5, we have
(26)
M
ε
s
=
|
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
μ
M
|
=
|
∑
j
=
1
M
σ
ξ
j
+
μ
-
μ
M
|
=
∑
j
=
1
M
σ
ξ
j
M
⇀
σ
λ
,
where λ is N
(
0,1
).Motivated by central limit theorem (CLT), that is,(27)
P
(
|
M
ε
s
σ
|
≤
c
0
)
≈
P
(
|
λ
|
≤
c
0
)
⟶
2
Φ
(
c
0
)
-
1
,
hhhhhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
,
with Φ
(
x
)
=
∫
-
∞
x
‍
(
exp
⁡
⁡
(
-
x
2
/
2
)
/
2
π
)
d
x, we can write the statistical error as
(28)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
≤
c
0
std
[
Q
(
u
h
)
]
M
=
c
0
σ
M
.
The exact variance σ is in practice approximated by the sample variance:
(29)
σ
^
2
=
1
M
-
1
∑
j
=
1
M
(
Y
j
-
∑
m
=
1
M
Y
m
M
)
2
.Thus, we have the following CLT error bound:(30)
ε
s
≤
c
0
σ
^
M
,
where c
0 is an appropriate constant.


## 2.1. Karhunen-Loève Expansion and Finite Dimensional Noise

The Karhunen-Loève (KL) expansion [7] is one of the most widely used techniques for dimension reduction in representing a random process. It is based on the expansion of the correlation function of the process. Consider that the stochastic process a
(
x
,
ω
) in (1) is a stationary random field with continuous covariance function Cov
[
a
]
:
D
¯
×
D
¯
→
R. For well-posedness we require a
(
x
,
ω
)
≥
a
min
⁡
>
0. We expand a
(
x
,
ω
) using truncated KL expansion with N terms
(5)
a
(
x
,
ω
)
=
E
[
a
(
x
,
ω
)
]
+
∑
n
=
1
N
λ
n
z
n
(
ω
)
b
n
(
x
)
,
where E
[
a
(
x
,
ω
)
] is the mean of the input process a
(
x
,
ω
), {
z
n
} is a set of independent and identically distributed random variables, and λ
n, b
n denote the eigenvalues and eigenfunctions of the correlation function C with respect to the eigenvalue problem:
(6)
∫
D
C
(
x
,
y
)
e
n
(
y
)
d
y
=
λ
n
e
n
(
x
)
.
The summation in (5) is truncated at finite number N, and the number of terms N is determined by the decay of eigenvalues to ensure the approximation error is sufficiently small.The solution corresponding to the stochastic partial differential equation (1) can be described by just a finite number of random variables; that is, u
(
x
,
t
,
ω
)
=
u
(
x
,
t
,
z
1
(
ω
)
,
z
2
(
ω
)
,
…
,
z
N
(
ω
)
). We denote by Γ
n the bounded image set of the random variables z
n
(
Ω
), Γ
=
Γ
1
×
Γ
2
×
⋯
Γ
N, and we assume that the random vector z has a joint probability density function ρ
:
Γ
→
R
+ that factorizes as ρ
(
z
)
=
∏
n
=
1
N
‍
ρ
n
(
z
n
), for all z
∈
Γ
⊂
R
N. Then we can rewrite (1)–(3) as follows with an N-dimensional parameter:
(7)
u
t
(
x
,
t
,
z
)
+
c
(
x
)
∇
u
(
x
,
t
,
z
)
-
∇
·
(
a
(
x
,
z
)
∇
u
(
x
,
t
,
z
)
)
=
f
(
x
,
t
)
,
D
×
(
0
,
T
]
×
Γ
,
u
(
x
,
0
,
z
)
=
u
0
(
x
,
z
)
,
D
×
Γ
,
u
(
x
,
t
,
z
)
=
0
,
∂
D
×
(
0
,
T
]
×
Γ
.


## 2.2. Finite Volume Element Approximation in the Physical Space

The physical domainD can be decomposed into a grid T
h with a set of spaced nodes; accordingly we place dual grid T
h
* [11]. Select trial function space U
h
(
D
) as the finite element space of piecewise continuous polynomials on T
h, with dimension I
=
dim
⁡
(
U
h
(
D
)
) and basis {
ϕ
i
}
i
=
1
I. The test function space V
h
(
D
) corresponding to T
h
* is taken as the piecewise constant function space, with the same dimension as U
h
(
D
) and basis {
ψ
i
}
i
=
1
I. We introduce the space L
ρ
2
(
Γ
) of square integrable functions on Γ with respect to the measure ρ
(
z
)
d
z; that is, L
ρ
2
(
Γ
)
=
{
v
:
Γ
→
R
:
∥
v
∥
L
ρ
2
=
(
∫
Γ
‍
v
2
(
z
)
ρ
(
z
)
d
z
)
1
/
2
<
∞
}. Then we introduce the semidiscrete problem: find u
h
∈
U
h
(
D
)
⊗
L
ρ
2
(
Γ
) such that for all v
h
∈
V
h
(
D
)
(8)
∫
D
u
h
t
(
x
,
t
,
z
)
v
h
(
x
)
d
x
+
∫
D
[
c
(
x
)
∇
u
h
(
x
,
t
,
z
)
v
h
(
x
)
hhhhhh
+
a
(
x
,
z
)
∇
u
h
(
x
,
t
,
z
)
·
∇
v
h
(
x
)
]
d
x
=
∫
D
‍
f
(
x
,
t
)
v
h
(
x
)
d
x
.Problem (8) admits a unique solution for every z
∈
Γ. Then the semidiscrete solution can be expressed as u
h
(
x
,
t
,
z
)
=
∑
i
=
1
I
u
h
,
i
(
t
,
z
)
ϕ
i
(
x
). Denoting by U
(
t
,
z
)
=
[
u
h
,
1
(
t
,
z
)
,
u
h
,
2
(
t
,
z
)
,
…
,
u
h
,
I
(
t
,
z
)
] the vector of nodal values as functions of the random variables z and t, problem (8) can be written in algebraic form as
(9)
∂
U
(
t
,
z
)
∂
t
+
A
U
(
t
,
z
)
+
B
(
z
)
U
(
t
,
z
)
=
F
,
where A is deterministic matrix with respect to c
(
x
), matrix B
(
z
) is a stochastic stiffness matrix with respect to the random diffusion coefficient a
(
x
,
z
) which is an affine function of the random variables z
n, and F
i
=
∫
D
‍
f
(
x
,
t
)
ψ
i
(
x
)
d
x is a deterministic right hand side. The discretization in time may be performed by usual methods, for example, explicit or implicit Euler schemes and Runge-Kutta schemes.


## 2.3. Monte Carlo Finite Volume Element Method

In this section we describe the use of the Monte Carlo finite volume element method (MCFVEM) to construct approximations of the expected value of the solution and quantity of interest.Formulation of the MCFVEMStep 1.
Give a number of realizations,M, trial function space U
h
(
D
), and test function space V
h
(
D
), as defined in Section 2.2.Step 2.
For eachi
=
1,2
,
…
,
M, sample independent and identically distributed realizations of the diffusion coefficient a
(
ω
i
,
·
) based on realizations of {
z
n
}
n
=
1
N, and find a corresponding approximation u
h
(
ω
i
,
·
) which satisfies (9).Step 3.
Finally, use the sample average(
1
/
M
)
∑
i
=
1
M
‍
u
h
(
z
(
ω
i
)
) to approximate E
[
u
] and (
1
/
M
)
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
) to approximate E
[
Q
(
u
)
].


## 2.4.  MCFVEM Error Analysis

After the PDE is discretized by finite volume element method, we can compute discrete solutionsu
h
(
z
(
ω
i
)
), i
=
1,2
,
…
,
M. Our goal is to compute E
[
Q
(
u
(
z
)
)
], and then our estimator is
(10)
E
[
Q
(
u
(
z
)
)
]
≈
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
.The computational error naturally separates into two parts:(11)
|
E
[
Q
(
u
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
|
E
[
Q
(
u
(
z
)
)
]
-
E
[
Q
(
u
h
(
z
)
)
]
|
+
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
ε
h
+
ε
s
.
The size of the spatial discretization schemes controls the space discretization error ε
h, while the number of realizations, M of u
h, controls the statistical error ε
s.


### 2.4.1. Discretization Error

Assumption 1.
The functionalQ
:
V
→
R, with Q
(
0
)
=
0, is globally Lipchitz; that is,
(12)
∃
C
Q
>
0
s
.
t
.
|
Q
(
u
)
-
Q
(
u
′
)
|
≤
C
Q
∥
u
-
u
′
∥
V
,
∀
u
,
u
′
∈
V
.Assumption 2.
There existα
>
0 and C
u
(
z
)
>
0, with ∫
Γ
‍
C
u
(
z
)
p
ρ
(
z
)
d
z
<
∞ for some p
>
1, such that
(13)
∥
u
(
z
)
-
u
h
(
z
)
∥
V
≤
C
u
(
z
)
h
α
,
∀
z
∈
Γ
,
0
<
h
<
h
0
.We apply Assumptions1 and 2 to estimate the discretization error
(14)
ε
h
=
|
E
[
Q
(
u
(
z
)
)
]
-
E
[
Q
(
u
h
(
z
)
)
]
|
=
|
∫
Γ
(
Q
(
u
(
z
)
)
-
Q
(
u
h
(
z
)
)
)
ρ
(
z
)
d
z
|
≤
C
Q
∫
Γ
∥
u
(
z
)
-
u
h
(
z
)
∥
V
ρ
(
z
)
d
z
≤
C
Q
h
α
∫
Γ
C
u
(
z
)
ρ
(
z
)
d
z
≤
C
Q
∥
C
u
∥
L
ρ
p
h
α
,
where ∥
C
u
∥
L
ρ
p
=
(
∫
‍
C
u
p
(
z
)
ρ
(
z
)
d
z
)
1
/
p and α depends on the physical space we choose. Once the trial function space and the test function space are certain, α is a fixed constant, the estimation for α; see [11].


### 2.4.2. Statistical Error

Theorem 3 (the weak law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are independent, identically distributed (i.i.d) random variables and E
[
Y
i
]
=
μ
<
∞. Then
(15)
∑
j
=
1
M
Y
j
M
→
p
μ
,
a
s
M
⟶
∞
,
where →
p denotes convergence in probability; that is, the convergence (15) means P
(
|
∑
j
=
1
M
‍
(
Y
j
/
M
)
-
E
[
Y
]
|
>
ϵ
)
→
0 for all ϵ
>
0.Theorem 4 (the strong law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are i.i.d random variables with E
[
|
Y
i
|
]
<
∞ and E
[
Y
i
]
=
μ. Then
(16)
∑
j
=
1
M
Y
j
M
→
a
.
s
.
μ
,
a
s
M
⟶
∞
,
where →
a
.
s
. denotes almost sure convergence; that is, the convergence (16) means P
(
∑
j
=
1
M
‍
(
Y
j
/
M
)
→
μ
)
=
1.Theorem 5 (the central limit theorem).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(17)
∑
j
=
1
M
ξ
j
M
⇀
λ
,
where λ is N
(
0,1
) and ⇀ denotes convergence of the distributions, also called weak convergence; that is, the convergence (17) means that the following limit holds as M
→
∞:
(18)
E
[
g
(
1
M
∑
j
=
1
M
ξ
j
)
]
⟶
E
(
g
(
λ
)
)
for all bounded and continuous function g.Proof.
See [9, page 23].Theorem 6 (law of the iterated logarithm).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(19)
sup
⁡
⁡
|
∑
m
=
1
M
‍
ξ
m
|
2
M
log
⁡
⁡
log
⁡
⁡
M
→
a
.
s
.
1
,
a
s
M
⟶
∞
.Let us analyze the statistics error assuming thatQ is a globally Lipschitz function:
(20)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
|
1
M
∑
i
=
1
M
{
E
[
Q
(
u
h
(
z
)
)
]
-
Q
(
u
h
(
z
(
ω
i
)
)
)
}
|
,
E
[
ε
s
(
M
)
]
=
0. Then
(21)
Var
⁡
[
ε
s
(
M
)
]
=
1
M
Var
⁡
[
Q
(
u
h
(
z
)
)
]
and we can estimate
(22)
Var
⁡
[
Q
(
u
h
(
z
)
)
]
=
E
[
Q
(
u
h
(
z
)
)
2
]
≤
C
Q
2
∥
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
∥
u
-
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
)
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
h
2
α
∥
C
u
∥
L
ρ
2
(
Γ
)
2
)
<
∞
,
where ∥
u
∥
L
ρ
2
(
Γ
;
V
)
=
(
∫
Γ
∥
u
(
z
)
∥
V
2
ρ
(
z
)
d
z
)
1
/
2. We conclude that for Lipschitz functionals Q we have
(23)
Var
⁡
[
ε
s
(
M
)
]
≤
C
M
⟶
0
as
M
⟶
∞
with constant C uniformly bounded with respect to h. Since Var
⁡
[
Q
(
u
h
(
z
)
)
]
≤
C, we can apply the law of large numbers, law of iterated logarithms, and the central limit theorem. In particular, the previous result implies, via the Chebyshev inequality, convergence in probalility, that is, for any given ɛ
>
0
(24)
P
(
|
ε
s
(
M
)
|
>
ɛ
)
≤
Var
⁡
[
ε
s
(
M
)
]
ɛ
2
≤
C
M
ɛ
2
⟶
0
hhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
.LetY
=
Q
(
u
h
(
z
)
), Y
j
=
Q
(
u
h
(
z
(
ω
j
)
)
), Y
j, j
=
1
,
2
,
3
,
… be i.i.d random variables with E
[
Y
j
]
=
μ
<
∞ and Var
⁡
[
Y
j
]
=
σ
<
∞. We have
(25)
ε
s
(
M
)
=
|
1
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
.Letξ
j
=
(
Y
j
-
μ
)
/
σ, then ξ
j, j
=
1
,
2
,
… are i.i.d, and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then by Theorem 5, we have
(26)
M
ε
s
=
|
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
μ
M
|
=
|
∑
j
=
1
M
σ
ξ
j
+
μ
-
μ
M
|
=
∑
j
=
1
M
σ
ξ
j
M
⇀
σ
λ
,
where λ is N
(
0,1
).Motivated by central limit theorem (CLT), that is,(27)
P
(
|
M
ε
s
σ
|
≤
c
0
)
≈
P
(
|
λ
|
≤
c
0
)
⟶
2
Φ
(
c
0
)
-
1
,
hhhhhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
,
with Φ
(
x
)
=
∫
-
∞
x
‍
(
exp
⁡
⁡
(
-
x
2
/
2
)
/
2
π
)
d
x, we can write the statistical error as
(28)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
≤
c
0
std
[
Q
(
u
h
)
]
M
=
c
0
σ
M
.
The exact variance σ is in practice approximated by the sample variance:
(29)
σ
^
2
=
1
M
-
1
∑
j
=
1
M
(
Y
j
-
∑
m
=
1
M
Y
m
M
)
2
.Thus, we have the following CLT error bound:(30)
ε
s
≤
c
0
σ
^
M
,
where c
0 is an appropriate constant.


## 2.4.1. Discretization Error

Assumption 1.
The functionalQ
:
V
→
R, with Q
(
0
)
=
0, is globally Lipchitz; that is,
(12)
∃
C
Q
>
0
s
.
t
.
|
Q
(
u
)
-
Q
(
u
′
)
|
≤
C
Q
∥
u
-
u
′
∥
V
,
∀
u
,
u
′
∈
V
.Assumption 2.
There existα
>
0 and C
u
(
z
)
>
0, with ∫
Γ
‍
C
u
(
z
)
p
ρ
(
z
)
d
z
<
∞ for some p
>
1, such that
(13)
∥
u
(
z
)
-
u
h
(
z
)
∥
V
≤
C
u
(
z
)
h
α
,
∀
z
∈
Γ
,
0
<
h
<
h
0
.We apply Assumptions1 and 2 to estimate the discretization error
(14)
ε
h
=
|
E
[
Q
(
u
(
z
)
)
]
-
E
[
Q
(
u
h
(
z
)
)
]
|
=
|
∫
Γ
(
Q
(
u
(
z
)
)
-
Q
(
u
h
(
z
)
)
)
ρ
(
z
)
d
z
|
≤
C
Q
∫
Γ
∥
u
(
z
)
-
u
h
(
z
)
∥
V
ρ
(
z
)
d
z
≤
C
Q
h
α
∫
Γ
C
u
(
z
)
ρ
(
z
)
d
z
≤
C
Q
∥
C
u
∥
L
ρ
p
h
α
,
where ∥
C
u
∥
L
ρ
p
=
(
∫
‍
C
u
p
(
z
)
ρ
(
z
)
d
z
)
1
/
p and α depends on the physical space we choose. Once the trial function space and the test function space are certain, α is a fixed constant, the estimation for α; see [11].


## 2.4.2. Statistical Error

Theorem 3 (the weak law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are independent, identically distributed (i.i.d) random variables and E
[
Y
i
]
=
μ
<
∞. Then
(15)
∑
j
=
1
M
Y
j
M
→
p
μ
,
a
s
M
⟶
∞
,
where →
p denotes convergence in probability; that is, the convergence (15) means P
(
|
∑
j
=
1
M
‍
(
Y
j
/
M
)
-
E
[
Y
]
|
>
ϵ
)
→
0 for all ϵ
>
0.Theorem 4 (the strong law of large numbers).
AssumeY
j, j
=
1
,
2
,
3
,
… are i.i.d random variables with E
[
|
Y
i
|
]
<
∞ and E
[
Y
i
]
=
μ. Then
(16)
∑
j
=
1
M
Y
j
M
→
a
.
s
.
μ
,
a
s
M
⟶
∞
,
where →
a
.
s
. denotes almost sure convergence; that is, the convergence (16) means P
(
∑
j
=
1
M
‍
(
Y
j
/
M
)
→
μ
)
=
1.Theorem 5 (the central limit theorem).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(17)
∑
j
=
1
M
ξ
j
M
⇀
λ
,
where λ is N
(
0,1
) and ⇀ denotes convergence of the distributions, also called weak convergence; that is, the convergence (17) means that the following limit holds as M
→
∞:
(18)
E
[
g
(
1
M
∑
j
=
1
M
ξ
j
)
]
⟶
E
(
g
(
λ
)
)
for all bounded and continuous function g.Proof.
See [9, page 23].Theorem 6 (law of the iterated logarithm).
Assumeξ
j, j
=
1
,
2
,
… are i.i.d and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then
(19)
sup
⁡
⁡
|
∑
m
=
1
M
‍
ξ
m
|
2
M
log
⁡
⁡
log
⁡
⁡
M
→
a
.
s
.
1
,
a
s
M
⟶
∞
.Let us analyze the statistics error assuming thatQ is a globally Lipschitz function:
(20)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
=
|
1
M
∑
i
=
1
M
{
E
[
Q
(
u
h
(
z
)
)
]
-
Q
(
u
h
(
z
(
ω
i
)
)
)
}
|
,
E
[
ε
s
(
M
)
]
=
0. Then
(21)
Var
⁡
[
ε
s
(
M
)
]
=
1
M
Var
⁡
[
Q
(
u
h
(
z
)
)
]
and we can estimate
(22)
Var
⁡
[
Q
(
u
h
(
z
)
)
]
=
E
[
Q
(
u
h
(
z
)
)
2
]
≤
C
Q
2
∥
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
∥
u
-
u
h
∥
L
ρ
2
(
Γ
;
V
)
2
)
≤
2
C
Q
2
(
∥
u
∥
L
ρ
2
(
Γ
;
V
)
2
+
h
2
α
∥
C
u
∥
L
ρ
2
(
Γ
)
2
)
<
∞
,
where ∥
u
∥
L
ρ
2
(
Γ
;
V
)
=
(
∫
Γ
∥
u
(
z
)
∥
V
2
ρ
(
z
)
d
z
)
1
/
2. We conclude that for Lipschitz functionals Q we have
(23)
Var
⁡
[
ε
s
(
M
)
]
≤
C
M
⟶
0
as
M
⟶
∞
with constant C uniformly bounded with respect to h. Since Var
⁡
[
Q
(
u
h
(
z
)
)
]
≤
C, we can apply the law of large numbers, law of iterated logarithms, and the central limit theorem. In particular, the previous result implies, via the Chebyshev inequality, convergence in probalility, that is, for any given ɛ
>
0
(24)
P
(
|
ε
s
(
M
)
|
>
ɛ
)
≤
Var
⁡
[
ε
s
(
M
)
]
ɛ
2
≤
C
M
ɛ
2
⟶
0
hhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
.LetY
=
Q
(
u
h
(
z
)
), Y
j
=
Q
(
u
h
(
z
(
ω
j
)
)
), Y
j, j
=
1
,
2
,
3
,
… be i.i.d random variables with E
[
Y
j
]
=
μ
<
∞ and Var
⁡
[
Y
j
]
=
σ
<
∞. We have
(25)
ε
s
(
M
)
=
|
1
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
.Letξ
j
=
(
Y
j
-
μ
)
/
σ, then ξ
j, j
=
1
,
2
,
… are i.i.d, and E
[
ξ
j
]
=
0, E
[
ξ
j
2
]
=
1. Then by Theorem 5, we have
(26)
M
ε
s
=
|
M
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
E
[
Y
]
M
|
=
|
∑
j
=
1
M
Y
j
-
μ
M
|
=
|
∑
j
=
1
M
σ
ξ
j
+
μ
-
μ
M
|
=
∑
j
=
1
M
σ
ξ
j
M
⇀
σ
λ
,
where λ is N
(
0,1
).Motivated by central limit theorem (CLT), that is,(27)
P
(
|
M
ε
s
σ
|
≤
c
0
)
≈
P
(
|
λ
|
≤
c
0
)
⟶
2
Φ
(
c
0
)
-
1
,
hhhhhhhhhhhhhhhhhhhhhhhhhhh
as
M
⟶
∞
,
with Φ
(
x
)
=
∫
-
∞
x
‍
(
exp
⁡
⁡
(
-
x
2
/
2
)
/
2
π
)
d
x, we can write the statistical error as
(28)
ε
s
(
M
)
=
|
E
[
Q
(
u
h
(
z
)
)
]
-
1
M
∑
i
=
1
M
Q
(
u
h
(
z
(
ω
i
)
)
)
|
≤
c
0
std
[
Q
(
u
h
)
]
M
=
c
0
σ
M
.
The exact variance σ is in practice approximated by the sample variance:
(29)
σ
^
2
=
1
M
-
1
∑
j
=
1
M
(
Y
j
-
∑
m
=
1
M
Y
m
M
)
2
.Thus, we have the following CLT error bound:(30)
ε
s
≤
c
0
σ
^
M
,
where c
0 is an appropriate constant.


## 3. Monte Carlo and Quasi-Monte Carlo

The convergence of MC for numerical integration can often be improved by replacing random sequences with more uniformly distributed sequences known as quasirandom sequences. Quasirandom sequences are the low-discrepancy sequences and have greater uniformity than random sequences. They combine the advantage of a random sequence with the advantage of a lattice. Examples of such sequences include the Halton sequence, Sobol sequence, and Faure sequence [12, 13]. We will use QMC [10, 14, 15] method with Sobol sequence in the following numerical simulation for better performance.Figure1 shows the two-dimensional scatter for MC points and QMC points over the unit square. The points distributed on the left are intense at some place and sparse at other places, while the points on the right are more uniformly distributed.Two-dimensional scatter plot for MC points (a) and QMC points (b).
(a)
(b)We calculate the quantityI
K
(
f
)
=
∫
[
0,1
]
K
‍
f
(
x
)
d
x using an equal weight cubature formula of the form
(31)
I
K
,
M
(
f
)
=
1
M
∑
i
=
1
M
f
(
x
i
)
,
where x is a sequence of K random points in the K dimensional cube [
0,1
]
K, and f
:
[
0,1
]
K
→
R is a continuous function.To estimate the error in a QMC integration, a shifted averaged lattice rule [16–20] will be employed. Let η
~
U
(
[
0,1
]
K
) be a uniformly distributed random vector (shift) and
(32)
I
K
,
M
(
f
,
η
)
=
1
M
∑
i
=
1
M
f
(
{
x
i
+
η
}
)
the QMC formula with the random shift η. Observe that the shift preserves the lattice structure and E
[
f
(
{
x
i
+
η
}
)
]
=
I
K
(
f
), for all i
=
1
,
2
,
…
,
M. Take s i.i.d shifts η
1
,
…
,
η
s
~
U
(
[
0,1
]
K
). Then the estimate of the integral is
(33)
I
K
,
M
(
f
)
=
1
s
∑
j
=
1
s
I
K
,
M
(
f
,
η
j
)
=
1
s
M
∑
j
=
1
s
∑
i
=
1
M
f
(
{
x
i
+
η
j
}
)
.
Estimate of QMC error based on the η-sample standard deviation, that is,
(34)
e
K
,
M
(
f
)
=
c
0
s
1
(
s
-
1
)
∑
j
=
1
s
(
I
K
,
M
(
f
,
η
j
)
-
I
K
,
M
(
f
)
)
2
.
In the following, we will use MC and QMC to numerically evaluate the integrals.Model 1.
Consider the following:f
(
x
)
=
cos
⁡
⁡
(
2
π
w
1
+
∑
n
=
1
K
‍
c
n
x
n
), with c
n
=
9
/
K, w
1
=
1
/
2. The exact solution is
(35)
∫
[
0,1
]
k
f
(
x
)
d
x
=
R
(
e
i
2
π
w
1
∏
n
=
1
K
1
i
c
n
(
e
i
c
n
-
1
)
)
,
where i denotes the imaginary unit and R
(
z
) is the real part of z
∈
C.Model 2.
Consider the following:(36)
f
(
x
)
=
{
0
,
if
x
1
>
w
1
or
x
2
>
w
2
,
exp
⁡
⁡
(
∑
n
=
1
K
c
n
x
n
)
,
otherwise
,
with c
i
=
4.3
/
K, w
1
=
π
/
4, and w
2
=
π
/
5. The exact solution is
(37)
∫
[
0,1
]
k
f
(
x
)
d
x
=
1
∏
n
=
1
K
‍
c
n
(
e
c
1
w
1
-
1
)
(
e
c
2
w
2
-
1
)
∏
n
=
3
K
(
e
c
n
-
1
)
.
To investigate the error convergence rate of MC and QMC on evaluating the integral, we will compare the CLT error bound and the exact error. The CLT error of MC is defined in (30) and the CLT error of QMC is defined in (34). For both MC and QMC, the exact error is defined as the absolute error between the exact solution and the numerical solution.
It seems plausible that if the random sequence is replaced by Sobol sequence, smaller error may result. The errors are far smaller in QMC than MC, which can be confirmed by Figures2 and 3. As the convergence rate, the CLT error of MC is near to −0.5, which corresponds to the theoretical error bound. For QMC, the convergence rate is larger; in other words, QMC error converges faster than MC, which is a very important part of the motivation of QMC. Therefore, QMC is a valid technique to accelerate convergence and improve efficiency, and it will be used in the following stochastic partial differential equations.Figure 2
The exact error and CLT error estimation for the integral of Model1, with K
=
10.Figure 3
The exact error and CLT error estimation for the integral of Model2, with K
=
10.


## 4. Numerical Experiments

In this section, we will numerically study examples of stochastic partial differential equation with random coefficient, using the methods described above.To keep the stochastic processa
(
x
,
ω
) in (1) strictly positive, we assume a
(
x
,
ω
)
=
e
k
(
x
,
ω
), where k
(
x
,
ω
) is a stationary random field with the Matérn covariance function:
(38)
C
(
x
,
y
)
=
γ
1
Γ
(
ν
)
2
ν
-
1
(
2
ν
∥
x
-
y
∥
ρ
)
ν
K
ν
(
2
ν
∥
x
-
y
∥
ρ
)
,
where Γ is the gamma function and K
ν is the modified Bessel function of the second kind. The following special cases of C are considered:
(39)
ν
=
0.5
,
C
(
x
,
y
)
=
γ
e
(
-
∥
x
-
y
∥
/
ρ
)
,
ν
=
1.5
,
C
(
x
,
y
)
=
γ
(
1
+
3
d
ρ
)
e
(
-
3
∥
x
-
y
∥
/
ρ
)
,
ν
=
2.5
,
C
(
x
,
y
)

=

γ
(
1

+
5
d
ρ
+

3
d
2
ρ
2
)
e
(
-
5
∥
x
-
y
∥
/
ρ
)

,

ν
⟶
∞
,
C
(
x
,
y
)
=
γ
e
(
-
∥
x
-
y
∥
/
2
ρ
2
)
,
where ρ scales the correlation length, d denotes the distance between two spatial points, and γ is a constant. In the following, we choose γ
=
2 and ρ
=
0.1.We expandk
(
x
,
ω
) using truncated KL expansion with N terms:
(40)
k
(
x
,
ω
)
=
∑
n
=
1
N
λ
n
z
n
(
ω
)
e
n
(
x
)
,
where {
z
k
} is a set of i.i.d uniform random variables over [
-
3
,
3
] with mean value zero and unit variance. To find the eigenvalues and eigenfunctions of C we solve the eigenvalue problem:
(41)
∫
D
C
(
x
,
y
)
e
n
(
y
)
d
y
=
λ
n
e
n
(
x
)
.
As no analytical solutions are available for the eigenvalue problem (41), a Galerkin numerical eigenvalue solver is employed [21]. We solve (41) by discretizing C as a matrix, and the discretization corresponds to a piecewise constant FEM approximation to (41), and then MATLAB’s function eig
(
) is used. For other numerical methods, see [22–24].The spectral contribution ofC to the KL expansion is plotted in Figure 4, which shows the decay of the eigenvalues with different choice of parameter ν. From it we can see that when ν increases from 0.5 to ∞, the spectral convergence becomes growingly faster. In other words, the decay of eigenvalues is faster with larger choice of ν. Therefore, for large ν, the number of truncated terms in KL expansion can be chosen smaller, which means the dimension of the problem can be reduced. In the following problems, the number of truncated terms in (40) is N
=
20.Figure 4
This figure shows spectral contribution ofC to the KL expansion; that is, λ
n
∥
e
n
∥
∞, and x
∈
(
0,1
).


### 4.1. Steady State Diffusion Problem

We first consider the following problem in one spatial dimension (d
=
1) and N
(
=
20
) random dimensions:
(42)
-
∂
∂
x
[
a
(
x
,
z
)
∂
u
∂
x
(
x
,
z
)
]
=
cos
⁡
⁡
(
2
π
x
)
,
x
∈
(
0,1
)
,
z
∈
R
N
,
(43)
u
(
0
,
·
)
=
u
(
1
,
·
)
=
0
.The region of interest[
0,1
] can be decomposed into a uniform grid T
h with nodes 0
=
x
0
<
x
1
<
⋯
<
x
I
=
1 and h
=
x
i
-
x
i
-
1
=
1
/
I, i
=
1
,
2
,
…
,
I. To obtain quadratic basis functions, we will use the midpoints x
i
-
1
/
2
=
(
x
i
+
x
i
-
1
)
/
2 as the interpolation points. Accordingly we place a dual grid T
h
* with nodes 0
=
x
0
<
x
1
/
4
<
x
3
/
4
<
⋯
<
x
I
-
3
/
4
<
x
I
-
1
/
4
<
x
I
=
1, where x
i
-
k
/
4
=
x
i
-
(
k
/
4
)
h, (k
=
1
,
3; i
=
1
,
2
,
…
,
I). Select trial function space U
h as the quadratic element space of Lagrange type with respect to T
h, and the test function space V
h corresponding to T
h
* is taken as the piecewise constant function space. Then the numerical solution u
h for (42) can be uniquely written as
(44)
u
h
(
x
,
z
)
=
∑
i
=
1
I
[
u
h
,
i
(
z
)
ϕ
i
(
x
)
+
u
h
,
i
-
1
/
2
(
z
)
ϕ
i
-
1
/
2
(
x
)
]
,
which is a piecewise quadratic FVEM approximation. We obtain the discrete schemes:
(45)
-
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
-
1
/
2
(
z
)
+
2
a
(
x
j
-
1
/
4
,
z
)
+
2
a
(
x
j
+
1
/
4
,
z
)
h
u
h
,
j
(
z
)
-
2
a
(
x
j
+
1
/
4
,
z
)
h
u
h
,
j
+
1
/
2
(
z
)
=
∫
x
j
-
1
/
4
x
j
+
1
/
4
cos
⁡
⁡
(
2
π
x
)
d
x
,
j
=
1,2
,
…
,
I
-
1
,
-
2
a
(
x
j
-
3
/
4
,
z
)
h
u
h
,
j
-
1
(
z
)
+
2
a
(
x
j
-
3
/
4
,
z
)
+
2
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
-
1
/
2
(
z
)
-
2
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
(
z
)
=
∫
x
j
-
3
/
4
x
j
-
1
/
4
‍
cos
⁡
⁡
(
2
π
x
)
d
x
,
j
=
1,2
,
…
,
I
,
and they yield a tridiagonal linear system for the nodal values:
(46)
A
(
z
)
U
h
(
z
)
=
F
.The value of QoI:Q
(
u
(
z
)
)
=
∫
0
1
‍
u
(
x
,
z
)
d
x can be computed exactly by a Simpson formula, yielding
(47)
Q
(
u
(
z
)
)
=
∑
i
=
1
I
h
3
(
u
h
,
i
-
1
/
2
+
u
h
,
i
+
u
h
,
i
+
1
/
2
)
.The numerical approximations of QoI using both MC method and QMC method are presented in Figures5 and 6 with different choice of ν. We can see from both of them that when the number of samples becomes larger and larger, the approximations gradually tend to a fixed value. Besides, approximations using QMC method in Figure 6 are more smooth than that in Figure 5, which means QMC method is more accurate than MC method regardless of the value of ν.Figure 5
MC approximation of QoI with different choice ofν.Figure 6
QMC approximation of QoI with different choice ofν.Figures7 and 8 show error estimates of MC method and QMC method when ν is chosen as 1.5 and 2.5, respectively. It is obvious that the QMC error is smaller than the MC error. Their slopes are also estimated as a comparison, and we can see that the slopes of QMC are bigger than that of MC. That means, with the same number of samples, estimate of QMC with Sobol sequence converges faster than that of MC method.Figure 7
Error estimate of MC (dotted line) and QMC (solid line) withν
=
1.5.Figure 8
Error estimate of MC (dotted line) and QMC (solid line) withν
=
2.5.


### 4.2. Convection-Diffusion Problem

Finally we apply the MCFVEM algorithm to the stochastic convection-diffusion equation:(48)
∂
u
∂
t
+
c
(
x
)
∂
u
∂
x
-
∂
∂
x
[
a
(
x
,
z
)
∂
u
∂
x
]
=
0
,
hhhhh
x
∈
(
0,1
)
,
z
∈
R
N
,
t
∈
[
0
,
T
]
,
(49)
u
(
0
,
t
,
·
)
=
1
,
u
(
1
,
t
,
·
)
=
0
,
u
(
x
,
0
,
·
)
=
u
0
(
x
,
·
)
.
The exact solution of (48) is
(50)
u
(
x
,
t
,
z
)
=
1
2
{

[
x
+
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
)
erfc
⁡
[
x
-
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
]
hh
+
exp
⁡
⁡
(
c
x
a
(
x
,
z
)
)
erfc
⁡
[
x
+
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
]
}
.
In this case we are interested in computing the expected value of the following QoI:
(51)
Q
(
u
(
x
,
T
,
ω
)
)
=
∫
0
1
u
(
x
,
T
,
ω
)
d
x
.
It can be computed by a Simpson formula similarly as in (47).Use the same physical space as in Section4.1, and employ implicit Euler scheme to discretize the first-order time derivative. We obtain the following fully discrete quadratic finite volume element schemes:
(52)
1
Δ
t
(
u
h
,
j
n
+
1
(
z
)
-
u
h
,
j
n
(
z
)
)
+
2
h
c
j
n
+
1
(
1
8
u
h
,
j
-
1
n
+
1
(
z
)
-
3
4
u
h
,
j
-
1
/
2
n
+
1
(
z
)
hhhhhhhhh
+
3
4
u
h
,
j
+
1
/
2
n
+
1
(
z
)
-
1
8
u
h
,
j
+
1
n
+
1
(
z
)
)
-
2
h
2
(
2
a
j
+
1
/
4
n
+
1
(
z
)
u
h
,
j
+
1
/
2
n
+
1
(
z
)
hhhhhh
-
2
(
a
j
+
1
/
4
n
+
1
(
z
)
+
a
j
-
1
/
4
n
+
1
(
z
)
)
u
h
,
j
n
+
1
(
z
)
hhhhhh
+
2
a
j
-
1
/
4
n
+
1
(
z
)
u
h
,
j
-
1
/
2
n
+
1
(
z
)
)
=
0
,
j
=
1,2
,
…
,
I
-
1
,
1
Δ
t
(
u
h
,
j
-
1
/
2
n
+
1
(
z
)
-
u
h
,
j
-
1
/
2
n
(
z
)
)
+
2
h
c
j
-
1
/
2
n
+
1
(
1
2
u
h
,
j
n
+
1
(
z
)
-
1
2
u
h
,
j
-
1
n
+
1
(
z
)
)
-
2
h
2
(
2
a
j
-
1
/
4
n
+
1
(
z
)
u
h
,
j
n
+
1
(
z
)
hhhhhh
-
2
(
a
j
-
1
/
4
n
+
1
(
z
)
+
a
j
-
3
/
4
n
+
1
(
z
)
)
u
h
,
j
-
1
/
2
n
+
1
(
z
)
hhhhhh
+
2
a
j
-
3
/
4
n
+
1
(
z
)
u
h
,
j
-
1
n
+
1
(
z
)
)
=
0
,
j
=
1,2
,
…
,
I
,
where Δ
t is time step size, u
h
,
j
n
=
u
h
(
x
j
,
n
Δ
t
), and the transport velocity will be fixed as c
(
x
)
=
1.We perform a simulation atT
=
2 with Δ
t
/
h
2
=
1 for the comparison of exact solution and numerical solution of (48), which are presented in Figure 9, with M
=
500 and ν
=
0.5. From it we can see that numerical solutions are consistent well with the exact solution, and the numerical solution obtained by QMC is more accurate than that of MC. Figures 10 and 11 correspond to the numerical approximations of QoI with different choice of ν. It is easy to find that the fluctuation of QMC is less obvious than that of MC, regardless of the value of ν. Figure 12 shows error estimates of MC and QMC when ν is chosen as 0.5. We can see from it that QMC presents a faster convergence rate than MC. It can be concluded that, to get the same results, we only need a smaller number of QMC samples than that of MC, which increases the computation efficiency and is good news for stochastic computation.Figure 9
Approximate the exact solution using MC method and QMC method, respectively, withν
=
0.5 and M
=
500.Figure 10
MC approximation of QoI with different choice ofν.Figure 11
QMC approximation of QoI with different choice ofν.Figure 12
Error estimate of MC (solid line) and QMC (dashed line), withν
=
0.5.


## 4.1. Steady State Diffusion Problem

We first consider the following problem in one spatial dimension (d
=
1) and N
(
=
20
) random dimensions:
(42)
-
∂
∂
x
[
a
(
x
,
z
)
∂
u
∂
x
(
x
,
z
)
]
=
cos
⁡
⁡
(
2
π
x
)
,
x
∈
(
0,1
)
,
z
∈
R
N
,
(43)
u
(
0
,
·
)
=
u
(
1
,
·
)
=
0
.The region of interest[
0,1
] can be decomposed into a uniform grid T
h with nodes 0
=
x
0
<
x
1
<
⋯
<
x
I
=
1 and h
=
x
i
-
x
i
-
1
=
1
/
I, i
=
1
,
2
,
…
,
I. To obtain quadratic basis functions, we will use the midpoints x
i
-
1
/
2
=
(
x
i
+
x
i
-
1
)
/
2 as the interpolation points. Accordingly we place a dual grid T
h
* with nodes 0
=
x
0
<
x
1
/
4
<
x
3
/
4
<
⋯
<
x
I
-
3
/
4
<
x
I
-
1
/
4
<
x
I
=
1, where x
i
-
k
/
4
=
x
i
-
(
k
/
4
)
h, (k
=
1
,
3; i
=
1
,
2
,
…
,
I). Select trial function space U
h as the quadratic element space of Lagrange type with respect to T
h, and the test function space V
h corresponding to T
h
* is taken as the piecewise constant function space. Then the numerical solution u
h for (42) can be uniquely written as
(44)
u
h
(
x
,
z
)
=
∑
i
=
1
I
[
u
h
,
i
(
z
)
ϕ
i
(
x
)
+
u
h
,
i
-
1
/
2
(
z
)
ϕ
i
-
1
/
2
(
x
)
]
,
which is a piecewise quadratic FVEM approximation. We obtain the discrete schemes:
(45)
-
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
-
1
/
2
(
z
)
+
2
a
(
x
j
-
1
/
4
,
z
)
+
2
a
(
x
j
+
1
/
4
,
z
)
h
u
h
,
j
(
z
)
-
2
a
(
x
j
+
1
/
4
,
z
)
h
u
h
,
j
+
1
/
2
(
z
)
=
∫
x
j
-
1
/
4
x
j
+
1
/
4
cos
⁡
⁡
(
2
π
x
)
d
x
,
j
=
1,2
,
…
,
I
-
1
,
-
2
a
(
x
j
-
3
/
4
,
z
)
h
u
h
,
j
-
1
(
z
)
+
2
a
(
x
j
-
3
/
4
,
z
)
+
2
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
-
1
/
2
(
z
)
-
2
a
(
x
j
-
1
/
4
,
z
)
h
u
h
,
j
(
z
)
=
∫
x
j
-
3
/
4
x
j
-
1
/
4
‍
cos
⁡
⁡
(
2
π
x
)
d
x
,
j
=
1,2
,
…
,
I
,
and they yield a tridiagonal linear system for the nodal values:
(46)
A
(
z
)
U
h
(
z
)
=
F
.The value of QoI:Q
(
u
(
z
)
)
=
∫
0
1
‍
u
(
x
,
z
)
d
x can be computed exactly by a Simpson formula, yielding
(47)
Q
(
u
(
z
)
)
=
∑
i
=
1
I
h
3
(
u
h
,
i
-
1
/
2
+
u
h
,
i
+
u
h
,
i
+
1
/
2
)
.The numerical approximations of QoI using both MC method and QMC method are presented in Figures5 and 6 with different choice of ν. We can see from both of them that when the number of samples becomes larger and larger, the approximations gradually tend to a fixed value. Besides, approximations using QMC method in Figure 6 are more smooth than that in Figure 5, which means QMC method is more accurate than MC method regardless of the value of ν.Figure 5
MC approximation of QoI with different choice ofν.Figure 6
QMC approximation of QoI with different choice ofν.Figures7 and 8 show error estimates of MC method and QMC method when ν is chosen as 1.5 and 2.5, respectively. It is obvious that the QMC error is smaller than the MC error. Their slopes are also estimated as a comparison, and we can see that the slopes of QMC are bigger than that of MC. That means, with the same number of samples, estimate of QMC with Sobol sequence converges faster than that of MC method.Figure 7
Error estimate of MC (dotted line) and QMC (solid line) withν
=
1.5.Figure 8
Error estimate of MC (dotted line) and QMC (solid line) withν
=
2.5.


## 4.2. Convection-Diffusion Problem

Finally we apply the MCFVEM algorithm to the stochastic convection-diffusion equation:(48)
∂
u
∂
t
+
c
(
x
)
∂
u
∂
x
-
∂
∂
x
[
a
(
x
,
z
)
∂
u
∂
x
]
=
0
,
hhhhh
x
∈
(
0,1
)
,
z
∈
R
N
,
t
∈
[
0
,
T
]
,
(49)
u
(
0
,
t
,
·
)
=
1
,
u
(
1
,
t
,
·
)
=
0
,
u
(
x
,
0
,
·
)
=
u
0
(
x
,
·
)
.
The exact solution of (48) is
(50)
u
(
x
,
t
,
z
)
=
1
2
{

[
x
+
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
)
erfc
⁡
[
x
-
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
]
hh
+
exp
⁡
⁡
(
c
x
a
(
x
,
z
)
)
erfc
⁡
[
x
+
c
t
(
4
a
(
x
,
z
)
t
)
1
/
2
]
}
.
In this case we are interested in computing the expected value of the following QoI:
(51)
Q
(
u
(
x
,
T
,
ω
)
)
=
∫
0
1
u
(
x
,
T
,
ω
)
d
x
.
It can be computed by a Simpson formula similarly as in (47).Use the same physical space as in Section4.1, and employ implicit Euler scheme to discretize the first-order time derivative. We obtain the following fully discrete quadratic finite volume element schemes:
(52)
1
Δ
t
(
u
h
,
j
n
+
1
(
z
)
-
u
h
,
j
n
(
z
)
)
+
2
h
c
j
n
+
1
(
1
8
u
h
,
j
-
1
n
+
1
(
z
)
-
3
4
u
h
,
j
-
1
/
2
n
+
1
(
z
)
hhhhhhhhh
+
3
4
u
h
,
j
+
1
/
2
n
+
1
(
z
)
-
1
8
u
h
,
j
+
1
n
+
1
(
z
)
)
-
2
h
2
(
2
a
j
+
1
/
4
n
+
1
(
z
)
u
h
,
j
+
1
/
2
n
+
1
(
z
)
hhhhhh
-
2
(
a
j
+
1
/
4
n
+
1
(
z
)
+
a
j
-
1
/
4
n
+
1
(
z
)
)
u
h
,
j
n
+
1
(
z
)
hhhhhh
+
2
a
j
-
1
/
4
n
+
1
(
z
)
u
h
,
j
-
1
/
2
n
+
1
(
z
)
)
=
0
,
j
=
1,2
,
…
,
I
-
1
,
1
Δ
t
(
u
h
,
j
-
1
/
2
n
+
1
(
z
)
-
u
h
,
j
-
1
/
2
n
(
z
)
)
+
2
h
c
j
-
1
/
2
n
+
1
(
1
2
u
h
,
j
n
+
1
(
z
)
-
1
2
u
h
,
j
-
1
n
+
1
(
z
)
)
-
2
h
2
(
2
a
j
-
1
/
4
n
+
1
(
z
)
u
h
,
j
n
+
1
(
z
)
hhhhhh
-
2
(
a
j
-
1
/
4
n
+
1
(
z
)
+
a
j
-
3
/
4
n
+
1
(
z
)
)
u
h
,
j
-
1
/
2
n
+
1
(
z
)
hhhhhh
+
2
a
j
-
3
/
4
n
+
1
(
z
)
u
h
,
j
-
1
n
+
1
(
z
)
)
=
0
,
j
=
1,2
,
…
,
I
,
where Δ
t is time step size, u
h
,
j
n
=
u
h
(
x
j
,
n
Δ
t
), and the transport velocity will be fixed as c
(
x
)
=
1.We perform a simulation atT
=
2 with Δ
t
/
h
2
=
1 for the comparison of exact solution and numerical solution of (48), which are presented in Figure 9, with M
=
500 and ν
=
0.5. From it we can see that numerical solutions are consistent well with the exact solution, and the numerical solution obtained by QMC is more accurate than that of MC. Figures 10 and 11 correspond to the numerical approximations of QoI with different choice of ν. It is easy to find that the fluctuation of QMC is less obvious than that of MC, regardless of the value of ν. Figure 12 shows error estimates of MC and QMC when ν is chosen as 0.5. We can see from it that QMC presents a faster convergence rate than MC. It can be concluded that, to get the same results, we only need a smaller number of QMC samples than that of MC, which increases the computation efficiency and is good news for stochastic computation.Figure 9
Approximate the exact solution using MC method and QMC method, respectively, withν
=
0.5 and M
=
500.Figure 10
MC approximation of QoI with different choice ofν.Figure 11
QMC approximation of QoI with different choice ofν.Figure 12
Error estimate of MC (solid line) and QMC (dashed line), withν
=
0.5.


## 5. Conclusions

The paper utilizes the Karhunen-Loève expansion of the coefficient of convection-diffusion equation and transforms the initial stochastic problem into a deterministic one albeit with a parameter in high dimensions. We design finite volume element method in the physical space, which has reasonable accuracy and high efficiency, and we employ Monte Carlo method in the random space. Monte Carlo methods are both general and simple to code, and they are naturally suited for parallelization when the deterministic code already exists. A quasi-Monte Carlo technique is also used to accelerate the Monte Carlo method. Numerical experiments indicate that these methods are effective and they can be extended to spatially multidimension stochastic partial differential equation directly.



---

*Source: 642470-2014-05-22.xml*