Convergence Analysis and Cost Estimate of an MLMC-HDG Method for Elliptic PDEs with Random Coefficients

We considered an hybridizable discontinuous Galerkin (HDG) method for discrete elliptic PDEs with random coefficients. By an approach of projection, we obtained the error analysis under the assumption that a(ω, x) is uniformly bounded. Together with the HDG method, we applied a multilevel Monte Carlo (MLMC) method (MLMC-HDG method) to simulate the random elliptic PDEs. We derived the overall convergence rate and total computation cost estimate. Finally, some numerical experiments are presented to confirm the theoretical results.


Introduction
In this paper, we focus on the elliptic partial differential equations (PDEs) with random coefficients which arise in random vibrations, seismic activity and oil reservoir management (see, e.g., [1][2][3][4]). For elliptic stochastic PDEs, there are many numerical methods, such as intrusive and non-intrusive methods. The stochastic collocation method and stochastic Galerkin method are important intrusive methods (see, e.g., [5][6][7][8][9]). The MC method is a non-intrusive method. It is easy to implement, since the deterministic numerical method programs can be used once the samples are given (see, e.g., [10]). However, the MC method has a slow rate of convergence and its computational cost is also very expensive. To improve the computational efficiency, many methods have been proposed. The MLMC method is one of the most important methods. It has been widely used to simulate the stochastic PDEs (see, e.g., [11][12][13][14]).
In [15], the HDG method was originally proposed by Cockburn, Gopalakrishnan and Lazarov. It was then used to solve the second-order elliptic problems with remarkable convergence properties in [16]. The error analysis of elliptic PDEs was given by introducing a projection (see, e.g., [17] and references therein). To date, the HDG method has been successfully applied to parabolic, convection diffusion problems, phase flows and optimal control problem (see, e.g., [18][19][20][21][22]). However, we did not find theoretical or numerical analysis works on the HDG method for PDEs with random coefficients.
In this work, we first considered the HDG method to investigate elliptic PDEs with random coefficients under the assumption of a uniform bound. By a projection approach, we derived the convergence analysis of the HDG method. Then, we achieved rigorous bounds of the error and complexity for the MLMC method based on the HDG method.
The outline of the paper is as follows. In Section 2, the problem setting and HDG method are given. In Section 3, we deduce the error analysis of the velocity and the pressure. In Section 4, we use the error results in Section 3 to derive the complexity analysis of the MLMC method for the elliptic PDEs with random coefficients. In Section 5, some numerical experiments are presented to verify the effectiveness of the proposed method.

Problem Setting
We consider the following elliptic problem with random coefficients: find a random function pair (q, u) : (Ω × D) 2 → R satisfying almost surely (a.s.): where D ⊂ R 2 is a Lipschitz polyhedral domain, (Ω, F , P) is the complete probability space with the set of outcomes Ω, σ-algebra F and probability measure P. The a(ω, x) is the random coefficient, f (ω, x) is the force term, and g(x) is the deterministic Dirichlet boundary condition. Let Y be a random variable in (Ω, F , P). Denote by E[Y] the expected value which is defined by E[Y] = Ω Y(ω)dP(ω). We make the following assumptions to the random coefficient a(ω, x), force term f (ω, x) and boundary condition g(x).

HDG Method
Let D h be a shape regular triangulation of the domain D (see, e.g., [23]). Denote by h K the diameter of simplex K ∈ D h and let h : where E i h is the set of the interior faces, E ∂ h is the set of the boundary faces. For these stochastic Sobolev spaces L 2 (Ω, L 2 (K)), L 2 (Ω, H m (K)), the corresponding norms are defined by Define the following finite dimensional spaces V, W, M by where P k (K) is the set of polynomials of degree at most k on the element K. We write: where (·, ·) is the L 2 inner product on K, and ·, · is the L 2 inner product on ∂K. Given a sample ω ∈ Ω, the HDG numerical method of (1)-(3) is to find (q ω,h , u ω,h ,û ω,h ) ∈ V × W × M such that: for all (v, w, µ) ∈ V × W × M. The numerical traceq ω,h is defined bŷ where τ > 0 is a constant on ∂D h in this paper.

HDG Projection
Denote by (Π V q, Π W u) the HDG projection function, where Π V q and Π W u are components of V and W, respectively. On each simplex K ∈ D h , the projection (Π V q, Π W u) satisfy the following equations: Then, we give an important result for the projection (Π V , Π W ) in the following lemma. The proof is similar to Theorem 2.1 in [17] and hence is omitted. Lemma 1. Give a sample ω ∈ Ω. Suppose the τ is a positive constant on ∂K. Then, the system (8)-(10) is uniquely solvable. Moreover, we have the following estimate for q , u ∈ [0, k] where C is a generic constant independent of h K .
where P M is the L 2 -orthogonal projection onto M which satisfies: Then, we have the following error equations.

Flux Estimate
and the assumption A1 holds. We have: where the constant C is independent of h.
Proof. Give a sample ω ∈ Ω. (14) and µ = −ρû ω,h in (15) and adding these equations together, we obtain: Applying the integration by parts to the left-hand-side (LHS) terms of the above equation, we have where the last identity follows from (16). Hence: Considering the assumption A1 and integrating in the probability space, we obtain: By Cauchy-Schwarz inequality and the fact τ > 0, we have: The triangle inequality implies: By Lemma 1, we complete the proof.

Pressure Estimate
In this subsection, we use the Aubin-Nitsche duality argument to give the error estimate of u − u h L 2 (Ω;D h ) . We introduce a dual problem for any given Θ ∈ L 2 (Ω; L 2 (D)) by finding Φ(ω, x) such that: We assume that the boundary value problem admits the regularity estimate: where C reg depends on σ − and σ + . Then, we have the following theorem.
and the assumption A1 holds. We have: where the constant C is independent of h.
Proof. Give a sample ω ∈ Ω. Considering the properties of (8) and (10), we have: (12), we obtain: where the last equation is due to the fact that Φ(ω, x) · n is continuity and ρû ω,h = 0 on ∂D.
Applying Lemma 1 and Theorem 1, we have: The regularity assumption (20) implies: Hence: Together with the Lemma 1, we obtain: We complete the proof.

The MC Method
Let u ∈ L 2 (Ω; L 2 (D h )) be a random field. The expectation E[u] is approximated by the sample average E M [u], which is defined by where u ω i := u(ω i , ·), i = 1, . . . , M are independent identically distributed (i.i.d.) realizations of the random field u. We give statistical error of the E M [u] in the following lemma.

Lemma 3.
Let u ∈ L 2 (Ω; L 2 (D h )). Then, we have for any M ∈ N: Proof. Considering that u ω i := u(ω i , ·), i = 1, . . . , M are i.i.d. samples of the random field u, we have: Taking the square root of both sides of the inequality, we complete the proof.
In practice, it is difficult to take samples from the random field u, since we do not know it most of time. To overcome this difficulty, we choose samples from the HDG approximation u h L with h L = 1/(2 2L ), where L ∈ N. We define the the classical MC estimator: where u ω i ,h L := u h L (ω i , ·), i = 1, . . . , M are i.i.d. realizations of the random field u h L .

Theorem 3.
Suppose assumptions A1-A2 hold, then: where the constant C is independent of h L and M.
Proof. By the triangle inequality, we have: Using the Cauchy-Schwarz inequality, we have: where the third inequality follows from Theorem 2.
applying Lemma 3, we obtain: Hence: Similarly, we can obtain: The proof is complete.
Theorem 3 implies that the total error result can be decomposed into two parts: statistical error with order M −1/2 and discretization error with order h k+1 L . For a fixed error, the optimal number of samples M should be equilibrated with the mesh size h L , meaning that: where N L = 2 2L is the degrees of freedom of the HDG method. Hence, the total computational cost is:

The MLMC Method
For the MLMC method, the random field u h L can be written as where u h 0 = 0, l = 0, . . . , L. The linearity of the expectation operator yields: We approximate E[u h l − u h l−1 ] by the MC estimator with M l i.i.d. samples on sub-level l. Hence, we estimate E[u] by where the samples over all levels are independent of each other. Theorem 4. Let assumptions A1-A2 hold, then: where the constant C is independent of h l and M, l = 1, . . . , L.
Proof. By the triangle inequality, we obtain: For estimating I, similarly to the proof of Theorem 2, we have: For I I, using the triangle inequality, we obtain: Combing the estimates of I and I I, we obtain the desired result: Analogously, we have: The proof is complete.
Then, we give the error bounds of MLMC-HDG approximation for any distribution {M l } L l=1 over the mesh levels. Just like the single level MC-HDG approximation, we are interested in the optimal ratio of sample size versus grid size in every level. To achieve the overall convergence rate, how to select the sampling number M l is given in the following theorem.
Theorem 5. Suppose assumptions A1-A2 hold. If the sampling number: on sub-level l, l = 1, 2, . . . , L, then we have the error bounds: where where the constant C depends on but is independent of level L.
Proof. We only give the detailed proof of: is similar.
To obtain the overall convergence rate O(h k+1 L ), we choose: for some > 0. Using Theorem 4, we obtain: Similar to the standard finite element, the computational cost of the HDG method is of linear complexity in the number N l on every sub-level l (see, e.g., [24]). For M l samples, the cost is O(M l · N l ). Hence, we have the following bound of overall cost at level L: We complete the proof.

Numerical Experiments
To conform the results of the approximation error and computational cost given in Theorem 5, we present a numerical example. The diffusion coefficient and exact solution are selected as follows: where ω obeys a uniform distribution on [0, 1] and x = (x 1 , x 2 ) ∈ [0, 1] 2 . The homogeneous Dirichlet boundary condition and source term are chosen to match the exact solution. The expectation of the solution is E[u] = 17 2 sin(2πx 1 ) sin(2πx 2 ). We take k = 1 and apply the MLMC-HDG method to simulate the example. All numerical experiments are carried out by using MATLAB R2018b software on an Intel Core i5 machine with 8 GB of memory.
We choose the number of samples with M l = O(l (2+2 ) 2 4(L−l) ) on every sub-level l. The convergence rate of the MLMC-HDG approximation depends on the sub-level l, which is displayed in Figure 1. Figure 2 shows the total CPU-time of E L (u) and E L (q) for different levels L. The CPU-time of sub-level l is plotted in Figure 3, l = 1, . . . , L.
We can see in Figure 1 that E u − u h L 2 (D h ) and E q − q h L 2 (D h ) converge with h 2 . We can observe that the computational cost is Cost(L) ≤ N 2 L in Figure 2. It shows that the CPU-time is O(l 2 2 −4l ) on each level l from Figure 3.
The study of HDG methods is active, since it can capture the discontinuity property well and keep the conservation property. The MLMC method is a good variance reduction technique. The combination of HDG and MLMC methods is a good choice for numerical approximation on PDEs with random coefficients. The MLMC-HDG method can be used in the numerical simulation of evolution equation with random diffusion coefficients. Furthermore, the combination of HDG and multilevel quasi-Monte Carlo methods is also a good topic for stochastic PDEs.

Informed Consent Statement: Not applicable.
Data Availability Statement: The study did not report any data.