# A Priori and a Posteriori Error Estimates of a WOPSIP DG Method for the Heat Equation

Authors: Yuping Zeng; Kunwen Wen; Fen Liang; Huijian Zhu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

We introduce and analyze a weakly overpenalized symmetric interior penalty method for solving the heat equation. We first provide optimal a priori error estimates in the energy norm for the fully discrete scheme with backward Euler time-stepping. In addition, we apply elliptic reconstruction techniques to derive a posteriori error estimators, which can be used to design adaptive algorithms. Finally, we present two numerical experiments to validate our theoretical analysis.


---

## Body

## 1. Introduction

LetΩ⊂ℝ2 be a bounded polygonal domain with Lipschitz boundary ∂Ω, then we consider the following heat equation:(1)∂tu−Δu=f,inΩ×0,T,u=0,on∂Ω×0,T,u⋅,0=u0,where T>0 is finite time, f∈L20,T;L2Ω is the source term, and u0∈L2Ω stands for the initial data. There has been much research on the a priori and a posteriori error estimates of finite element methods (FEM) for parabolic equations (see [1–12] and the references therein). However, most of the literature on this subject considers only conforming (or nonconforming) FEM in space. The error estimate of discontinuous Galerkin (DG) methods for such problems is still very rare. In the context of DG methods for space variable, see [13–15] for a priori error analysis, and see [16–21] for a posteriori error analysis. Recently, Ern et al. [22, 23] have developed a posteriori error estimates for the parabolic problem with DG discretization in time (see [24] for a posteriori error estimates of nonconforming Crouzeix–Raviart FEM for the heat equation). The object of the present work is to investigate a weakly overpenalized symmetric interior penalty (WOPSIP) method in space for problem (1), combined with an implicit Euler scheme in time.The WOPSIP method is a kind of nonconsistent discontinuous Galerkin (DG) scheme, which was initially proposed in [25] for solving the second order elliptic equation, therein a priori error estimates were obtained. In recent years, DG methods have received significant attention since they are suitable for hp-adaptive computations. Besides this, they can also deal with nonhomogeneous boundary conditions and curved boundaries easily and allow for meshes with hanging nodes. Compared with the well-known DG methods in [26], the WOPSIP method has some advantages. For example, it has less computational complexity and thus is easy to implement [25]. Additionally, a high intrinsic parallelism property of the WOPSIP method was investigated in [27]. Therefore, the WOPSIP methods have been further developed to solve biharmonic problems [28] and Stokes equations [29], Reissner–Mindlin plate equations [30, 31], non-self-adjoint and indefinite problems [32], and variational inequalities [33]. In the present work, we shall extend the results of elliptic equations in [25] to the parabolic case. More precisely, the space variable is approximated by the WOPSIP method, and time variable is discretized by the backward Euler scheme. We shall give a detailed a priori and a posteriori error estimates. In this case, one may come across a difficulty that stems from the nonconsistency of the numerical method (for more details, see (34) in Theorem 2). On the contrary, a posteriori error analysis for parabolic problems is more involved than that for elliptic equations, since they involve both spatial error and temporal error. According to a framework stated in [18], we derive a posteriori error estimates which rely on the available estimates for elliptic problems [32, 34].The rest of our paper is organized as follows. In Section2, we state some notations and the numerical scheme. We establish a priori error estimates in the energy norm in Section 3. Section 4 is devoted to a posteriori error analysis which is based on the work [18]. We make some conclusions in Section 5. Finally, we provide some numerical results to validate theoretical analysis of a priori error estimates.


## 2. Preliminaries and WOPSIP Method

Let us first give some notations. For a bounded subdomainD⊂Ω, we denote by HmDm≥0 the standard Sobolev space, associated with norm ⋅m,D and seminorm ⋅m,D. When m=0, H0D is the standard Lebesgue space L2D, with the inner product defined by ⋅,⋅D. When D=Ω, we will omit the index Ω.To deal with functions of time and space, we also introduce the standard Bochner spaceLp0,T;HmD, which consists of all measurable functions u:0,T⟶HmD with norm(2)uLp0,T;HmD=∫0Tusm,Dpds1/p,for 1≤p<∞.The weak formulation of (1) reads: find u∈L20,T;H01Ω with ∂tu∈L20,T;H−1Ω such that(3)∂tu,v+au,v=f,v,∀v∈H01Ω,t∈0,T,u⋅,0=u0,with au,v=∫Ω∇u⋅∇vdx.LetTh be a family of conforming shape-regular meshes which decompose Ω into triangle elements K. Set hK=diamK and h=maxK∈ThhK. Denote by ℰh the set of all edges. Furthermore, ℰh=ℰhI∪ℰh∂, where ℰhI is the set of interior edges and ℰh∂ is the set of edges on ∂Ω. In what follows, he stands for the length of the edge e. Given D∈ℝ2, we let ℙrD be the space of polynomials of degree at most r on D. Moreover, we associate a fixed unit normal n with each edge e∈ℰh such that for edges on the boundary ∂Ω, n is the exterior unit normal vector.Lete be an interior edge in ℰhI shared by elements K+ and K−. For v:Ω⟶ℝ, set v±=ve∩∂K±, we define the following quantities:(4)v=v+−v−,v=12v++v−.Ife∈ℰh∂, set v=v and v=v. Furthermore, for m≥1, we also define(5)HmTh=v∈L2Ω:vK∈HmK,∀K∈Th.Consider the discontinuousℙ1 finite element space:(6)Vh=v∈L2Ω:vK∈ℙ1K,∀K∈Th.The bilinear form of the WOPSIP DG method is defined by (see [25])(7)ahw,v=∑K∈Th∫K∇w⋅∇vdx+∑e∈ℰhhe−2Πe0wΠe0v,where Πe0v stands for the mean of v over e, that is,(8)Πe0v=1he∫evdA.For the time discretization, we introduce the uniform partition0=t0<t1<⋯<tN=T of 0,T, with time step τ=T/N and tn=nτn=1,…,N. The backward Euler WOPSIP DG method for solving heat equation (1) is to find uhn∈Vh for n=1,…,N such that(9)uhn−uhn−1τ,vh+ahuhn,vh=fn,vh,∀vh∈Vh,uh0=uh0,where uh0 is a projection of u0 onto Vh which will be specified later in Sections 3 and 4.


## 3. Stability and A Priori Error Analysis

We begin by defining the mesh-dependent norm⋅h on H2Th as(10)vh=∑K∈Th∇v0,T2+∑e∈ℰhhe−2Πe0v21/2.For the subsequent analysis, we need the following Poincaré–Friedrichs inequality (see [13, 35]):(11)v02≤C∑K∈Th∇v0,T2+∑e∈ℰhhe−1v0,e2,∀v∈H2Th.Here and hereafter, we useC to denote a positive constant which is independent of h and τ, but may have different values at different places.Moreover, from Lemma 3.1 in [25] we have(12)∑e∈ℰhhe−1v0,e2≤Cvh2,∀v∈H2Th.Combining (11) and (12) gives(13)v0≤Cvh,∀v∈H2Th.We first prove that the numerical solutions satisfy the following stability result.Theorem 1.
Letuhnn≥0 be the solution of (9). It holds that, for all m>0,(14)uhm02+τ∑n=1muhnh2≤Cuh002+τ∑n=1mfn02.Proof.
Testingvh=uhn in (9) and using the definition of ⋅h, we have(15)1τuhn−uhn−1,uhn+uhnh2=fn,uhn.
Then, using the equationaa−b=1/2a2−b2+1/2a−b2 and Cauchy–Schwarz inequality gives(16)12τuhn02−uhn−102+12τuhn−uhn−102+uhnh2≤fn0uhn0,which together with the inequality (13) and Young’s inequality implies that(17)12τuhn02−uhn−102+12τuhn−uhn−102+uhnh2≤Cfn0uhnh≤12uhnh2+Cfn02.
Thus,(18)12τuhn02−uhn−102+12τuhn−uhn−102+12uhnh2≤Cfn02.
Multiplying (18) by 2τ and summing from n=1 to n=m, we arrive at(19)uhm02−uh002+∑n=1muhn−uhn−102+τ∑n=1muhnh2≤Cτ∑n=1mfn02.
Noting that∑n=1muhn−uhn−102≥0, we get(20)uhm02+τ∑n=1muhnh2≤uh002+Cτ∑n=1mfn02.
The conclusion (14) follows immediately.To carry out a detailed error analysis, we introduceℐhu which denotes the continuous linear interpolation of u. We then have the following standard estimates (see [36]):(21)u−ℐhum,K≤Ch2−mu2,K,m=0,1.In addition, we introduce the following trace inequality:(22)he−1w0,e2≤ChK−2w0,K2+∇w0,K2,∀w∈H1K,with e being an edge of K.Also, we need the elliptic projection operatorPh:H2Ω⟶Vh defined by(23)ahv−Phv,wh=0,∀wh∈Vh.It is well known that the projection operatorPh satisfies the following estimates (see [11]).Lemma 1.
For anyv∈H2Ω, it holds that(24)v−Phv0,Ω+hv−Phvh≤Ch2v2.For convenience, we use the following notation for any functiongt,x: at each time step tn=nτ,n=1,…,N,gn=gtn,x,∀x∈Ω. In addition, Taylor expansion yields(25)gn−gn−1Δt=∂tgn+1Δt∫tn−1tntn−1−s∂ttgsds.Now, we are in a position to state a priori error estimates, which is the main result of this section.Theorem 2.
Letu and uhnn≥0 be the solutions of (1) and (9), respectively. Assume that ∂tu∈L20,T;H2Ω, ∂ttu∈L20,T;L2Ω, and uh0 satisfying(26)u0−uh00≤Ch2u02.Then, for allm>0 it holds that(27)τ∑n=1mun−uhnh2≤Ch2τ∑n=1mun22+h2u022+h2∫0tm∂tus22ds+Cτ2∫0tm∂ttus02ds.Proof.
Integrating by parts and using the Taylor formulation (25), we have(28)un−un−1τ,vh+ahun,vh=fn,vh+Rn,vh+∑e∈ℰh∫e∇un⋅nvhdA,∀vh∈Vh,where(29)Rn=−1τ∫tn−1tns−tn−1∂ttusds.
We then split the errorun−uhn into(30)un−uhn=un−Phun+Phun−uhn=θn+ρn.
Thus, subtracting (9) from (28) and using the definition of θn and ρn gives(31)ρn−ρn−1τ,vh+ahρn,vh=∑e∈ℰh∫e∇un⋅nvhdA−θn−θn−1τ,vh+Rn,vh.
Testingvh=ρn in (31), using the definition of ⋅h and the formula aa−b=1/2a2−b2+1/2a−b2, we obtain(32)12τρn02−ρn−102+12τρn−ρn−102+ρnh2=∑e∈ℰh∫e∇un⋅nρndA+Rn,ρn−θn−θnτ,ρn≡S1+S2+S3.
We now bound the first termS1 of the above equation. We then employ the definition of Πe0 to find that(33)S1=∑e∈ℰh∫e∇un⋅nρndA=∑e∈ℰh∫e∇un−ℐhun⋅nρndA+∑e∈ℰh∫e∇ℐhun⋅nΠe0ρndA≡S11+S12.
Furthermore, applying Cauchy–Schwarz inequality and Young’s inequality and using (11), (21), and (22), we can estimate the terms S11 and S12 as follows:(34)S11=∑e∈ℰh∫e∇un−ℐhun⋅nρndA≤∑e∈ℰhhe∇un−ℐhun0,e21/2∑e∈ℰhhe−1ρn0,e21/2≤Chun2ρnh≤Ch2un22+18ρnh2,S12=∑e∈ℰh∫e∇ℐhun⋅nΠe0ρndA≤∑e∈ℰhhe3∇ℐhun0,e21/2∑e∈ℰhhe−3Πe0ρn0,e21/2=∑e∈ℰhhe3∇ℐhun0,e21/2∑e∈ℰhhe−2Πe0ρn21/2≤∑e∈ℰhhe3∇un−ℐhun0,e21/2∑e∈ℰhhe−2Πe0ρn21/2+∑e∈ℰhhe3∇un0,e21/2∑e∈ℰhhe−2Πe0ρn21/2≤Ch2un2ρnh+Chun1,Ω+h2un2,Ωρnh≤Chun2ρnh≤Ch2un22+18ρnh2.
Next, we give bounds forS2. Applying Cauchy–Schwarz inequality, the inequality (13) and Young’s inequality yields(35)S2=Rn,ρn≤Rn0ρn0≤CRn0ρnh≤CRn02+18ρnh2.
Similarly, it holds that(36)S3=−θn−θn−1τ,ρn≤θn−θn−1τ0ρn0≤Cθn−θn−1τ0ρnh≤Cθn−θn−1τ02+18ρnh2.
Combining (32)–(36), we obtain(37)12τρn02−ρn−102+12τρn−ρn−102+12ρnh2≤Ch2un22+Rn02+θn−θn−1τ02.
Moreover, applying Cauchy–Schwarz inequality forRn gives(38)Rn02≤1τ2∫tn−1tns−tn−12ds∫tn−1tnutts02ds=τ3∫tn−1tn∂ttus02ds.
On the contrary,(39)θn−θn−1=∫tn−1tn∂tθsds.
Thus,(40)θn−θn−1τ02=1τ2∫tn−1tn∂tθsds02≤1τ2∫tn−1tn12dt∫tn−1tn∂tθs02ds=1τ∫tn−1tn∂tθs02ds.
Plugging (38) and (40) into (37) and then multiplying the result inequality by 2τ and summing from n=1 to n=m, we obtain(41)ρm02−ρ002+∑n=1mρn−ρn−102+τ∑n=1mρnh2≤Ch2τ∑n=1mun22+τ2∫0tm∂ttus02ds+∫0tm∂tθs02ds.
Noting that∑n=1mρn−ρn−102≥0 and τ∑n=1mρnh2≥0, and using the estimate (24), we infer that(42)τ∑n=1mρnh2≤ρ002+Ch2τ∑n=1mun22+τ2∫0tm∂ttus02ds+h4∫0tm∂tus22ds.
Moreover, it follows from (24) and (26) that(43)ρ00=Phu0−uh00≤Phu0−u00+u0−uh00≤Ch2u02.
Plugging (43) into (42) and then combining (30), the triangle inequality yields the desired estimate in (27).Remark 1.
For simplicity, in the present work, we only consider the lowest linear finite element, and the extension of the results to high order methods [37] can be derived straightforwardly.Remark 2.
In this work, we only consider the constant coefficient parabolic equation, and the extension to more practical problems such as equations with variable coefficient can be derived by using the techniques developed in Section6 in [25]. Additionally, the present work only addresses the conforming mesh, and the extension to meshes with hanging nodes can be also obtained by utilizing approaches stated in Section 6 in [25].


## 4. A Posteriori Error Analysis

We begin by recalling an a posteriori error estimator for the stationary problem (see [32]), which is a key step in the subsequent error analysis.Theorem 3.
Letz∈H01Ω be the solution of(44)az,v=f,v,∀v∈H01Ω,and let Z∈Vh be the solution satisfying(45)ahZ,vh=f,vh,∀vh∈Vh.Then, we have the following a posteriori error estimates:(46)z−Zh≤CEZ,f,Th,with(47)ℰZ,f,Th=∑K∈ThhKf0,K+∑e∈ℰhIhe∇Z⋅n0,e2+∑e∈ℰhhe−2Πe0Z2+he−1Z0,e21/2.Remark 3.
It is worth mentioning that, the authors in [34] have proposed a different error estimator, which can also be applied directly in the forthcoming analysis.
The following useful fact is a standard result in studying a posteriori error estimates of DG methods. It shows that, given any discontinuous functionvh∈Vh, there exists a continuous polynomial function to approximate it (see [34, 38]).Lemma 2.
For anyvh∈Vh, there exits an decomposition vh=vhc+vhd such that(48)∑K∈Th∇vh−vhc0,K2+hK−2vh−vhc0,K2≤C∑e∈ℰhhe−1vh0,e2,where vhc∈H01Ω∩Vh and vhd∈Vh.To proceed, we reformulate the DG method (9) to make it be suitable for adaptive computations. At each time step n>0, we assume that a mesh Tn which can be obtained from Tn−1 by locally refining and coarsening Tn−1. Analogous to Section 3, we denote by ℰn=ℰnI∪ℰn∂ the set of edge of Tn, and use Vhn to stand for the finite element space with respect to Tn. Let In:Vhn−1⟶Vhn be a general data transfer operator, set τn=tn−tn−1, the backward Euler WOPSIP DG method for approximating (1) can be reformulated as: find uhn∈Vhnn=1,…,N such that(49)uhn−Inuhn−1τn,vhn+Bnuhn,vhn=fn,vhn,∀vhn∈Vhn,uh0=Ph0u0,where(50)Bnw,v=∑K∈Tn∫K∇w⋅∇vdx+∑e∈ℰnhe−2Πe0wΠe0v,and Ph0:L2Ω⟶Vh0 is the L2-projection operator onto Vh0. The corresponding mesh-dependent norm is defined by vTn2=Bnv,v. In view of uhn, we define a function which is piecewise linear continuous in time:(51)uh0=uh0,uht=ln−1tuhn−1+lntuhn,for t∈tn−1,tnn=1,…,N, where(52)lnt=t−tn−1τn1tn−1,tn+tn+1−tτn+11tn,tn+1,are the Lagrange basic functions.As in Lemma2, at each time step tnn=1,…,N, we can decompose uhn into(53)uhn=ucn+udn,with respect to the meshTˇn,(54)uhn=uc+n+ud+n,with respect to the meshTˇn+1,where Tˇn refers to the coarsest common refinement. More precisely, it is the coarsest triangulation which satisfies Tn≤Tˇn and Tn−1≤Tˇn. Here, we write Tn≤Tˇn to mean that Tˇn is a refinement of Tn. In addition, we define uhct=ln−1tuc+n−1+lntucn and uhdt=ln−1tud+n−1+lntudn.Moreover, we introduce the discrete elliptic operatorAn:Vhn⟶Vhn by(55)Anη,θ=Bnη,θ,∀θ∈Vhn,with Anη∈Vhn. Additionally, we define some elliptic reconstructions which are used commonly in the a posterior error analysis for parabolic equations [6].Definition 1.
The elliptic reconstructionwn∈H01Ω of uhn is defined as the solution of the elliptic problem:(56)Bnwn,v=Anuhn,v,∀v∈H01Ω.
Similarly,wn−1∈H01Ω satisfies(57)Bnwn−1,v=AnInuhn−1,v,∀v∈H01Ω.Remark 4.
We should point that the DG solution ofwn, denoted by whn, is also the DG solution uhn. Indeed, we have Bnwhn,vh=Anuhn,vh=Bnuhn,vh, for all vh∈Vhn; thus, whn=uhn.
Similar touht stated in (51), we define wht by(58)wt=ln−1twn−1+lntwn.
Now, the errore can be decomposed into(59)e=uh−u=φ−χ,φ=w−u,χ=w−uh,where φ is called the elliptic error and χ is the parabolic error. Furthermore, we define ec=uhc−u and χc=w−uhc. Thus, e=ec+uhd and χc=χ+uhd.
Then, we can obtain the following result, which is a key step to prove the main result in Theorem4. We remark that similar result can be found in Lemma 5.3 in [39].Lemma 3.
For eacht∈tn−1,tnn=1,2,…,N, we have(60)∂te,v+Bnφ,v=Inuhn−1−uhn−1τn,v+ln−1tAnInuhn−1−Anuhn,v+Phnfn−f,v,for all v∈H01Ω.Proof.
Since∂tuh=uhn−uhn−1/τn, this combined with (49), (56), and (57) yields(61)∂tuh,v=I−Phn∂tuh,v+Phn∂tuh,v=I−Phn∂tuh,v+Phn∂tuh,Phnv=I−Phn∂tuh,v+∂tuh,Phnv=I−Phn∂tuh,v+fn,Phnv+Inuhn−1−uhn−1τn,Phnv−Anuhn,Phnv=I−Phn∂tuh,v+Phnfn,v+Inuhn−1−uhn−1τn,Phnv−Bnwn,v,for all v∈H01Ω. Here, Phn stands for the L2-projection operator onto Vhn. Thus,(62)∂te,v+Bnφ,v=∂tuh−u,v+Bnw−u,v=∂tuh,v+Bnw,v−∂tu,v−Bnu,v=∂tuh,v+Bnw,v−f,v=I−Phn∂tuh,v+Phnfn,v+Inuhn−1−uhn−1/τn,Phnv−Bnwn,v+Bnw,v−f,v,where in the last line we have used (61).
In addition, for eacht∈tn−1,tnn=1,2,…,N, it follows from (56)–(58) that(63)Bnw−wn,v=Bnlntwn+ln−1twn−1−wn,v=t−tnτnBnwn−1−wn,v=t−tnτnAnInuhn−1−Anuhn,v.
On the contrary, we have(64)I−Phn∂tuh,v=I−Phnuhn−uhn−1τn,v=−1τnI−Phnuhn−1,v=−1τnuhn−1,v+1τnPhnuhn−1,v.
We also have(65)Inuhn−1−uhn−1τn,Phnv−v=−1τnuhn−1,Phnv−v=−1τnuhn−1,Phnv+1τnuhn−1,v.
Substituting (63)–(65) into (62) gives the desired result.
We then define the error estimators as follows:(1)
We define the time-stepping estimator as(66)αn=CP3Inuhn−1−uhnTn,withCP satisfying(67)PhnvTn≤CPvTn,for allv∈H01Ω.(2)
Set the data approximation estimator in time to be(68)βn=∫tn−1tnPhnfn−fsH−1Ω2τnds1/2.(3)
The estimator with respect to mesh-change is given by(69)ζn=Inuhn−1−uhn−1H−1Ωτn.(4)
The nonconforming part of parabolic estimator is defined as(70)ξn=udn−ud+n−1H−1Ωτn,and the nonconforming part of elliptic estimator is given by(71)θn=udnTn2+ud+n−1Tn21/2.(5)
The space error estimator is(72)ρn=ℰuhn,Anuhn,Tn,whereℰ is stated in (47) in Theorem 3. Similarly, we define(73)ρn−1=ℰInuhn−1,AnInuhn−1,Tn.(6)
We also define(74)ϵn=udn−ud+n0τn,which can be understood as the nonconforming parabolic part estimator of higher order.
We now state the main result of this section, which is largely based on the work [18]. For the sake of completeness, we sketch a proof.Theorem 4.
Letu and uht be solutions of (1) and (51), respectively. For each m=1,2,…,N, we have(75)∫0tmus−uhsTn2ds1/2≤u0−uhc00+2ηe,m+3ηp,m+12∑n=1mθn2τn1/2+3∑n=1m−1ϵnτn.Here,ηe,m refers to the elliptic error estimator which is defined as(76)ηe,m=∑n=1mρn−12+ρn2τn1/2,and ηp,m is parabolic error estimator defined by(77)ηp,m=∑n=1mαn+βn+ζn+ξn2τn1/2.Proof.
Selectingv=ec in (60) and noting that e=ec+uhd and φ=ec+χc, we deduce that(78)12ddtec02+φTn2=∂tec,ec+Bnφ,φ=−∂tuhd,ec+Bnφ,χc+ln−1tAnInuhn−1−Anuhn,ec+Inuhn−1−uhn−1τn+Phnfn−f,ec.
Integrating (78) on 0,tm implies that(79)12ectm02+12∫0tmφsTn2ds=12ec002+∫0tmBnφ,χcds+∑n=1m∫tn−1tnln−1sAnInuhn−1−Anuhn,ecds+∑n=1m∫tn−1tnInuhn−1−uhn−1τn+Phnfn−f,ecds−∫0tm∂tuhd,ecds+12∑n=1m−1utn−uc+n02−ectn02≡M0+M1+M2+M3+M4+M5.
The first termM1 can be bounded by(80)M1=∫0tmBnφ,χcds≤∫0tmφsTnχcsTnds.
To estimateM2, it follows from (67) that(81)AnInuhn−1−Anuhn,ec=AnInuhn−1−Anuhn,Phnec=BnInuhn−1−uhn,Phnec≤Inuhn−1−uhnTnPhnecTn≤CPInuhn−1−uhnTnecTn,this together with ec=φ−χc implies that(82)M2≤∑n=1m∫tn−1tnln−1sCPInuhn−1−uhnTnecTnds≤∑n=1m∫tn−1tnln−12sCP2Inuhn−1−uhnTn2ds1/2∫tn−1tnecTn2ds1/2=∑n=1mαnτn∫tn−1tnecTn2ds1/2≤∑n=1mαnτn∫tn−1tnφTn2ds1/2+∫tn−1tnχcTn2ds1/2.
For the termM3, we first have(83)M3=∑n=1m∫tn−1tnInuhn−1−uhn−1τn+Phnfn−f,ecds≤∑n=1m∫tn−1tnInuhn−1−uhn−1τnH−1Ω∇ec0ds+∑n=1m∫tn−1tnPhnfn−fH−1Ω∇ec0ds≡M31+M32.
We then estimateM31 as follows:(84)M31=∑n=1m∫tn−1tnInuhn−1−uhn−1τnH−1Ω∇ec0ds≤∑n=1m∫tn−1tnInuhn−1−uhn−1H−1Ω2τn2ds1/2∫tn−1tn∇ec02ds1/2=∑n=1mζnτn∫tn−1tn∇ec02ds1/2.
Similarly, the termM32 can be bounded by(85)M32=∑n=1m∫tn−1tnτn−1/2Phnfn−fH−1Ωτn1/2∇ec0ds≤∑n=1m∫tn−1tnPhnfn−fsH−1Ω2τnds1/2τn∫tn−1tn∇ec02ds1/2=∑n=1mβnτn∫tn−1tn∇ec02ds1/2.
Substituting (84) and (85) into (83) yields(86)M3≤∑n=1mζn+βnτn∫tn−1tn∇φ02ds1/2+∫tn−1tn∇χc02ds1/2.
For the termM4, we have(87)M4≤∑n=1m∫tn−1tn∂tuhdH−1Ω∇ec0ds≤∑n=1mξnτn∫tn−1tn∇φ02ds1/2+∫tn−1tn∇χc02ds1/2.
Observing that(88)utn−uc+n02−ectn02=udn−ud+n02+2udn−ud+n,ectn,we then have(89)M5=12∑n=1m−1utn−uc+n02−ectn02≤12∑n=1m−1ϵn2τn2+2ϵnτnmax1≤k≤m−1ectk0≤32∑n=1m−1ϵnτn2+14max1≤k≤m−1ectk02.
Adapting similar techniques as in the proof of Theorem 5.7 in [18], we can obtain the desired estimate (75).


## 5. Numerical Experiments

In this section, we give some numerical tests to validate our theoretical analysis in Theorem2. The numerical scheme (9) shows that, at each time step 1≤n≤N, we shall solve the following linear system: Given uh0=Ph0u0, find uhn∈Vh such that(90)τahuhn,vh+uhn,vh=τfn,vh+uhn−1,vh,∀vh∈Vh.In view of this expression, we then implemented our numerical experiments by MATLAB package.Example 1.
We consider problem (1) with Ω=0,1×0,1, T=1. We choose the source term f such that the exact solution is given by ux,t=sinπt/2x1−xy1−y.
The numerical results for the WOPSIP method are stated in Table1 and Figure 1. The convergence orders shown in Table 1 (which are also the slopes of the line in Figure 1) are consistent with the theoretical results in Theorem 2.Table 1
Numerical results fort = 1 with τ=h2, for Example 1.
h
∇
h
u
n
−
u
h
n
0
Order
u
n
−
u
h
n
0
Order
1
/
2
9.013994e − 2
—
1.833267e − 2
—
1
/
4
4.754053e − 2
0.9230
5.642120e − 3
1.7001
1
/
8
2.375353e − 2
1.0010
1.396086e − 3
2.0148
1
/
16
1.184181e − 2
1.0043
3.369024e − 4
2.0510
1
/
32
5.915072e − 3
1.0014
8.210111e − 05
2.0369Figure 1
The convergence rate fort = 1 with τ=h2 for Example 1.Example 2.
We consider problem (1) with Ω=0,1×0,1, T=1. The exact solution is ux,t=t2sinπxsinπy, and the source term f can be computed accordingly.
The convergence orders stated in Table2 (see also the slopes of the line in Figure 2) agree with the theoretical analysis in Theorem 2.Table 2
Numerical results fort = 1 with τ=h2 for Example 2.
h
∇
h
u
n
−
u
h
n
0
Order
u
n
−
u
h
n
0
Order
1
/
2
1.259907
—
1.386570e − 1
—
1
/
4
0.665702
0.9204
4.937369e − 2
1.4897
1
/
8
0.328806
1.0176
1.259566e − 2
1.9708
1
/
16
0.163091
1.0116
3.053313e − 3
2.0445
1
/
32
0.081347
1.0035
7.4407982e − 4
2.0368Figure 2
The convergence rate fort = 1 with τ=h2 for Example 2.


## 6. Conclusion and Future Work

A weakly overpenalized symmetric interior penalty method is proposed and analyzed for solving the heat equation. Optimal a priori error estimates in the energy norm are established. Moreover, we derive a posteriori error estimators which are key indicators to design adaptive algorithms. We only present some numerical tests to validate the a priori error estimate. A posteriori error estimators stated in Section4 include elliptic and parabolic terms. In particular, the error estimators with respect to the mesh-change operator In, such as αn and ςn, that need to be carefully investigated in the implementation procedure. Thus, the implementation of adaptive algorithms based on the proposed error estimators will be addressed in the further work.



---

*Source: 7525676-2020-04-17.xml*