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

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 theMonte 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.


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][2][3][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 convectiondiffusion equation with homogeneous Dirichlet boundary conditions.Let  be a convex bounded polygonal domain in R  ,  = 1, 2, 3, and let (Ω, F, ) be a complete probability space.Here Ω is the set of outcomes, F ∈ 2 Ω is the -algebra of events, and  : F → [0, 1] is a probability measure.Consider the following stochastic linear convection-diffusion problem: find a stochastic function,  :  × Ω → R, such that the following equation holds: (x, , ) = 0,  × (0, ] × Ω. Here x = ( 1 ,  2 , . ..,   ) ∈  ⊂ R  ,  = 1, 2, 3, are spatial coordinates, :  × Ω → R is a stochastic function with continuous and bounded covariance function, and  is the velocity of flow.Our goal is to compute the quantity of interest (QoI): 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 Section 2. 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.

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 (x, ) in ( 1) is a stationary random field with continuous covariance function Cov[] :  ×  → R. For well-posedness we require (x, ) ≥  min > 0. We expand (x, ) using truncated KL expansion with  terms where [(x, )] is the mean of the input process (x, ), {  } is a set of independent and identically distributed random variables, and   ,   denote the eigenvalues and eigenfunctions of the correlation function  with respect to the eigenvalue problem: The summation in ( 5) is truncated at finite number , and the number of terms  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, (x, , ) = (x, ,  1 (),  2 (), . . .,   ()).We denote by Γ  the bounded image set of the random variables   (Ω), Γ = Γ 1 × Γ 2 × ⋅ ⋅ ⋅ Γ  , and we assume that the random vector z has a joint probability density function  : Γ → R + that factorizes as (z) = ∏  =1   (  ), for all z ∈ Γ ⊂ R  .Then we can rewrite (1)-(3) as follows with an -dimensional parameter:

Finite Volume Element Approximation in the Physical
Space.The physical domain  can be decomposed into a grid  ℎ with a set of spaced nodes; accordingly we place dual grid  * ℎ [11].Select trial function space  ℎ () as the finite element space of piecewise continuous polynomials on  ℎ , with dimension  = dim( ℎ ()) and basis {  }  =1 .The test function space  ℎ () corresponding to  * ℎ is taken as the piecewise constant function space, with the same dimension as  ℎ () and basis {  }  =1 .We introduce the space  2  (Γ) of square integrable functions on Γ with respect to the measure (z)z; that is, Problem ( 8) admits a unique solution for every z ∈ Γ.Then the semidiscrete solution can be expressed as Denoting by (, z) = [ ℎ,1 (, z),  ℎ,2 (, z), . . .,  ℎ, (, z)] the vector of nodal values as functions of the random variables z and , problem (8) can be written in algebraic form as where  is deterministic matrix with respect to (x), matrix (z) is a stochastic stiffness matrix with respect to the random diffusion coefficient (x, z) which is an affine function of the random variables   , and   = ∫  (x, )  (x)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.

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 MCFVEM
Step 1. Give a number of realizations, , trial function space  ℎ (), and test function space  ℎ (), as defined in Section 2.2.
Step 3. Finally, use the sample average

MCFVEM Error Analysis.
After the PDE is discretized by finite volume element method, we can compute discrete solutions  ℎ (z(  )),  = 1, 2, . . ., .Our goal is to compute [((z))], and then our estimator is The computational error naturally separates into two parts: The size of the spatial discretization schemes controls the space discretization error  ℎ , while the number of realizations,  of  ℎ , controls the statistical error   .

Discretization Error
Assumption 1.The functional  :  → R, with (0) = 0, is globally Lipchitz; that is, Assumption 2. There exist  > 0 and   (z) > 0, with We apply Assumptions 1 and 2 to estimate the discretization error where where  is (0, 1) and ⇀ denotes convergence of the distributions, also called weak convergence; that is, the convergence (17) means that the following limit holds as  → ∞: for all bounded and continuous function .
Proof.See [ Let us analyze the statistics error assuming that  is a globally Lipschitz function: and we can estimate where where  is (0, 1).Motivated by central limit theorem (CLT), that is, (28) The exact variance  is in practice approximated by the sample variance: Thus, we have the following CLT error bound: where  0 is an appropriate constant.

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 lowdiscrepancy 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.
Figure 1 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.
We calculate the quantity   () = ∫ [0,1]  (x)x using an equal weight cubature formula of the form where x is a sequence of  random points in the  dimensional cube [0, 1]  , and  : [0, 1]  → R is a continuous function.
To estimate the error in a QMC integration, a shifted averaged lattice rule [16][17][18][19][20] will be employed.Let  ∼ ([0, 1]  ) be a uniformly distributed random vector (shift) and  Estimate of QMC error based on the -sample standard deviation, that is, In the following, we will use MC and QMC to numerically evaluate the integrals.
Model 2. Consider the following: with   = 4.3/,  1 = /4, and  2 = /5.The exact solution is 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 Figures 2 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.

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 process (x, ) in (1) strictly positive, we assume (x, ) =  (x,) , where (x, ) is a stationary random field with the Matérn covariance function: where Γ is the gamma function and  ] is the modified Bessel function of the second kind.The following special cases of  are considered:  where  scales the correlation length,  denotes the distance between two spatial points, and  is a constant.In the following, we choose  = 2 and  = 0.1.
We expand (x, ) using truncated KL expansion with  terms: where {  } 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  we solve the eigenvalue problem: As no analytical solutions are available for the eigenvalue problem (41), a Galerkin numerical eigenvalue solver is employed [21].We solve (41) by discretizing  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][23][24].
The spectral contribution of  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  = 20.

+ 2𝑎 (𝑥
and they yield a tridiagonal linear system for the nodal values: The value of QoI: ((z)) = ∫ The numerical approximations of QoI using both MC method and QMC method are presented in Figures 5 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 ].
Figures 7 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.
It can be computed by a Simpson formula similarly as in (47).
Use the same physical space as in Section 4.1, and employ implicit Euler scheme to discretize the first-order time derivative.We obtain the following fully discrete quadratic finite volume element schemes:

Mathematical Problems in Engineering 5 2000Figure 1 :
Figure 1: Two-dimensional scatter plot for MC points (a) and QMC points (b).

Figure 2 :
Figure 2: The exact error and CLT error estimation for the integral of Model 1, with  = 10.

Figure 3 :Figure 4 :
Figure 3: The exact error and CLT error estimation for the integral of Model 2, with  = 10.

Figure 5 :
Figure 5: MC approximation of QoI with different choice of ].

Figure 6 :
Figure 6: QMC approximation of QoI with different choice of ].