Mixed Weak Galerkin Method for Heat Equation with Random Initial Condition

(is paper is devoted to the numerical analysis of weak Galerkin mixed finite element method (WGMFEM) for solving a heat equation with random initial condition. To set up the finite element spaces, we choose piecewise continuous polynomial functions of degree j + 1 with j≥ 0 for the primary variables and piecewise discontinuous vector-valued polynomial functions of degree j for the flux ones. We further establish the stability analysis of both semidiscrete and fully discrete WGMFE schemes. In addition, we prove the optimal order convergence estimates in L2 norm for scalar solutions and triple-bar norm for vector solutions and statistical variance-type convergence estimates. Ultimately, we provide a few numerical experiments to illustrate the efficiency of the proposed schemes and theoretical analysis.


Introduction
In this paper, we apply the recently developed weak Galerkin mixed finite element method to solve the following heat equation with random initial condition: u(x, t, ω) � 0, for(x, t, ω) ∈ zD ×(0, T) × Ω, where D is an open-bounded polygonal or polyhedral domain in R 2 or R 3 with boundary zD, (Ω, Σ, P) is a probability space, f is a given deterministic function, and u 0 is a random initial input. In this paper, we shall concentrate on two-dimensional problem and the result can be applied to the three-dimensional case. An alternative choice for solving (1)-(3) is to introduce a vector-valued random process q � − ∇u in order that the problem (1)-(3) is converted into a system of two first-order linear equations, that is to say, to seek a vector-valued random process q and a scalar random process u such that q + ∇u � 0, for(x, t, ω) ∈ D ×(0, T) × Ω, u(x, 0, ω) � u 0 (x, ω), for(x, ω) ∈ D × Ω, u(x, t, ω) � 0, for(x, t, ω) ∈ zD ×(0, T) × Ω.
e parabolic problems (1)- (3) have been analyzed by using various numerical schemes. In [1], the stability criterion for the finite element solutions is given and it is demonstrated that the stability criterion is valid for the practical calculation of statistical quantities of the finite element solutions. Tasaka [2] presents an order-of-convergence result for a finite element solution of the stochastic heat equation with random initial value and boundary conditions. In [3], probabilistic convergence of random Galerkin approximations of a heat equation with random initial conditions is studied, and both the continuous and discrete time (backward difference) approximations are considered. e authors study the one-dimensional homogeneous heat equation with random initial conditions and consider the minimal contrast estimates of parameters of random fields which are normalized solutions in [4]. Kozachenko and Veresh [5] study the conditions for the applicability of the Fourier method to parabolic equations with random initial conditions from Orlicz spaces of random variables and obtain estimates for the distribution of the supremum. In [6], the author constructs a solution of a parabolic stochastic partial differential equation with random initial conditions by using Kolmogorov's criterion. Additionally, by virtue of the weak regularity of solutions of stochastic problems, the idea of using discontinuous elements to approximate true solutions of stochastic problems is considered in [7]. In view of this, the goal of this paper is to extend the recently developed weak Galerkin (WG for short) mixed finite element method to approximate equations (1)- (3). e WG method is first introduced in [8] for the secondorder elliptic problem based on a definition of discrete weak gradient operator arising from RT element [9] and BDM element [10]. Its central idea is to approximate the differential operators by weak forms as distributions over the space of discontinuous functions including boundary information.
en, the WG method is successfully applied to elliptic interface problems [11,12] and linear parabolic problems [13][14][15][16] and further developed for other applications, such as biharmonic problems [17][18][19][20][21][22], Cahn-Hilliard problems [23], and Stokes problems [24][25][26][27][28][29]. To ensure the method is highly flexible in element construction and mesh generation, the idea of a stabilization term is introduced in [30], which allows for arbitrary piecewise polynomial shape functions in deformed and honeycomb meshes. e use of parameter-free stabilization term is the key point to the high flexibility of WG to ensure weak continuities for approximating functions. Similarly, the WGMFEM with a stabilization term is introduced in [31] based on the definition of discrete weak divergence operator and applicable for general finite element partitions consisting of regular polygon shapes in 2D or polyhedra in 3D. Moreover, WGMFEM is developed in second-order elliptic equations with Robin boundary conditions [32], parabolic differential equations with memory [33], Helmholtz equations with large wave numbers [34], and biharmonic equations [35]. In addition to solving deterministic problems, there has also been some progress in the WG method for stochastic problems, such as stochastic Brinkman problems [36], elliptic problems with stochastic jump coefficients [37], and random interface grating problems [38].
is paper is organized as follows. In Section 2, we present some necessary notations and definitions and then establish the semidiscrete and fully discrete WGMFE schemes for parabolic problem with random initial conditions (4)-(7). In Section 3, we first provide the stability analysis for the semidiscrete and fully discrete WGMFE schemes. And then, we derive the optimal order convergence estimates for the scalar solution in L 2 norm and the vectorvalued solution in triple-bar norm. Finally, we study the statistical variance-type convergence estimates for the scalar and vector-valued solutions, respectively. Ultimately in Section 4, we provide a few numerical experiments to illustrate the behavior of the proposed schemes and confirm our theoretical analysis.

Weak Galerkin Mixed Finite Element Method
In this section, we firstly introduce some notations and definitions of spaces and corresponding norms for describing the vector-valued random process q and scalar random process u. en, we set up a pair of WG finite element spaces consisting of weak functions which are discontinuous and defined in a piecewise manner on the interior of any partition element and its boundary. We also define a weak divergence operator on such weak functions as an extension of the standard divergence operator defined on smooth functions. While we use such weak divergence operator to discretize equations (4)-(7), we recommend a stabilization term to patch the different parts of the weak functions on different elements, which also ensures the uniqueness of the WG approximations. Finally, we present the semidiscrete and fully discrete WGMFE schemes.

Notations and Variational Formulation.
Let H � L 2 (D) be the standard square integral space with norm ‖·‖ and inner product (·, ·) and H m be the Sobolev space with norm ‖·‖ m for m ≥ 1. Furthermore, denote by H m 0 the subspace of H m whose element vanishes at the boundary. Additionally, we define a space with a norm ‖v‖ H(div;D) � ‖v‖ 2 +‖∇ · v‖ 2 1/2 .

Weak Galerkin Approximations.
ere are many numerical methods to solve equations (1)-(3), such as stochastic collocation methods [39], discontinuous Galerkin methods [40], and finite element methods [2]. e standard finite element approach is to properly choose two finite element spaces W h ⊆ H 1 0 and V h ⊆ H(div; D), respectively, and then project equations (11) and (12) into the space W h × V h . erefore, we obtain finite dimensional equations whose solutions are the approximations 2 Mathematical Problems in Engineering of (11) and (12).
e recently developed weak Galerkin finite element method adopts the piecewisely defined discontinuous functions to construct the finite element spaces W h and V h , respectively, and extend the definition of the divergence operator to such discontinuous functions. Based on these settings, we can set up the weak Galerkin finite element approximations together with an extra stabilization term.
Let D h be a partition of the domain D satisfying the shape regularity requirements A1-A4 in [31], h K be the partition diameter of the element K ∈ D h , and h � max K h K . For any element K ∈ D h , let K 0 and zK denote its interior and boundary, respectively. Denote by E h the set of all edges in D h , and let E 0 h � E h \zD be the set of all interior edges.
Denote by n e a fixed unit normal vector on E h satisfying when e ∈ E 0 h , n e is normal to e, and when e ∈ E h ∩ zD, n e is outward normal to e.
For each element K ∈ D h , let P j (K) denote the set of polynomials defined on K with degree no more than j ≥ 0. Similarly, we can define a polynomial space P j (e), for each e ∈ E h . Now, we define three weak Galerkin finite element spaces by Remark 1. e functions in V h are totally discontinuous. Such functions may contain independent information on any K ∈ D h and any e ∈ E h . Next, we extend the divergence operator where n is a unit outward normal vector to zK. Remark 2. Let K 1 and K 2 be two adjacent elements in D h sharing a common interior edge e. n 1 and n 2 are unit outward normal vectors of K 1 and K 2 to e, respectively. Let n e be the pre-fixed unit normal vector to e. As can be seen from Figure 1, we shall obtain n e � n 1 and n e � − n 2 , or vice versa.
In order to simplify the notations, we introduce some bilinear forms: We also define two norms where E denotes the expectation with respect to the probability measure P. From [30,31], we know that the bilinear form a s (·, ·) is bounded, symmetric, and positive definite in V h , and the bilinear form b(·, ·) is bounded in V h × W h and they also satisfy the standard inf-sup condition. Now, we can formulate the semidiscrete weak Galerkin mixed finite element scheme. For any t ∈ (0, T), finding (19) and satisfying the initial condition where P h is an L 2 projection operator from H 1 0 onto W 0 h . According to [30,31], for almost every ω ∈ Ω, the solutions of problems (18)- (20) exist and are unique. Now, we turn our attention to the fully discrete numerical scheme. We introduce a time step size κ � T/N for some positive integer N and t n � nκ for n � 0, 1, . . . , N. By q n and u n , we denote the approximation of q(t n ) and u(t n ), respectively. We use the backward Euler method to approximate the time derivative in (19); then, the fully discrete scheme reads seeking (q n , u n ) ∈ L 2 (Ω; Figure 1: Relationship between n e and n.

Mathematical Problems in Engineering
and satisfying the initial condition where z u n � (u n − u n− 1 )/κ.

Numerical Analysis
In this section, we first present the stability analysis for the semidiscrete and fully discrete WGMFE schemes. And then, we derive an optimal order convergence estimate for the scalar solution in L 2 norm ‖·‖ L 2 (Ω;H) . We also analyze the convergence order for the vector-valued solution in the norm |‖·‖| L 2 (Ω;V h ) . Finally, we study the statistical variancetype convergence estimates for both the scalar and vectorvalued solutions, respectively.

Stability Analysis.
In this section, we shall derive stability estimates in both ‖·‖ L 2 (Ω;H) and |‖·‖| L 2 (Ω;V h ) norms for the semidiscrete and fully discrete WGMFE schemes. ese will be the bases of our later error analysis.

Theorem 1.
For the numerical solutions of (18) and (19) with the random initial setting (20), the stability holds true: where C > 0 is a constant independent of the mesh size h.
Proof. In equations (18) and (19) Adding (26) to (27), we have Because of Cauchy-Schwarz inequality and |‖q h (ω)‖| 2 ≥ 0, we obtain Hence, Squaring both sides of (30) and using Cauchy-Schwarz inequality, where C is a constant independent of the mesh size h. Taking expectation on both sides of (31), then (24) is obtained. We differentiate (18) with respect to t, Taking v h � q h and χ h � u h,t in (32) and (19) and summing up, we have 1 2 Noticing us, Taking expectation on both sides of (35), and combining with the definitions (16) and (17), we obtain (25), which completes the proof.

□
Similarly, the fully discrete analysis is given below.

Theorem 2.
For the WGMFE solutions of (21) and (22) with the initial condition (23), the stability follows where C > 0 is a constant independent of the mesh size. (21) and (22) and sum up these two equations, we have Making use of Cauchy-Schwarz inequality and |‖q n (ω)‖| 2 ≥ 0, we obtain and by recurrence, we obtain Squaring and taking expectation on both sides of (40), we derive (36).
Using the backward Euler method to discretize the time derivative in (32), Taking v h � q n and χ h � z u n in (41) and (22) and summing up, z u n (ω)‖ 2 + a s z q n , q n ) � f t n , z u n (ω)).
Since a s z q n , q n ) � 1 2 z a s q n , q n + κ 2 a s z q n , z q n ), we find By virtue of Cauchy-Schwarz inequality and the iteration method, we obtain Taking expectation on both sides of (45), and making use of the definitions (16) and (17), the proof is completed. □

Optimal Order Error Estimates of WGMFEM.
In this section, we shall establish convergence analysis for both semidiscrete and fully discrete weak Galerkin methods, respectively. Specifically, we firstly study the optimal order convergence estimates for the scalar solution. And then, we derive the convergence order for the vector-valued solution q h . Finally, we analyze the statistical variance-type convergence estimates for both the scalar and vector-valued solutions, respectively. First, we introduce two new projection operators into the finite element spaces V h and W h by using local L 2 projections. For each element K ∈ D h , denote by Q 0 the projection operator from (L 2 (K)) 2 onto (P j (K)) 2 . Similarly, for each edge e ∈ E h , let Q b be the projection operator from L 2 (e) onto P j (e). For any

Semidiscrete Error Analysis.
For simplicity, we introduce two new notations: In [41], WGMFEM is applied to study the following deterministic problem: seeking a flux function q and a scalar function u satisfying q + ∇u � 0, and the initial-boundary condition According to the results of [41], we have the following useful estimates.

Mathematical Problems in Engineering
with Var[u] denoting the variance of the random variable u.
Proof. For m h ∈ L 2 (Ω; W 0 h ), taking account of (50), we obtain, for almost ω, By virtue of the definition of m h (t, ω), the error estimate of L 2 projection, and (54), we have Squaring both sides of (55) and utilizing Cauchy-Schwarz inequality, it follows that Take expectation on both sides of (56), which leads to (52).
By virtue of Cauchy-Schwarz inequality, we supply where Cov represents the covariance in regard to the probability measure P. Combining the fact and Cauchy-Schwarz inequality, we provide Additionally, Because of (58)-(60), we find According to (24), (52), the random initial condition (20), and the regularity condition of f, we complete the proof.
□ Before the next theorem, we shall introduce the definition of the trace operator Tr(·): for a square matrix A, Tr(A) is defined as the sum of the elements on the main diagonal of A.

Theorem 4. Under the assumption of eorem 3, we have
Proof. For θ h ∈ L 2 (Ω; V h ), considering (51), we obtain, for almost ω, Because of the definition of θ h (t, ω), Cauchy-Schwarz inequality, the error estimate of L 2 projection, and (63), it follows that Taking expectation on both sides of (64), together with the definitions (16) and (17), the first inequality of (62) is provided.
Next, we will derive the second result of (62). In view of the two-dimensional problem, for any t ∈ (0, T), q ∈ L 2 (Ω; [H j+1 ] 2 ) and q h � q 0 , q b ∈ V h ; assuming that q � (q 1 , q 2 ) T and q 0 � (q 01 , q 02 ) T , then Making use of (57), we get Together with Cauchy-Schwarz inequality, (66), and (67), we obtain Coupled with the first inequality of (62), (25), the definition of |‖ · ‖| L 2 (Ω;V h ) , and the regularity for f, the proof is completed. □ Remark 3. In view of (65) in eorem 4, for two-dimensional vector-valued random process q � (q 1 , q 2 ) T and q 0 � (q 01 , q 02 ) T , it can be seen that the results of Var[q] and Var[q 0 ] are two 2 × 2 square matrixes. In analogy with the deterministic case, we focus more on the information of principle diagonal parts of two square matrixes, i.e., the difference between sum of components variance of true solution q and numerical solution q 0 . is is the reason why we use trace operator here.

Fully Discrete Error Analysis. Define
Similar to the semidiscrete case, we introduce the fully discrete error estimates of (47)-(49) in [41] as follows.

Numerical Experiments
In this section, we will provide two examples to illustrate our theoretical results. Let D � (0, 1) × (0, 1). We take two different probability spaces in our calculation. One is a standard normal distribution N(0, 1), and the other one is a uniform distribution of the interval (0, 1). We use 1000 sample points ω generated by legitimate 64-bit Matlab 2017b to simulate the randomness in terms of these two random variables, respectively. Consider the following equations: We utilize a uniform triangular mesh D h to partition D and conduct numerical experiments with the fully discrete weak Galerkin mixed finite element schemes (21)-(23) by choosing j � 0 in W 0 h and V h . e Monte Carlo sampling method is applied to evaluate the expectation and variance appearing in eorem 5 and eorem 6.
For simplicity, we firstly introduce several abbreviation notations in the following tables for presenting the numerical results: e n u ≔ u t n − u n , e n q ≔ q t n − q n , r n u ≔ Var u t n − Var u n , r n q ≔ Tr Var q t n − Tr Var q n 0 . (80)

(86)
In this test, Tables 5 and 6 illustrate all error estimate results with the use of a normal distribution as the probability space of equations (76)-(79). Likewise, a uniform distribution is utilized as the probability space to perform the numerical simulations, and the corresponding results are shown in Tables 7 and 8, respectively.

Data Availability
e Matlab data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.