# Mixed Weak Galerkin Method for Heat Equation with Random Initial Condition

Authors: Chenguang Zhou; Yongkui Zou; Shimin Chai; Fengshan Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This 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 degreej+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.


---

## Body

## 1. 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:(1)ut−Δu=f,forx,t,ω∈D×0,T×Ω,(2)ux,0,ω=u0x,ω,forx,ω∈D×Ω,(3)ux,t,ω=0,forx,t,ω∈∂D×0,T×Ω,where D is an open-bounded polygonal or polyhedral domain in ℝ2 or ℝ3 with boundary ∂D, Ω,Σ,ℙ is a probability space, f is a given deterministic function, and u0 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(4)q+∇u=0,forx,t,ω∈D×0,T×Ω,(5)ut+∇⋅q=f,forx,t,ω∈D×0,T×Ω,(6)ux,0,ω=u0x,ω,forx,ω∈D×Ω,(7)ux,t,ω=0,forx,t,ω∈∂D×0,T×Ω.The 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. The 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).The WG method is first introduced in [8] for the second-order 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. Then, the WG method is successfully applied to elliptic interface problems [11, 12] and linear parabolic problems [13–16] and further developed for other applications, such as biharmonic problems [17–22], Cahn–Hilliard problems [23], and Stokes problems [24–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. The 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].This paper is organized as follows. InSection 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 L2 norm and the vector-valued 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.


## 2. 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 processq and scalar random process u. Then, 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.


### 2.1. Notations and Variational Formulation

LetH=L2D be the standard square integral space with norm ⋅ and inner product ⋅,⋅ and Hm be the Sobolev space with norm ⋅m for m≥1. Furthermore, denote by H0m the subspace of Hm whose element vanishes at the boundary. Additionally, we define a space(8)Hdiv;D=vv∈L2D2,∇⋅v∈L2D,with a norm(9)vHdiv;D=v2+∇⋅v21/2.Next, we introduce a generalH-valued random variable space as(10)L2Ω;H=v:Ω⟶H∫Ωv⋅,ωH2dℙω<∞,with norm ⋅L2Ω;H2=∫Ω⋅H2dℙω.Now, we can define the variational equations of the stochastic parabolic problems (4)–(7) as follows: seeking ut,qt∈L2Ω;H01×L2Ω;Hdiv;D, such that(11)q,v−∇⋅v,u=0,∀v∈Hdiv;D,t∈0,T,a.s.,(12)ut,χ+∇⋅q,χ=f,χ,∀χ∈L2D,t∈0,T,a.s.,and satisfying the initial condition (6).


### 2.2. Weak Galerkin Approximations

There are many numerical methods to solve equations (1)–(3), such as stochastic collocation methods [39], discontinuous Galerkin methods [40], and finite element methods [2]. The standard finite element approach is to properly choose two finite element spaces Wh⊆H01 and Vh⊆Hdiv;D, respectively, and then project equations (11) and (12) into the space Wh×Vh. Therefore, we obtain finite dimensional equations whose solutions are the approximations of (11) and (12). The recently developed weak Galerkin finite element method adopts the piecewisely defined discontinuous functions to construct the finite element spaces Wh and Vh, 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.LetDh be a partition of the domain D satisfying the shape regularity requirements A1–A4 in [31], hK be the partition diameter of the element K∈Dh, and h=maxKhK. For any element K∈Dh, let K0 and ∂K denote its interior and boundary, respectively. Denote by ℰh the set of all edges in Dh, and let ℰh0=ℰh\∂D be the set of all interior edges.Denote byne a fixed unit normal vector on ℰh satisfying when e∈ℰh0, ne is normal to e, and when e∈ℰh∩∂D, ne is outward normal to e.For each elementK∈Dh, let PjK denote the set of polynomials defined on K with degree no more than j≥0. Similarly, we can define a polynomial space Pje, for each e∈ℰh.Now, we define three weak Galerkin finite element spaces by(13)Wh=χh∈L2DχhK∈Pj+1K,Wh0=χh∈Whχhvanishes on∂D,Vh=vh=v0,vbv0K∈PjK2,vbe=vbne,vb∈Pje,∀K∈Dh,∀e∈ℰh.Remark 1.
The functions inVh are totally discontinuous. Such functions may contain independent information on any K∈Dh and any e∈ℰh.
Next, we extend the divergence operator∇⋅:Hdiv;D⟶H to the space pair Vh,Wh. Define a discrete weak divergence operator ∇w⋅ from Vh to Wh: for any vh=v0,vb∈Vh, ∇w⋅vh∈Wh is defined as a functional on Wh which satisfies(14)∇w⋅vh,φK=−v0,∇φK+vb⋅n,φ∂K,∀φ∈Wh,∀K∈Dh,where n is a unit outward normal vector to ∂K.Remark 2.
LetK1 and K2 be two adjacent elements in Dh sharing a common interior edge e. n1 and n2 are unit outward normal vectors of K1 and K2 to e, respectively. Let ne be the pre-fixed unit normal vector to e. As can be seen from Figure 1, we shall obtain ne=n1 and ne=−n2, or vice versa.
In order to simplify the notations, we introduce some bilinear forms: forηh=η0,ηb∈Vh, ζh=ζ0,ζb∈Vh and ξh∈Wh, such that(15)aηh,ζh=∑K∈Dhη0K,ζ0K,sηh,ζh=∑K∈DhhKη0−ηb⋅n,ζ0−ζb⋅n∂K,asηh,ζh=aηh,ζh+sηh,ζh,bηh,ξh=∇w⋅ηh,ξh.
We also define two norms(16)ζh≔asζh,ζh1/2,for anyζh∈Vh,(17)ζhL2Ω;Vh≔Eζh21/2,for anyζh∈L2Ω;Vh,where E denotes the expectation with respect to the probability measure ℙ. From [30, 31], we know that the bilinear form as⋅,⋅ is bounded, symmetric, and positive definite in Vh, and the bilinear form b⋅,⋅ is bounded in Vh×Wh and they also satisfy the standard inf-sup condition.
Now, we can formulate the semidiscrete weak Galerkin mixed finite element scheme. For anyt∈0,T, finding qh=q0,qb∈L2Ω;Vh and uh∈L2Ω;Wh0 such that(18)asqh,vh−bvh,uh=0,∀vh=v0,vb∈Vh,a.s.,(19)uh,t,χh+bqh,χh=f,χh,∀χh∈Wh,a.s.,and satisfying the initial condition(20)uh⋅,0,ω=Phu0⋅,ω,inD×Ω,where Ph is an L2 projection operator from H01 onto Wh0. 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 tn=nκ for n=0,1,…,N. By qn and un, we denote the approximation of qtn and utn, respectively. We use the backward Euler method to approximate the time derivative in (19); then, the fully discrete scheme reads seeking qn,un∈L2Ω;Vh×L2Ω;Wh0 such that(21)asqn,vh−bvh,un=0,∀vh=v0,vb∈Vh,a.s.,(22)∂¯un,χh+bqn,χh=ftn,χh,∀χh∈Wh,a.s.,and satisfying the initial condition(23)u0⋅,0,w=Phu0⋅,w,inD×Ω,where ∂¯un=un−un−1/κ.Figure 1
Relationship betweenne and n.


## 2.1. Notations and Variational Formulation

LetH=L2D be the standard square integral space with norm ⋅ and inner product ⋅,⋅ and Hm be the Sobolev space with norm ⋅m for m≥1. Furthermore, denote by H0m the subspace of Hm whose element vanishes at the boundary. Additionally, we define a space(8)Hdiv;D=vv∈L2D2,∇⋅v∈L2D,with a norm(9)vHdiv;D=v2+∇⋅v21/2.Next, we introduce a generalH-valued random variable space as(10)L2Ω;H=v:Ω⟶H∫Ωv⋅,ωH2dℙω<∞,with norm ⋅L2Ω;H2=∫Ω⋅H2dℙω.Now, we can define the variational equations of the stochastic parabolic problems (4)–(7) as follows: seeking ut,qt∈L2Ω;H01×L2Ω;Hdiv;D, such that(11)q,v−∇⋅v,u=0,∀v∈Hdiv;D,t∈0,T,a.s.,(12)ut,χ+∇⋅q,χ=f,χ,∀χ∈L2D,t∈0,T,a.s.,and satisfying the initial condition (6).


## 2.2. Weak Galerkin Approximations

There are many numerical methods to solve equations (1)–(3), such as stochastic collocation methods [39], discontinuous Galerkin methods [40], and finite element methods [2]. The standard finite element approach is to properly choose two finite element spaces Wh⊆H01 and Vh⊆Hdiv;D, respectively, and then project equations (11) and (12) into the space Wh×Vh. Therefore, we obtain finite dimensional equations whose solutions are the approximations of (11) and (12). The recently developed weak Galerkin finite element method adopts the piecewisely defined discontinuous functions to construct the finite element spaces Wh and Vh, 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.LetDh be a partition of the domain D satisfying the shape regularity requirements A1–A4 in [31], hK be the partition diameter of the element K∈Dh, and h=maxKhK. For any element K∈Dh, let K0 and ∂K denote its interior and boundary, respectively. Denote by ℰh the set of all edges in Dh, and let ℰh0=ℰh\∂D be the set of all interior edges.Denote byne a fixed unit normal vector on ℰh satisfying when e∈ℰh0, ne is normal to e, and when e∈ℰh∩∂D, ne is outward normal to e.For each elementK∈Dh, let PjK denote the set of polynomials defined on K with degree no more than j≥0. Similarly, we can define a polynomial space Pje, for each e∈ℰh.Now, we define three weak Galerkin finite element spaces by(13)Wh=χh∈L2DχhK∈Pj+1K,Wh0=χh∈Whχhvanishes on∂D,Vh=vh=v0,vbv0K∈PjK2,vbe=vbne,vb∈Pje,∀K∈Dh,∀e∈ℰh.Remark 1.
The functions inVh are totally discontinuous. Such functions may contain independent information on any K∈Dh and any e∈ℰh.
Next, we extend the divergence operator∇⋅:Hdiv;D⟶H to the space pair Vh,Wh. Define a discrete weak divergence operator ∇w⋅ from Vh to Wh: for any vh=v0,vb∈Vh, ∇w⋅vh∈Wh is defined as a functional on Wh which satisfies(14)∇w⋅vh,φK=−v0,∇φK+vb⋅n,φ∂K,∀φ∈Wh,∀K∈Dh,where n is a unit outward normal vector to ∂K.Remark 2.
LetK1 and K2 be two adjacent elements in Dh sharing a common interior edge e. n1 and n2 are unit outward normal vectors of K1 and K2 to e, respectively. Let ne be the pre-fixed unit normal vector to e. As can be seen from Figure 1, we shall obtain ne=n1 and ne=−n2, or vice versa.
In order to simplify the notations, we introduce some bilinear forms: forηh=η0,ηb∈Vh, ζh=ζ0,ζb∈Vh and ξh∈Wh, such that(15)aηh,ζh=∑K∈Dhη0K,ζ0K,sηh,ζh=∑K∈DhhKη0−ηb⋅n,ζ0−ζb⋅n∂K,asηh,ζh=aηh,ζh+sηh,ζh,bηh,ξh=∇w⋅ηh,ξh.
We also define two norms(16)ζh≔asζh,ζh1/2,for anyζh∈Vh,(17)ζhL2Ω;Vh≔Eζh21/2,for anyζh∈L2Ω;Vh,where E denotes the expectation with respect to the probability measure ℙ. From [30, 31], we know that the bilinear form as⋅,⋅ is bounded, symmetric, and positive definite in Vh, and the bilinear form b⋅,⋅ is bounded in Vh×Wh and they also satisfy the standard inf-sup condition.
Now, we can formulate the semidiscrete weak Galerkin mixed finite element scheme. For anyt∈0,T, finding qh=q0,qb∈L2Ω;Vh and uh∈L2Ω;Wh0 such that(18)asqh,vh−bvh,uh=0,∀vh=v0,vb∈Vh,a.s.,(19)uh,t,χh+bqh,χh=f,χh,∀χh∈Wh,a.s.,and satisfying the initial condition(20)uh⋅,0,ω=Phu0⋅,ω,inD×Ω,where Ph is an L2 projection operator from H01 onto Wh0. 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 tn=nκ for n=0,1,…,N. By qn and un, we denote the approximation of qtn and utn, respectively. We use the backward Euler method to approximate the time derivative in (19); then, the fully discrete scheme reads seeking qn,un∈L2Ω;Vh×L2Ω;Wh0 such that(21)asqn,vh−bvh,un=0,∀vh=v0,vb∈Vh,a.s.,(22)∂¯un,χh+bqn,χh=ftn,χh,∀χh∈Wh,a.s.,and satisfying the initial condition(23)u0⋅,0,w=Phu0⋅,w,inD×Ω,where ∂¯un=un−un−1/κ.Figure 1
Relationship betweenne and n.


## 3. 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 inL2 norm ⋅L2Ω;H. We also analyze the convergence order for the vector-valued solution in the norm ⋅L2Ω;Vh. Finally, we study the statistical variance-type convergence estimates for both the scalar and vector-valued solutions, respectively.


### 3.1. Stability Analysis

In this section, we shall derive stability estimates in both⋅L2Ω;H and ⋅L2Ω;Vh norms for the semidiscrete and fully discrete WGMFE schemes. These 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:(24)uhtL2Ω;H2≤Cuh0L2Ω;H2+C∫0tf2dτ,(25)qhtL2Ω;Vh2≤qh0L2Ω;Vh2+∫0tf2dτ,where C>0 is a constant independent of the mesh size h.Proof.
In equations (18) and (19), we take vh=qh∈L2Ω;Vh and χh=uh∈L2Ω;Wh0, then(26)qhω2−bqh,uh=0,(27)uh,tω,uhω+bqh,uh=f,uhω.
Adding (26) to (27), we have(28)qhω2+12ddtuhω2=f,uhω.
Because of Cauchy–Schwarz inequality andqhω2≥0, we obtain(29)12ddtuhω2=uhωddtuhω≤fuhω.
Hence,(30)uht,ω≤uh0,ω+∫0tfdτ.
Squaring both sides of (30) and using Cauchy–Schwarz inequality,(31)uht,ω2≤Cuh0,ω2+C∫0tf2dτ,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,(32)asqh,t,vh−bvh,uh,t=0.
Takingvh=qh and χh=uh,t in (32) and (19) and summing up, we have(33)12ddtqhω2+uh,tω2=f,uh,tω≤12f2+12uh,tω2.
Noticinguh,tω2≥0, we get(34)ddtqhω2≤f2.
Thus,(35)qht,ω2≤qh0,ω2+∫0tf2dτ.
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(36)unL2Ω;H2≤Cu0L2Ω;H2+Cκ2∑i=1nfti2,(37)qnL2Ω;Vh2≤q0L2Ω;Vh2+κ∑i=1nfti2,where C>0 is a constant independent of the mesh size.Proof.
Takevh=qn∈L2Ω;Vh and χh=un∈L2Ω;Wh0 in (21) and (22) and sum up these two equations, we have(38)qnω2+∂¯unω,unω=ftn,unω.
Making use of Cauchy–Schwarz inequality andqnω2≥0, we obtain(39)unω≤un−1ω+κftn,and by recurrence, we obtain(40)unω≤u0ω+κ∑i=1nfti.
Squaring and taking expectation on both sides of (40), we derive (36).
Using the backward Euler method to discretize the time derivative in (32),(41)as∂¯qn,vh−bvh,∂¯un=0.
Takingvh=qn and χh=∂¯un in (41) and (22) and summing up,(42)∂¯unω2+as∂¯qn,qn=ftn,∂¯unω.
Since(43)as∂¯qn,qn=12∂¯asqn,qn+κ2as∂¯qn,∂¯qn,we find(44)∂¯unω2+12∂¯asqn,qn+κ2as∂¯qn,∂¯qn=ftn,∂¯unω.
By virtue of Cauchy–Schwarz inequality and the iteration method, we obtain(45)qnω2≤q0ω2+κ∑i=1nfti2.
Taking expectation on both sides of (45), and making use of the definitions (16) and (17), the proof is completed.


### 3.2. 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 solutionqh. 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 spacesVh and Wh by using local L2 projections. For each element K∈Dh, denote by Q0 the projection operator from L2K2 onto PjK2. Similarly, for each edge e∈ℰh, let Qb be the projection operator from L2e onto Pje. For any v∈Hdiv;D, denote by Qhv=Q0v,Qbv⋅nene the L2 projection of v in Vh. Furthermore, denote by Lh the L2 projection operator from L2K onto Pj+1K for each K∈Dh.


#### 3.2.1. Semidiscrete Error Analysis

For simplicity, we introduce two new notations:(46)θh=θ0,θb≔Qhq−qh,mh≔Lhu−uh.In [41], WGMFEM is applied to study the following deterministic problem: seeking a flux function q and a scalar function u satisfying(47)q+∇u=0,ut+∇⋅q=f,inD,t∈0,T,and the initial-boundary condition(48)u=−g,on∂D,t∈0,T,(49)u⋅,0=u0,inD.According to the results of [41], we have the following useful estimates.Lemma 1 (see [41], Theorem 3.9 and Theorem 3.10).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qh,uh be the semidiscrete solution of WGMFE scheme of (47)–(49). Then, we have(50)mht≤mh0+Chj+2u0j+2+q0j+1+utj+2+qtj+1+∫0tutj+2+qtj+1dτ,(51)θht2≤4θh02+Ch2j+1u0j+2+q0j+12+utj+2+qtj+12+Ch2j+2∫0tutj+2+qtj+12dτ.Now, we shall present convergence analysis for semidiscrete weak Galerkin methods in the next two theorems.Theorem 3.
For anyt∈0,T and j≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qht,uht∈L2Ω;Vh×L2Ω;Wh0 be the solution of (18) and (19). Assuming that f∈L20,T;Hj+2, then(52)ut−uhtL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,(53)∫DVarut−Varuhtdxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ,with Varu denoting the variance of the random variable u.Proof.
Formh∈L2Ω;Wh0, taking account of (50), we obtain, for almost ω,(54)mht,ω≤mh0,ω+Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
By virtue of the definition ofmht,ω, the error estimate of L2 projection, and (54), we have(55)ut,ω−uht,ω≤ut,ω−Lhut,ω+mht,ω≤Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
Squaring both sides of (55) and utilizing Cauchy–Schwarz inequality, it follows that(56)ut,ω−uht,ω2≤Ch2j+2u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+∫0tutωj+22+qtωj+12dτ.
Take expectation on both sides of (56), which leads to (52).
By virtue of Cauchy–Schwarz inequality, we supply(57)Covu−uh,uh≤Varu−uhVaruh,where Cov represents the covariance in regard to the probability measure ℙ. Combining the fact(58)∫DVaruhdxdy≤∫D∫Ωuh2dℙωdxdy=uhL2Ω;H2and Cauchy–Schwarz inequality, we provide(59)2∫DCovu−uh,uhdxdy≤u−uhL2Ω;H2+uhL2Ω;H2.
Additionally,(60)Varu−Varuh≤Varu+Varuh=Varu−uh+uh+Varuh≤Varu−uh+2Covu−uh,uh+3Varuh.
Because of (58)–(60), we find(61)∫DVaru−Varuhdxdy≤Cu−uhL2Ω;H2+CuhL2Ω;H2.
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 operatorTr⋅: for a square matrix A, TrA is defined as the sum of the elements on the main diagonal of A.Theorem 4.
Under the assumption of Theorem3, we have(62)qt−qhtL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,∫DTrVarqt−TrVarq0tdxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ.Proof.
Forθh∈L2Ω;Vh, considering (51), we obtain, for almost ω,(63)θht,ω2≤4θh0,ω2+Ch2j+1u0,ωj+2+q0,ωj+12+ut,ωj+2+qt,ωj+12+Ch2j+2∫0tutωj+2+qtωj+12dτ.
Because of the definition ofθht,ω, Cauchy–Schwarz inequality, the error estimate of L2 projection, and (63), it follows that(64)qt,ω−qht,ω2≤Cqt,ω−Qhqt,ω2+Cθht,ω2≤Ch2j+1u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+Ch2j+2∫0tutωj+22+qtωj+12dτ.
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∈L2Ω;Hj+12 and qh=q0,qb∈Vh; assuming that q=q1,q2T and q0=q01,q02T, then(65)Varq−q0=Eq−q0−Eq−q0q−q0−Eq−q0T.Thus,(66)TrVarq−q0=Varq1−q01+Varq2−q02.
Utilizing (58) and the definition of ⋅L2Ω;Vh, we arrive at(67)∫DTrVarq−q0dxdy=∫DVarq1−q01dxdy+∫DVarq2−q02dxdy≤q1−q01L2Ω;H2+q2−q02L2Ω;H2≤q−qhL2Ω;Vh2.
Making use of (57), we get(68)TrCovq−q0,q0≤Covq1−q01,q01+Covq2−q02,q02≤Varq1−q01Varq01+Varq2−q02Varq02.
Together with Cauchy–Schwarz inequality, (66), and (67), we obtain(69)2∫DTrCovq−q0,q0dxdy≤2∫DVarq1−q01Varq01dxdy+2∫DVarq2−q02Varq02dxdy≤∫DVarq1−q01dxdy+∫DVarq01dxdy+∫DVarq2−q02dxdy+∫DVarq02dxdy=∫DTrVarq−q0dxdy+∫DTrVarq0dxdy≤q−qhL2Ω;Vh2+qhL2Ω;Vh2.
Taking account of the fact(70)TrVarq−TrVarq0≤CTrVarq−q0+TrCovq−q0,q0+TrVarq0and applying (67) and (69), we present(71)∫DTrVarq−TrVarq0dxdy≤Cq−qhL2Ω;Vh2+CqhL2Ω;Vh2.
Coupled with the first inequality of (62), (25), the definition of ⋅L2Ω;Vh, and the regularity for f, the proof is completed.Remark 3.
In view of (65) in Theorem 4, for two-dimensional vector-valued random process q=q1,q2T and q0=q01,q02T, it can be seen that the results of Varq and Varq0 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 q0. This is the reason why we use trace operator here.


#### 3.2.2. Fully Discrete Error Analysis

Define(72)θn=θ0n,θbn≔Qhqtn−qn,mn≔Lhutn−un.Similar to the semidiscrete case, we introduce the fully discrete error estimates of (47)–(49) in [41] as follows.Lemma 2 (see [41], Theorem 3.11 and Theorem 3.12).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qn,un be the fully discrete solution of WGMFE scheme of (47)–(49). Assuming that utt∈L20,T;H, then(73)mn≤m0+Chj+2u0j+2+q0j+1+utnj+2+qtnj+1+∫0tnutj+2+qtj+1dτ+κ∫0tnuttdτθn2≤4θ02+Ch2j+1u0j+2+q0j+12+utnj+2+qtnj+12+Ch2j+2∫0tnutj+2+qtj+12dτ+κ2∫0tnutt2dτ.In a similar way, we can prove the following two theorems for the fully discrete weak Galerkin methods.Theorem 5.
Forj≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qn,un∈L2Ω;Vh×L2Ω;Wh0 be the solution of (21) and (22). Supposing that utt∈L20,T;L2Ω;H, then(74)utn−unL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ∫DVarutn−Varundxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ+∑i=1nfti2.Theorem 6.
Under the assumption of Theorem5, we have(75)qtn−qnL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ2∫0tnuttL2Ω;H2dτ∫DTrVarqtn−TrVarq0ndxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ∑i=1nfti2+κ2∫0tnuttL2Ω;H2dτ.


## 3.1. Stability Analysis

In this section, we shall derive stability estimates in both⋅L2Ω;H and ⋅L2Ω;Vh norms for the semidiscrete and fully discrete WGMFE schemes. These 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:(24)uhtL2Ω;H2≤Cuh0L2Ω;H2+C∫0tf2dτ,(25)qhtL2Ω;Vh2≤qh0L2Ω;Vh2+∫0tf2dτ,where C>0 is a constant independent of the mesh size h.Proof.
In equations (18) and (19), we take vh=qh∈L2Ω;Vh and χh=uh∈L2Ω;Wh0, then(26)qhω2−bqh,uh=0,(27)uh,tω,uhω+bqh,uh=f,uhω.
Adding (26) to (27), we have(28)qhω2+12ddtuhω2=f,uhω.
Because of Cauchy–Schwarz inequality andqhω2≥0, we obtain(29)12ddtuhω2=uhωddtuhω≤fuhω.
Hence,(30)uht,ω≤uh0,ω+∫0tfdτ.
Squaring both sides of (30) and using Cauchy–Schwarz inequality,(31)uht,ω2≤Cuh0,ω2+C∫0tf2dτ,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,(32)asqh,t,vh−bvh,uh,t=0.
Takingvh=qh and χh=uh,t in (32) and (19) and summing up, we have(33)12ddtqhω2+uh,tω2=f,uh,tω≤12f2+12uh,tω2.
Noticinguh,tω2≥0, we get(34)ddtqhω2≤f2.
Thus,(35)qht,ω2≤qh0,ω2+∫0tf2dτ.
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(36)unL2Ω;H2≤Cu0L2Ω;H2+Cκ2∑i=1nfti2,(37)qnL2Ω;Vh2≤q0L2Ω;Vh2+κ∑i=1nfti2,where C>0 is a constant independent of the mesh size.Proof.
Takevh=qn∈L2Ω;Vh and χh=un∈L2Ω;Wh0 in (21) and (22) and sum up these two equations, we have(38)qnω2+∂¯unω,unω=ftn,unω.
Making use of Cauchy–Schwarz inequality andqnω2≥0, we obtain(39)unω≤un−1ω+κftn,and by recurrence, we obtain(40)unω≤u0ω+κ∑i=1nfti.
Squaring and taking expectation on both sides of (40), we derive (36).
Using the backward Euler method to discretize the time derivative in (32),(41)as∂¯qn,vh−bvh,∂¯un=0.
Takingvh=qn and χh=∂¯un in (41) and (22) and summing up,(42)∂¯unω2+as∂¯qn,qn=ftn,∂¯unω.
Since(43)as∂¯qn,qn=12∂¯asqn,qn+κ2as∂¯qn,∂¯qn,we find(44)∂¯unω2+12∂¯asqn,qn+κ2as∂¯qn,∂¯qn=ftn,∂¯unω.
By virtue of Cauchy–Schwarz inequality and the iteration method, we obtain(45)qnω2≤q0ω2+κ∑i=1nfti2.
Taking expectation on both sides of (45), and making use of the definitions (16) and (17), the proof is completed.


## 3.2. 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 solutionqh. 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 spacesVh and Wh by using local L2 projections. For each element K∈Dh, denote by Q0 the projection operator from L2K2 onto PjK2. Similarly, for each edge e∈ℰh, let Qb be the projection operator from L2e onto Pje. For any v∈Hdiv;D, denote by Qhv=Q0v,Qbv⋅nene the L2 projection of v in Vh. Furthermore, denote by Lh the L2 projection operator from L2K onto Pj+1K for each K∈Dh.


### 3.2.1. Semidiscrete Error Analysis

For simplicity, we introduce two new notations:(46)θh=θ0,θb≔Qhq−qh,mh≔Lhu−uh.In [41], WGMFEM is applied to study the following deterministic problem: seeking a flux function q and a scalar function u satisfying(47)q+∇u=0,ut+∇⋅q=f,inD,t∈0,T,and the initial-boundary condition(48)u=−g,on∂D,t∈0,T,(49)u⋅,0=u0,inD.According to the results of [41], we have the following useful estimates.Lemma 1 (see [41], Theorem 3.9 and Theorem 3.10).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qh,uh be the semidiscrete solution of WGMFE scheme of (47)–(49). Then, we have(50)mht≤mh0+Chj+2u0j+2+q0j+1+utj+2+qtj+1+∫0tutj+2+qtj+1dτ,(51)θht2≤4θh02+Ch2j+1u0j+2+q0j+12+utj+2+qtj+12+Ch2j+2∫0tutj+2+qtj+12dτ.Now, we shall present convergence analysis for semidiscrete weak Galerkin methods in the next two theorems.Theorem 3.
For anyt∈0,T and j≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qht,uht∈L2Ω;Vh×L2Ω;Wh0 be the solution of (18) and (19). Assuming that f∈L20,T;Hj+2, then(52)ut−uhtL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,(53)∫DVarut−Varuhtdxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ,with Varu denoting the variance of the random variable u.Proof.
Formh∈L2Ω;Wh0, taking account of (50), we obtain, for almost ω,(54)mht,ω≤mh0,ω+Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
By virtue of the definition ofmht,ω, the error estimate of L2 projection, and (54), we have(55)ut,ω−uht,ω≤ut,ω−Lhut,ω+mht,ω≤Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
Squaring both sides of (55) and utilizing Cauchy–Schwarz inequality, it follows that(56)ut,ω−uht,ω2≤Ch2j+2u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+∫0tutωj+22+qtωj+12dτ.
Take expectation on both sides of (56), which leads to (52).
By virtue of Cauchy–Schwarz inequality, we supply(57)Covu−uh,uh≤Varu−uhVaruh,where Cov represents the covariance in regard to the probability measure ℙ. Combining the fact(58)∫DVaruhdxdy≤∫D∫Ωuh2dℙωdxdy=uhL2Ω;H2and Cauchy–Schwarz inequality, we provide(59)2∫DCovu−uh,uhdxdy≤u−uhL2Ω;H2+uhL2Ω;H2.
Additionally,(60)Varu−Varuh≤Varu+Varuh=Varu−uh+uh+Varuh≤Varu−uh+2Covu−uh,uh+3Varuh.
Because of (58)–(60), we find(61)∫DVaru−Varuhdxdy≤Cu−uhL2Ω;H2+CuhL2Ω;H2.
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 operatorTr⋅: for a square matrix A, TrA is defined as the sum of the elements on the main diagonal of A.Theorem 4.
Under the assumption of Theorem3, we have(62)qt−qhtL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,∫DTrVarqt−TrVarq0tdxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ.Proof.
Forθh∈L2Ω;Vh, considering (51), we obtain, for almost ω,(63)θht,ω2≤4θh0,ω2+Ch2j+1u0,ωj+2+q0,ωj+12+ut,ωj+2+qt,ωj+12+Ch2j+2∫0tutωj+2+qtωj+12dτ.
Because of the definition ofθht,ω, Cauchy–Schwarz inequality, the error estimate of L2 projection, and (63), it follows that(64)qt,ω−qht,ω2≤Cqt,ω−Qhqt,ω2+Cθht,ω2≤Ch2j+1u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+Ch2j+2∫0tutωj+22+qtωj+12dτ.
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∈L2Ω;Hj+12 and qh=q0,qb∈Vh; assuming that q=q1,q2T and q0=q01,q02T, then(65)Varq−q0=Eq−q0−Eq−q0q−q0−Eq−q0T.Thus,(66)TrVarq−q0=Varq1−q01+Varq2−q02.
Utilizing (58) and the definition of ⋅L2Ω;Vh, we arrive at(67)∫DTrVarq−q0dxdy=∫DVarq1−q01dxdy+∫DVarq2−q02dxdy≤q1−q01L2Ω;H2+q2−q02L2Ω;H2≤q−qhL2Ω;Vh2.
Making use of (57), we get(68)TrCovq−q0,q0≤Covq1−q01,q01+Covq2−q02,q02≤Varq1−q01Varq01+Varq2−q02Varq02.
Together with Cauchy–Schwarz inequality, (66), and (67), we obtain(69)2∫DTrCovq−q0,q0dxdy≤2∫DVarq1−q01Varq01dxdy+2∫DVarq2−q02Varq02dxdy≤∫DVarq1−q01dxdy+∫DVarq01dxdy+∫DVarq2−q02dxdy+∫DVarq02dxdy=∫DTrVarq−q0dxdy+∫DTrVarq0dxdy≤q−qhL2Ω;Vh2+qhL2Ω;Vh2.
Taking account of the fact(70)TrVarq−TrVarq0≤CTrVarq−q0+TrCovq−q0,q0+TrVarq0and applying (67) and (69), we present(71)∫DTrVarq−TrVarq0dxdy≤Cq−qhL2Ω;Vh2+CqhL2Ω;Vh2.
Coupled with the first inequality of (62), (25), the definition of ⋅L2Ω;Vh, and the regularity for f, the proof is completed.Remark 3.
In view of (65) in Theorem 4, for two-dimensional vector-valued random process q=q1,q2T and q0=q01,q02T, it can be seen that the results of Varq and Varq0 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 q0. This is the reason why we use trace operator here.


### 3.2.2. Fully Discrete Error Analysis

Define(72)θn=θ0n,θbn≔Qhqtn−qn,mn≔Lhutn−un.Similar to the semidiscrete case, we introduce the fully discrete error estimates of (47)–(49) in [41] as follows.Lemma 2 (see [41], Theorem 3.11 and Theorem 3.12).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qn,un be the fully discrete solution of WGMFE scheme of (47)–(49). Assuming that utt∈L20,T;H, then(73)mn≤m0+Chj+2u0j+2+q0j+1+utnj+2+qtnj+1+∫0tnutj+2+qtj+1dτ+κ∫0tnuttdτθn2≤4θ02+Ch2j+1u0j+2+q0j+12+utnj+2+qtnj+12+Ch2j+2∫0tnutj+2+qtj+12dτ+κ2∫0tnutt2dτ.In a similar way, we can prove the following two theorems for the fully discrete weak Galerkin methods.Theorem 5.
Forj≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qn,un∈L2Ω;Vh×L2Ω;Wh0 be the solution of (21) and (22). Supposing that utt∈L20,T;L2Ω;H, then(74)utn−unL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ∫DVarutn−Varundxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ+∑i=1nfti2.Theorem 6.
Under the assumption of Theorem5, we have(75)qtn−qnL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ2∫0tnuttL2Ω;H2dτ∫DTrVarqtn−TrVarq0ndxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ∑i=1nfti2+κ2∫0tnuttL2Ω;H2dτ.


## 3.2.1. Semidiscrete Error Analysis

For simplicity, we introduce two new notations:(46)θh=θ0,θb≔Qhq−qh,mh≔Lhu−uh.In [41], WGMFEM is applied to study the following deterministic problem: seeking a flux function q and a scalar function u satisfying(47)q+∇u=0,ut+∇⋅q=f,inD,t∈0,T,and the initial-boundary condition(48)u=−g,on∂D,t∈0,T,(49)u⋅,0=u0,inD.According to the results of [41], we have the following useful estimates.Lemma 1 (see [41], Theorem 3.9 and Theorem 3.10).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qh,uh be the semidiscrete solution of WGMFE scheme of (47)–(49). Then, we have(50)mht≤mh0+Chj+2u0j+2+q0j+1+utj+2+qtj+1+∫0tutj+2+qtj+1dτ,(51)θht2≤4θh02+Ch2j+1u0j+2+q0j+12+utj+2+qtj+12+Ch2j+2∫0tutj+2+qtj+12dτ.Now, we shall present convergence analysis for semidiscrete weak Galerkin methods in the next two theorems.Theorem 3.
For anyt∈0,T and j≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qht,uht∈L2Ω;Vh×L2Ω;Wh0 be the solution of (18) and (19). Assuming that f∈L20,T;Hj+2, then(52)ut−uhtL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,(53)∫DVarut−Varuhtdxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ,with Varu denoting the variance of the random variable u.Proof.
Formh∈L2Ω;Wh0, taking account of (50), we obtain, for almost ω,(54)mht,ω≤mh0,ω+Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
By virtue of the definition ofmht,ω, the error estimate of L2 projection, and (54), we have(55)ut,ω−uht,ω≤ut,ω−Lhut,ω+mht,ω≤Chj+2u0,ωj+2+q0,ωj+1+ut,ωj+2+qt,ωj+1+∫0tutωj+2+qtωj+1dτ.
Squaring both sides of (55) and utilizing Cauchy–Schwarz inequality, it follows that(56)ut,ω−uht,ω2≤Ch2j+2u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+∫0tutωj+22+qtωj+12dτ.
Take expectation on both sides of (56), which leads to (52).
By virtue of Cauchy–Schwarz inequality, we supply(57)Covu−uh,uh≤Varu−uhVaruh,where Cov represents the covariance in regard to the probability measure ℙ. Combining the fact(58)∫DVaruhdxdy≤∫D∫Ωuh2dℙωdxdy=uhL2Ω;H2and Cauchy–Schwarz inequality, we provide(59)2∫DCovu−uh,uhdxdy≤u−uhL2Ω;H2+uhL2Ω;H2.
Additionally,(60)Varu−Varuh≤Varu+Varuh=Varu−uh+uh+Varuh≤Varu−uh+2Covu−uh,uh+3Varuh.
Because of (58)–(60), we find(61)∫DVaru−Varuhdxdy≤Cu−uhL2Ω;H2+CuhL2Ω;H2.
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 operatorTr⋅: for a square matrix A, TrA is defined as the sum of the elements on the main diagonal of A.Theorem 4.
Under the assumption of Theorem3, we have(62)qt−qhtL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12dτ,∫DTrVarqt−TrVarq0tdxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utL2Ω;Hj+22+qtL2Ω;Hj+12+Ch2j+2∫0tutL2Ω;Hj+22+qtL2Ω;Hj+12+fj+22dτ.Proof.
Forθh∈L2Ω;Vh, considering (51), we obtain, for almost ω,(63)θht,ω2≤4θh0,ω2+Ch2j+1u0,ωj+2+q0,ωj+12+ut,ωj+2+qt,ωj+12+Ch2j+2∫0tutωj+2+qtωj+12dτ.
Because of the definition ofθht,ω, Cauchy–Schwarz inequality, the error estimate of L2 projection, and (63), it follows that(64)qt,ω−qht,ω2≤Cqt,ω−Qhqt,ω2+Cθht,ω2≤Ch2j+1u0,ωj+22+q0,ωj+12+ut,ωj+22+qt,ωj+12+Ch2j+2∫0tutωj+22+qtωj+12dτ.
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∈L2Ω;Hj+12 and qh=q0,qb∈Vh; assuming that q=q1,q2T and q0=q01,q02T, then(65)Varq−q0=Eq−q0−Eq−q0q−q0−Eq−q0T.Thus,(66)TrVarq−q0=Varq1−q01+Varq2−q02.
Utilizing (58) and the definition of ⋅L2Ω;Vh, we arrive at(67)∫DTrVarq−q0dxdy=∫DVarq1−q01dxdy+∫DVarq2−q02dxdy≤q1−q01L2Ω;H2+q2−q02L2Ω;H2≤q−qhL2Ω;Vh2.
Making use of (57), we get(68)TrCovq−q0,q0≤Covq1−q01,q01+Covq2−q02,q02≤Varq1−q01Varq01+Varq2−q02Varq02.
Together with Cauchy–Schwarz inequality, (66), and (67), we obtain(69)2∫DTrCovq−q0,q0dxdy≤2∫DVarq1−q01Varq01dxdy+2∫DVarq2−q02Varq02dxdy≤∫DVarq1−q01dxdy+∫DVarq01dxdy+∫DVarq2−q02dxdy+∫DVarq02dxdy=∫DTrVarq−q0dxdy+∫DTrVarq0dxdy≤q−qhL2Ω;Vh2+qhL2Ω;Vh2.
Taking account of the fact(70)TrVarq−TrVarq0≤CTrVarq−q0+TrCovq−q0,q0+TrVarq0and applying (67) and (69), we present(71)∫DTrVarq−TrVarq0dxdy≤Cq−qhL2Ω;Vh2+CqhL2Ω;Vh2.
Coupled with the first inequality of (62), (25), the definition of ⋅L2Ω;Vh, and the regularity for f, the proof is completed.Remark 3.
In view of (65) in Theorem 4, for two-dimensional vector-valued random process q=q1,q2T and q0=q01,q02T, it can be seen that the results of Varq and Varq0 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 q0. This is the reason why we use trace operator here.


## 3.2.2. Fully Discrete Error Analysis

Define(72)θn=θ0n,θbn≔Qhqtn−qn,mn≔Lhutn−un.Similar to the semidiscrete case, we introduce the fully discrete error estimates of (47)–(49) in [41] as follows.Lemma 2 (see [41], Theorem 3.11 and Theorem 3.12).
Letq,u∈H10,T;Hj+12×H10,T;Hj+2 be the solution of (47)–(49) and qn,un be the fully discrete solution of WGMFE scheme of (47)–(49). Assuming that utt∈L20,T;H, then(73)mn≤m0+Chj+2u0j+2+q0j+1+utnj+2+qtnj+1+∫0tnutj+2+qtj+1dτ+κ∫0tnuttdτθn2≤4θ02+Ch2j+1u0j+2+q0j+12+utnj+2+qtnj+12+Ch2j+2∫0tnutj+2+qtj+12dτ+κ2∫0tnutt2dτ.In a similar way, we can prove the following two theorems for the fully discrete weak Galerkin methods.Theorem 5.
Forj≥0, let q,u∈H10,T;L2Ω;Hj+12×H10,T;L2Ω;H0j+2 be the solution of (11) and (12) and qn,un∈L2Ω;Vh×L2Ω;Wh0 be the solution of (21) and (22). Supposing that utt∈L20,T;L2Ω;H, then(74)utn−unL2Ω;H2≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ∫DVarutn−Varundxdy≤Ch2j+2u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+Cκ2∫0tnuttL2Ω;H2dτ+∑i=1nfti2.Theorem 6.
Under the assumption of Theorem5, we have(75)qtn−qnL2Ω;Vh2≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ2∫0tnuttL2Ω;H2dτ∫DTrVarqtn−TrVarq0ndxdy≤Ch2j+1u0L2Ω;Hj+22+q0L2Ω;Hj+12+utnL2Ω;Hj+22+qtnL2Ω;Hj+12+Ch2j+2∫0tnutL2Ω;Hj+22+qtL2Ω;Hj+12dτ+κ∑i=1nfti2+κ2∫0tnuttL2Ω;H2dτ.


## 4. Numerical Experiments

In this section, we will provide two examples to illustrate our theoretical results. LetD=0,1×0,1. We take two different probability spaces in our calculation. One is a standard normal distribution N0,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:(76)q+∇u=0,forx,t,ω∈D×0,1×Ω,(77)ut+∇⋅q=fx,t,forx,t,ω∈D×0,1×Ω,(78)ux,0,ω=u0x,ω,forx,ω∈D×Ω,(79)ux,t,ω=gx,t,forx,t,ω∈∂D×0,1×Ω.We utilize a uniform triangular meshDh to partition D and conduct numerical experiments with the fully discrete weak Galerkin mixed finite element schemes (21)–(23) by choosing j=0 in Wh0 and Vh. The Monte Carlo sampling method is applied to evaluate the expectation and variance appearing in Theorem 5 and Theorem 6.For simplicity, we firstly introduce several abbreviation notations in the following tables for presenting the numerical results:(80)eun≔utn−un,eqn≔qtn−qn,run≔Varutn−Varun,rqn≔TrVarqtn−TrVarq0n.Example 1.
Consider equations (76)–(79) with(81)f=x2+y2sinxycost−sinxysint,u0=sinxy+ω,and the boundary condition(82)u=0,ony=0,t∈0,1,u=sinycost,onx=1,t∈0,1,u=sinxcost,ony=1,t∈0,1,u=0,onx=0,t∈0,1.
Then, the exact solution is(83)q=−ycosxy+ωcost,−xcosxy+ωcostT,u=sinxy+ωcost.
We carry out the numerical calculation with two different choices of the time stepsizeκ=h and κ=h2. In Tables 1 and 2, we use a normal distribution as the random input and we also present all errors and convergence rates. Similar results are shown in Tables 3 and 4, where we use a uniform distribution as the random input.Table 1
Error estimates of Example 1 with a normal distribution,κ=h, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
20
1.294e–04
1.847e–08
1.028e–02
1.153e–04
1
/
40
5.721e–05
1.178
3.586e–09
2.365
5.136e–03
1.002
5.758e–05
1.002
1
/
60
3.761e–05
1.035
1.541e–09
2.083
3.422e–03
1.001
3.836e–05
1.001
1
/
80
2.828e–05
0.991
8.651e–10
2.007
2.566e–03
1.001
2.877e–05
1.001
1
/
100
2.268e–05
0.989
5.552e–10
1.988
2.053e–03
1.000
2.301e–05
1.000Table 2
Error estimates of Example 1 with a normal distribution,κ=h2, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
4
2.477e–03
6.799e–06
5.272e–02
3.074e–03
1
/
8
6.083e–04
2.026
4.062e–07
4.065
2.636e–02
1.000
7.684e–04
2.000
1
/
16
1.508e–04
2.012
2.484e–08
4.031
1.318e–02
1.001
1.919e–04
2.002
1
/
32
3.765e–05
2.003
1.542e–09
4.010
6.582e–03
1.001
4.787e–05
2.003Table 3
Error estimates of Example 1 with the uniform distribution of0,1, κ=h, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
order
e
q
n
L
2
Ω
;
V
h
order
∫
D
r
q
n
d
x
d
y
Order
1
/
20
1.412e–04
2.209e–08
1.252e–02
1.730e–04
1
/
40
6.263e–05
1.173
4.327e–09
2.352
6.253e–03
1.001
8.641e–05
1.002
1
/
60
4.136e–05
1.023
1.877e–09
2.060
4.167e–03
1.001
5.759e–05
1.001
1
/
80
3.113e–05
0.987
1.060e–09
1.987
3.125e–03
1.000
4.318e–05
1.001
1
/
100
2.500e–05
0.983
6.807e–10
1.984
2.500e–03
1.000
3.455e–05
1.000Table 4
Error estimates of Example 1 with the uniform distribution of0,1, κ=h2, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
4
2.694e–03
8.047e–06
6.346e–02
4.461e–03
1
/
8
6.620e–04
2.025
4.844e–07
4.054
3.174e–02
1.000
1.116e–03
2.000
1
/
16
1.644e–04
2.009
2.983e–08
4.022
1.586e–02
1.001
2.786e–04
2.001
1
/
32
4.105e–05
2.002
1.853e–09
4.009
7.926e–03
1.001
6.953e–05
2.003Example 2.
In equations (76)–(79), let(84)f=4πtcos2πt2+1+π2sinπyexp−x+π2−1sin2πt2+1+π2sinπyexp−x,u0=sinπy+ωexp−x,and the boundary condition (85)u=0,ony=0,t∈0,1,u=sin2πt2+1+π2sinπyexp−1,onx=1,t∈0,1,u=0,ony=1,t∈0,1,u=sin2πt2+1+π2sinπy,onx=0,t∈0,1.
Then, the exact solution is(86)q=sin2πt2+1+π2sinπy+ωexp−x,−πsin2πt2+1+π2cosπy+ωexp−xT,u=sin2πt2+1+π2sinπy+ωexp−x.
In this test, Tables5 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.Table 5
Error estimates of Example 2 with a normal distribution,κ=h, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
20
3.427e–02
1.303e–03
1.604e–01
2.798e–02
1
/
40
1.860e–02
0.882
3.635e–04
1.841
8.390e–02
0.935
1.460-e–02
0.939
1
/
60
1.270e–02
0.941
1.678e–04
1.907
5.685e–02
0.960
9.886e–03
0.961
1
/
80
9.590e–03
0.975
9.538e–05
1.963
4.293e–02
0.976
7.457e–03
0.980
1
/
100
7.695e–03
0.987
6.119e–05
1.989
3.443e–02
0.989
5.977e–03
0.992Table 6
Error estimates of Example 2 with a normal distribution,κ=h2, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
4
3.336e–02
1.235e–03
4.411e–01
2.147e–01
1
/
8
8.940e–03
1.900
8.428e–05
3.873
2.091e–01
1.077
4.818e–02
2.156
1
/
16
2.241e–03
1.996
5.242e–06
4.007
1.020e–01
1.036
1.139e–02
2.081
1
/
32
5.583e–04
2.005
3.247e–07
4.013
5.052e–02
1.013
2.779e–03
2.035Table 7
Error estimates of Example 2 with the uniform distribution of0,1, κ=h, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
20
3.778e–02
1.583e–03
1.967e–01
4.290e–02
1
/
40
2.063e–02
0.873
4.686e–04
1.756
1.033e–01
0.929
2.244e–02
0.934
1
/
60
1.414e–02
0.931
2.186e–04
1.881
7.012e–02
0.956
1.521e–02
0.959
1
/
80
1.074e–02
0.959
1.252e–04
1.937
5.303e–02
0.971
1.149e–02
0.975
1
/
100
8.645e–03
0.970
8.088e–05
1.959
4.259e–02
0.981
9.215e–03
0.988Table 8
Error estimates of Example 2 with the uniform distribution of0,1, κ=h2, n=1/κ, and j=0.
h
e
u
n
L
2
Ω
;
H
Order
∫
D
r
u
n
d
x
d
y
Order
e
q
n
L
2
Ω
;
V
h
Order
∫
D
r
q
n
d
x
d
y
Order
1
/
4
3.681e–02
1.504e–03
5.456e–01
3.303e–01
1
/
8
9.906e–03
1.894
1.082e–04
3.797
2.594e–01
1.073
7.441e–02
2.150
1
/
16
2.493e–03
1.991
6.778e–06
3.996
1.265e–01
1.036
1.763e–02
2.078
1
/
32
6.225e–04
2.002
4.215e–07
4.007
6.273e–02
1.012
4.313e–03
2.031



---

*Source: 8796345-2020-01-28.xml*