A Weak Galerkin Method with Implicit \theta-schemes for Second-Order Parabolic Problems

# A Weak Galerkin Method with Implicit θ-schemes for Second-Order Parabolic Problems

Wenya Qi School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, PR China
###### Abstract

We introduce a new weak Galerkin finite element method whose weak functions on interior neighboring edges are double-valued for parabolic problems. Based on element, a fully discrete approach is formulated with implicit -schemes in time for , which include first-order backward Euler and second-order Crank-Nicolson schemes. Moreover, the optimal convergence rates in the and energy norms are derived. Numerical example is given to verify the theory.

###### keywords:
parabolic problem, weak Galerkin, error estimate, -schemes, backward Euler scheme, Crank-Nicolson scheme
journal: \newdefinition

exampleExample

## 1 Introduction

In this paper, an extension of weak Galerkin finite element method (WG) in Wang2013A () to parabolic problems will be introduced, and referred to over-penalized weak Galerkin finite element method (OPWG). Different from single-valued weak functions on interior edges in WG (Wang2013A (), Mu2015poly (), Mu2013Weakinterface1 ()), double-valued weak functions appeared in liu2018an () has been employed to strengthen flexibility of WG with element. For realizing weak continuity of WG, we naturally deal with jumps on the interior edges by the penalty terms. Importantly, penalized terms on weak functions will be analyzed with sharp penalized parameters explicitly given.

Let be an open and bounded polygonal or polyhedral domain. A linear parabolic model is listed as follows

 ut−∇⋅(A∇u) =f(x,t),          in Ω×(0,¯T], (1.1) u =g(x,t),          on ∂Ω×(0,¯T], (1.2) u(x,0) =φ(x),            in Ω, (1.3)

where the functions , and are known in some specific spaces for well-posedness. The coefficient matrix is symmetric positive, i.e., there exist two positive constants and such that for each

 (Aw,v) ≤β1∥w∥∥v∥, (1.4) (Av,v) ≥α1∥v∥2.

With different approximation spaces for weak gradient operator, WG with element and element were developed for the parabolic equations in Li2013Weak () and Gao2014ON (), respectively, first-order backward Euler full-discrete scheme being investigated. However, there are few publications on second-order fully discrete WG schemes. Based on element Raviart1977A () and -schemes, optimal convergence of the fully discrete OPWG approximations will be analyzed in this paper. Note that we concern about double-valued weak functions and if the jumps go to zero along the interior edges, the usual WG method can be recovered.

The paper is organized as follows. In Sec. , the semi-discrete and full-discrete OPWG schemes are introduced and the latter is unconditionally stable. In Sec. , optimal convergence analysis is presented including error estimates in the and energy norms. Finally, numerical results demonstrate the efficiency and feasibility of the new method.

Throughout this paper, we denote by an arbitrarily small positive constant, the -norm and with the spaces with respect to time where represents Sobolev space (see details in Adams2003Sobolev () or Quarteroni1994 ()). Moreover, we use for a positive constant independent of mesh size and time step .

## 2 OPWG schemes and stability

Let be a partition of domain satisfying shape regularity in Wang2014A (). For each element , is its diameter and is the mesh size of . Denote by the set of interior edges or flat faces, and the edges or flat faces of element . Let be the space of polynomials of degree less than or equal to in variables. The weak Galerkin finite element space for OPWG is defined as

 Vh:= {(v0,vb):v0|T∈Pk(T),T∈Th;vb|e∈Pk(e)×Pk(e),e∈EI; vb|e∈Pk(e),e∈∂Ω,k≥0},

in particular, . For each , we define a unique local weak gradient on each element satisfying

 (∇wv,q)T=−(v0,∇⋅q)T+⟨vb,q⋅n⟩∂T, ∀ q∈RTk(T).

In addition, we define several local projection operators onto space . For each and , let , be the projection operators to and , respectively. Denote by . Meantime, define the projection onto , and then one can obtain

 ∇w(Qhv)=Rh(∇v), ∀v∈H1(T).

Furthermore, we define a div projection for satisfying that and on each element , and (see Wang2013A () and Brezzi1991Mixed ())

 (∇⋅q,v0)T=(∇⋅Πhq,v0)T, ∀ v0∈Pk(T). (2.1)

Then, an approximation property of the projection is given.

###### Lemma 2.1.

Wang2013A () For , it holds

 ∥Πh(A∇u)−A∇w(Qhu)∥≤Chk+1∥u∥k+2. (2.2)

For the sake of achieving the scheme of OPWG for parabolic problem (1.1), it is necessary to define a weak bilinear form in the following equation

 aw(v,χ):=(A∇wv,∇wχ)+J0(v,χ), ∀ v,χ∈Vh,

where the penalty term is well defined as

 J0(v,χ):=∑e∈EI|e|−β0⟨\llbracketvb\rrbracket,\llbracketχb\rrbracket⟩e, β0≥1.

Let be shared by adjacent elements and , then we define the jump on by .

The semi-discrete OPWG scheme for (1.1)-(1.3) is to seek satisfying the boundary condition on and the initial condition such that

 ((u0)t,v0)+aw(uh,vh)=(f,v0),  ∀ vh∈V0h. (2.3)

Now, we define energy norm as for any

 \interleavev\interleave2:=aw(v,v).

The existence and uniqueness of semi-discrete solution of (2.3) are obtained from coercivity and continuity of (see Lemma in liu2018an ()).

Next, we present full-discrete OPWG schemes. The interval is divided into subintervals by time step uniformly, i.e. . With the -schemes applied, the full-discrete OPWG schemes are to seek satisfying the boundary condition on and the initial condition such that

 (¯∂un,v0)+aw(θun+(1−θ)un−1,vh) (2.4) =(θf(tn)+(1−θ)f(tn−1),v0),  ∀ vh∈V0h,

where the parameters vary in and difference quotient . For simplification, we denote . Moreover, when , (2.4) is backward Euler scheme, and Crank-Nicolson (CN) scheme is recovered if .

Let be a small subdomain. The flux in time interval holds

 ∫t+∇tt−∇t∫Kutdxdt+∫t+∇tt−∇t∫∂Kq⋅ndsdt=∫t+∇tt−∇t∫Kfdxdt,

where is the flow rate of heat energy. Multiplying a test function such that in and elsewhere in (2.3), we can obtain that

 ∫t+∇tt−∇t∫K(uh)tdxdt −∫t+∇tt−∇t∫∂KRh(A∇wuh)⋅ndsdt=∫t+∇tt−∇t∫Kfdxdt.

which implies mass conservation, by taking a numerical flux .

### 2.1 Stability of the full-discrete scheme

At first, we will give the following Poincar-type inequality between the norm and the energy norm.

###### Lemma 2.2.

For any , it holds

 ∥v0∥≤C\interleavev\interleave. (2.5)
###### Proof.

Based on Theorem in Chen1998 (), we know that the weak solution of elliptic problem satisfies -regularity, i.e. . Denote by and it is obvious that . With the use of the definitions of (2.1) and discrete weak gradient, trace inequality Wang2014A () and Cauchy-Schwarz inequality, then we can deduce that

 ∥v0∥2 =∑T∈Th(v0,∇⋅q)T=∑T∈Th(v0,∇⋅Πhq)T (2.6) =−∑T∈Th(∇wv,Πhq)T+∑e∈EI⟨\llbracketvb\rrbracket,Πhq⋅ne⟩e ≤∥∇wv∥∥Πhq∥+∑e∈EI∥\llbracketvb\rrbracket∥e∥Πhq∥e ≤∥∇wv∥∥Πhq∥+(∑e∈EI∥\llbracketvb\rrbracket∥e)h−12∥Πhq∥,

where is a unit normal on . From Lemma 2.1, it follows

 ∥Πhq∥ =∥Πh∇Ψ∥≤∥Πh∇Ψ−∇w(QhΨ)∥+∥∇w(QhΨ)∥ (2.7) ≤Ch∥Ψ∥2+∥Rh∇Ψ∥ ≤C∥v0∥.

Combining (2.6) with (2.7) leads to

 ∥v0∥ ≤C(∥∇wv∥+∑e∈EIh−12∥\llbracketvb\rrbracket∥e) ≤C(∥∇wv∥+hβ0(d−1)−12∑e∈EI|e|−β02∥\llbracketvb\rrbracket∥e) ≤C\interleavev\interleave.

###### Theorem 2.3.

Let be the numerical solution of (2.4). Assume i.e. the parabolic problem is homogeneous problem and is bounded in . Then there exists a positive constant such that

 ∥un∥≤∥u0∥+Csupt∈[0,¯T]∥f(t)∥.
###### Proof.

By Cauchy-Schwarz inequality and (2.5), taking

 vh=θun+(1−θ)un−1=(θ−12)(un−un−1)+12(un+un−1)

in (2.4) yields

 12∥un∥2−12∥un−1∥2+(θ−12)∥un−un−1∥2+τ\interleaveθun+(1−θ)un−1\interleave2 (2.8) =τ(θf(tn)+(1−θ)f(tn−1),θun+(1−θ)un−1) ≤τ4ε∥θf(tn)+(1−θ)f(tn−1)∥2+ετ\interleaveθun+(1−θ)un−1\interleave2,

where . Let in (2.8), then it follows

 ∥un∥2≤∥un−1∥2+Cτ∥θf(tn)+(1−θ)f(tn−1)∥2.

Therefore, summing the above inequality from to , and with the boundedness of source function , we obtain that

 ∥un∥2 ≤∥u0∥2+Cτn∑j=1∥θf(tj)+(1−θ)f(tj−1)∥2 ≤∥u0∥2+C¯Tsupt∈[0,¯T]∥f(t)∥2,

and then the conclusion follows. ∎

## 3 Optimal convergence orders

Thanks to elliptic projection, we will establish optimal convergence analysis of the fully discrete OPWG schemes.

For , we define an elliptic projection satisfying the following equation

 aw(Ehv,χ)=−(∇⋅A∇v,χ0),  ∀ χ∈V0h, (3.1)

where is the projection of the trace of on the boundary. Then is the OPWG approximation of the solution of the elliptic problem

 −∇⋅(A∇v) =f∗,  in Ω, v =g,    on ∂Ω.

The following error estimates for the elliptic projection will be used later (see liu2018an ()).

###### Lemma 3.1.

liu2018an () Let , then there exists a positive constant such that

 \interleaveQhu−Ehu\interleave ≤C(hk+1+hβ0(d−1)−12)∥u∥k+2, (3.2) ∥Qhu−Ehu∥ ≤C(hk+2+hβ0(d−1)+12+hβ0(d−1)−1)∥u∥k+2.

### 3.1 Convergence of the semi-discrete scheme

Denote the error of the semi-discrete scheme (2.3) by . With the use of Lemma 3.1, error estimates can be derived as follows.

###### Theorem 3.2.

Let and be the exact solution of (1.1)-(1.3) and the numerical solution of (2.3), respectively. Assume , and where , then there exists a positive constant such that

 ∥eh(t)∥≤C(hk+2+hβ0(d−1)+12+hβ0(d−1)−1)(∥φ∥k+2+∥u(t)∥k+2+∫t0∥ut∥k+2ds), (3.3)

and

 \interleaveeh(t)\interleave≤ C(hk+1+hβ0(d−1)−12)(∥φ∥k+2+∥u(t)∥k+2+∫t0∥ut∥k+2ds). (3.4)

Moreover, optimal convergence orders appear when the penalty parameter satisfies .

###### Proof.

It is necessary to decompose into two items

 eh=(Qhu−Ehu)+(Ehu−uh):=ρ+η.

From Lemma 3.1, we just need to estimate . On account of the semi-discrete scheme (2.3), the definition of and that of , we get the following identity for each (see Gao2014ON ())

 (ηt,v0)+aw(η,v) =(Ehut,v0)−((u0)t,v0)+aw(Ehu,v)−aw(uh,v) (3.5) =(Ehut,v0)−(f,v0)+aw(Ehu,v)=−(ρt,v0).

Thus, choosing in the above identity yields

 12ddt(η,η)+aw(η,η)=−(ρt,η).

With Cauchy-Schwarz inequality and (2.5), integrating the above equation over on the both sides shows that

 ∥η(t)∥2≤∥η(0)∥2+C∫t0∥ρt∥2ds.

Applying triangle inequality and Lemma 3.1 results in the estimate (3.3).

Moreover, we estimate . By taking in (3.5), one obtains

 ∥ηt∥2+aw(η,ηt)=−(ρt,ηt).

Notice that with Cauchy-Schwarz inequality, it is easy to get

 12ddtaw(η,η)≤12∥ρt∥2,

and then integrating the inequality on leads to

 \interleaveη(t)\interleave2≤\interleaveη(0)\interleave2+∫t0∥ρt∥2ds.

Consequently, with the use of Lemma 3.1, (3.4) is proved. ∎

### 3.2 Convergence of the full-discrete scheme

For each , we denote the error term of full-discrete schemes by

 en:=Qhu(tn)−un=(Qhu(tn)−Ehu(tn))+(Ehu(tn)−un):=ρn+ηn.

Fully discrete error estimates are given in following theorem.

###### Theorem 3.3.

Let and be the exact solution of (1.1)-(1.3) and the numerical solution of (2.4), respectively. Assume , and with .
When , assume , then there exists a positive constant such that

 ∥en∥≤ C(hk+2+hβ0(d−1)+12+hβ0(d−1)−1)(∥φ∥k+2+∥u(tn)∥k+2 (3.6) +∫tn0∥ut∥k+2ds)+CτM1,

and

 \interleaveen\interleave≤ C(hk+1+hβ0(d−1)−12)(∥φ∥k+2+∥u(tn)∥k+2+∫tn0∥ut∥k+2ds)+CτM1, (3.7)

where .

When , assume , then there exists a positive constant such that

 ∥en∥≤ C(hk+2+hβ0(d−1)+12+hβ0(d−1)−1)(∥φ∥k+2+∥u(tn)∥k+2 (3.8) +∫tn0∥ut∥k+2ds)+Cτ2M2,

and

 \interleaveen\interleave≤ C(hk+1+hβ0(d−1)−12)(∥φ∥k+2+∥u(tn)∥k+2+∫tn0∥ut∥k+2ds)+Cτ2M2, (3.9)

where .
Here, optimal convergence orders appear when the penalty parameter satisfies .

###### Proof.

Based on Lemma 3.1, it is required to estimate in the energy norm. For each , combining scheme (2.4) with the definition (3.1), we have the following error equation

 (¯∂ηn,v0)+aw(θηn+(1−θ)ηn−1,v) (3.10) =(¯∂Ehu(tn),v0)−(¯∂un,v0)+aw(θEhu(tn)+(1−θ)Ehu(tn−1),v) −aw(θun+(1−θ)un−1,v) =(¯∂Ehu(tn),v0)−(∇⋅A(∇(θu(tn)+(1−θ)u(tn−1))),v0) −(θf(tn)+(1−θ)f(tn−1),v0) =(¯∂Ehu(tn),v0)−(θut(tn)+(1−θ)ut(tn−1),v0) =−(¯∂ρn,v0)+(¯∂Qhu(tn)−(θut(tn)+(1−θ)ut(tn−1)),v0).

By integration by parts, we can deduce the following identities when ,

 ¯∂u(tn)−(θut(tn)+(1−θ)ut(tn−1))=−1τ∫tntn−1(s−(1−θ)tn−θtn−1)uttds, (3.11)

and when ,

 ¯∂u(tn)−12(ut(tn)+ut(tn−1))=12τ∫tntn−1(tn−s)(tn−1−s)utttds. (3.12)

(i) In the case . Taking and substituting (3.11) into (3.10), it holds

 12τ∥ηn∥2−12τ∥ηn−1∥2+1τ(θ−12)∥ηn−ηn−1∥2+\interleaveθηn+(1−θ)ηn−1\interleave2 ≤−(¯∂ρn,θηn+(1−θ)ηn−1)+τ12(∫tntn−1∥utt∥2ds)12∥θηn+(1−θ)ηn−1∥.

Furthermore, with the use of Cauchy-Schwarz inequality, the above inequality can be rewritten as

 ∥ηn∥2−∥ηn−1∥2+2(θ−12)∥ηn−ηn−1∥2+2τ\interleaveθηn+(1−θ)ηn−1\interleave2 (3.13) ≤2τ∥¯∂ρn∥∥θηn+(1−θ)ηn−1∥+2τ32(∫tntn−1∥utt∥2ds)12∥θηn+(1−θ)ηn−1∥ ≤τε∥¯∂ρn∥2+2ετ∥θηn+(1−θ)ηn−1∥2+τ2ε∫tntn−1∥utt∥2ds.

Therefore, by (2.5) and , choosing in (3.13) leads to

 ∥ηn∥2−∥ηn−1∥2≤Cτ∥¯∂ρn∥2+Cτ2∫tntn−1∥utt∥2ds. (3.14)

By using and summing (3.14) from to , it holds

 ∥ηn∥2≤∥η0∥2+C∫tn0∥ρt∥2+Cτ2∫tn0∥utt∥2ds.

The error estimate (3.6) in -norm follows owing to and Lemma 3.1.

Concerning about the error estimate in energy norm, taking in (3.10) leads to

 1τ∥ηn−ηn−1∥2+(θ−12)\interleaveηn−ηn−1\interleave2+12(\interleaveηn\interleave2−\interleaveηn−1\interleave2) =−(¯∂ρn,ηn−ηn−1)+(¯∂Qhu(tn)−(θut(tn)+(1−θ)ut(tn−1)),ηn−ηn−1) ≤Cτ∥¯∂ρn∥2+Cτ2∫tntn−1∥utt∥2ds+ετ∥ηn−ηn−1∥2.

With and appropriate choice , we show that

 \interleaveηn\interleave2≤\interleaveηn−1\interleave2+Cτ∥¯∂ρn∥2+Cτ2∫tntn−1∥utt∥2ds. (3.15)

Finally, owing to Lemma 3.1, the result (3.7) follows immediately.

(ii) In the case , by applying (3.12) to the above process, the last term at the right hand side of (3.14) and that of (3.15) become . Analogously, we can also prove the results (3.8) and (3.9). ∎

## 4 Numerical experiments

In this section, we will give an example in 2D to verify our theory. Let and . We consider the backward Euler and CN schemes on time discretization, respectively. The error estimates are established in time level . In the example, we set is an identity matrix, denote the convergence order by , and take the penalty parameter . The programming is implemented in Matlab while uniform triangular meshes are generated by Gmsh. {example} The exact solution is as the same as in Li2013Weak () and the initial condition, the Dirichlet boundary condition and the source function are determined by exact solution. We first consider the backward Euler OPWG scheme. Table 1 and Table 2 show the numerical convergence with respect to mesh sizes while the time steps are taken small enough. When taking in Table 1 for , the convergence orders are and in the energy norm and -norm, respectively, which are in agreement with our analysis completely. Moreover, by choosing in Table 2 for , the convergence orders are and in the energy norm and -norm, respectively. On the other hand, Table 3 presents convergence orders about time step for . When mesh size is fixed, the convergence orders on time are in both energy norm and -norm. Moreover, we illustrate the errors of Table 3 in Fig. 1 with functions. The least squares fitting method is used to get convergence rates.