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.
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. Let be a convex bounded polygonal domain in , , , , and let () be a complete probability space. Here is the set of outcomes, is the -algebra of events, and is a probability measure. Consider the following stochastic linear convection-diffusion problem: find a stochastic function, , such that the following equation holds: Here , , , , , , , are spatial coordinates, : 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.
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 in (1) is a stationary random field with continuous covariance function . For well-posedness we require . We expand using truncated KL expansion with terms where is the mean of the input process , 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, . We denote by the bounded image set of the random variables , , and we assume that the random vector has a joint probability density function that factorizes as , for all . Then we can rewrite (1)–(3) as follows with an -dimensional parameter:
2.2. 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 and basis . The test function space corresponding to is taken as the piecewise constant function space, with the same dimension as and basis . We introduce the space of square integrable functions on with respect to the measure ; that is, . Then we introduce the semidiscrete problem: find such that for all
Problem (8) admits a unique solution for every . Then the semidiscrete solution can be expressed as . Denoting by the vector of nodal values as functions of the random variables and , problem (8) can be written in algebraic form as where is deterministic matrix with respect to , matrix is a stochastic stiffness matrix with respect to the random diffusion coefficient which is an affine function of the random variables , and 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 MCFVEM
Step 1. Give a number of realizations, , trial function space , and test function space , as defined in Section 2.2.
Step 2. For each , sample independent and identically distributed realizations of the diffusion coefficient based on realizations of , and find a corresponding approximation which satisfies (9).
Step 3. Finally, use the sample average to approximate and to approximate .
2.4. MCFVEM Error Analysis
After the PDE is discretized by finite volume element method, we can compute discrete solutions , . Our goal is to compute , 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 .
2.4.1. Discretization Error
Assumption 1. The functional , with , is globally Lipchitz; that is,
Assumption 2. There exist and , with for some , such that
We apply Assumptions 1 and 2 to estimate the discretization error where 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). Assume , are independent, identically distributed (i.i.d) random variables and . Then where denotes convergence in probability; that is, the convergence (15) means for all .
Theorem 4 (the strong law of large numbers). Assume , are i.i.d random variables with and . Then where denotes almost sure convergence; that is, the convergence (16) means .
Theorem 5 (the central limit theorem). Assume , are i.i.d and , . Then where is 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 [9, page 23].
Theorem 6 (law of the iterated logarithm). Assume , are i.i.d and , . Then
Let us analyze the statistics error assuming that is a globally Lipschitz function: . Then and we can estimate where . We conclude that for Lipschitz functionals we have with constant uniformly bounded with respect to . Since , 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
Let , , , be i.i.d random variables with and . We have
Let , then , are i.i.d, and , . Then by Theorem 5, we have where is .
Motivated by central limit theorem (CLT), that is, with , we can write the statistical error as The exact variance is in practice approximated by the sample variance:
Thus, we have the following CLT error bound: where 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.
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.
(a)
(b)
We calculate the quantity using an equal weight cubature formula of the form where is a sequence of random points in the dimensional cube , and is a continuous function.
To estimate the error in a QMC integration, a shifted averaged lattice rule [16–20] will be employed. Let be a uniformly distributed random vector (shift) and the QMC formula with the random shift . Observe that the shift preserves the lattice structure and , for all . Take i.i.d shifts . Then the estimate of the integral is 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 1. Consider the following: , with , . The exact solution is where denotes the imaginary unit and is the real part of .
Model 2. Consider the following:
with , , and . 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.
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 process in (1) strictly positive, we assume , where 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 and .
We expand using truncated KL expansion with terms: where is a set of i.i.d uniform random variables over 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 is used. For other numerical methods, see [22–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 .
4.1. Steady State Diffusion Problem
We first consider the following problem in one spatial dimension () and random dimensions:
The region of interest can be decomposed into a uniform grid with nodes and , . To obtain quadratic basis functions, we will use the midpoints as the interpolation points. Accordingly we place a dual grid with nodes , where , (; ). Select trial function space as the quadratic element space of Lagrange type with respect to , and the test function space corresponding to is taken as the piecewise constant function space. Then the numerical solution for (42) can be uniquely written as which is a piecewise quadratic FVEM approximation. We obtain the discrete schemes: and they yield a tridiagonal linear system for the nodal values:
The value of QoI: can be computed exactly by a Simpson formula, yielding
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.
4.2. Convection-Diffusion Problem
Finally we apply the MCFVEM algorithm to the stochastic convection-diffusion equation: The exact solution of (48) is In this case we are interested in computing the expected value of the following QoI: 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: where is time step size, , and the transport velocity will be fixed as .
We perform a simulation at with for the comparison of exact solution and numerical solution of (48), which are presented in Figure 9, with and . 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.
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project is supported by the National Natural Science Foundation of China Grant no. 11071123, the Scientific Research Fund of Zhejiang Provincial Education Department (no. Y201327415), and public projects of Lishui Science and Technology Bureau (no. 2013JYZB18).