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

TL;DR

This paper presents a novel weak Galerkin finite element method with implicit $ heta$-schemes for solving second-order parabolic problems, achieving optimal convergence rates and verified through numerical examples.

Contribution

It introduces a new weak Galerkin method with double-valued functions on edges and combines it with implicit $ heta$-schemes, including backward Euler and Crank-Nicolson, for improved accuracy.

Findings

Achieves optimal convergence rates in $L^2$ and energy norms.

Includes numerical verification of theoretical results.

Applicable to a range of $ heta$-schemes from backward Euler to Crank-Nicolson.

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 $(P_k(T), P_{k}(e), RT_k(T))$ element, a fully discrete approach is formulated with implicit $\theta$-schemes in time for $\frac{1}{2}\leq\theta\leq 1$, which include first-order backward Euler and second-order Crank-Nicolson schemes. Moreover, the optimal convergence rates in the $L^2$ and energy norms are derived. Numerical example is given to verify the theory.

Tables

Table 1. Table 1: Convergence with respect to mesh size h ℎ h in backward Euler scheme

h	$\tau =h^2,\theta =1,k=0,{\beta }_0=3$
	$⫴e_h⫴$	$\gamma$	$\Vert e_0\Vert$	$\gamma$
1/8	2.4624e-01		7.8378e-02
1/16	1.2250e-01	1.0072	2.0368e-02	1.9441
1/32	6.1139e-02	1.0026	5.1424e-03	1.9857
1/64	3.0554e-02	1.0007	1.2888e-03	1.9964
1/128	1.5275e-02	1.0001	3.2239e-04	1.9991

Table 2. Table 2: Convergence with respect to mesh size h ℎ h in backward Euler scheme

h	$\tau =h^3,\theta =1,k=1,{\beta }_0=5$
	$⫴e_h⫴$	$\gamma$	$\Vert e_0\Vert$	$\gamma$
1/4	1.5624e-01		2.6917e-02
1/8	3.9689e-02	1.9769	2.4406e-03	3.4632
1/16	9.9530e-03	1.9955	2.6886e-04	3.1823
1/32	2.4911e-03	1.9983	3.2656e-05	3.0414
1/64	6.2309e-04	1.9992	4.0350e-06	3.0167

Table 3. Table 3: Convergence with respect to time step τ 𝜏 \tau in backward Euler scheme

$\tau$	$h=1/128,\theta =1,k=1,{\beta }_0=5$
	$⫴e_h⫴$	$\sigma$	$\Vert e_0\Vert$	$\sigma$
$1/32$	1.8566e-02		1.8572e-02
$1/64$	9.6512e-03	0.9438	9.7364e-03	0.9316
$1/128$	4.9120e-03	0.9743	4.9825e-03	0.9665
$1/256$	2.4872e-03	0.9817	2.5281e-03	0.9788
$1/512$	1.2522e-03	0.9900	1.2667e-03	0.9969

Table 4. Table 4: Convergence with respect to mesh size h ℎ h of CN scheme

h	$\tau =h^2,\theta =\frac{1}{2},k=0,{\beta }_0=3$
	$⫴e_h⫴$	$\gamma$	$\Vert e_0\Vert$	$\gamma$
1/8	8.1480e-02		7.8963e-02
1/16	2.1661e-02	1.9113	2.0195e-02	1.9671
1/32	5.6280e-03	1.9444	5.0780e-03	1.9916
1/64	1.4812e-03	1.9258	1.2714e-03	1.9978
1/128	4.0315e-04	1.8773	3.1795e-04	1.9995

Table 5. Table 5: Convergence with respect to mesh size h ℎ h of CN scheme

h	$\tau =h^3,\theta =\frac{1}{2},k=1,{\beta }_0=5$
	$⫴e_h⫴$	$\gamma$	$\Vert e_0\Vert$	$\gamma$
1/4	3.5257e-02		2.7750e-02
1/8	8.5109e-03	2.0505	2.2400e-03	3.6309
1/16	2.5646e-03	1.7305	2.2664e-04	3.3050
1/32	6.7612e-04	1.9233	2.6522e-05	3.0951
1/64	1.7167e-04	1.9776	3.2647e-06	3.0221

Table 6. Table 6: Convergence with respect to time step τ 𝜏 \tau of CN scheme

$\tau$	$h=1/128,\theta =\frac{1}{2},k=1,{\beta }_0=5$
	$⫴e_h⫴$	$\sigma$	$\Vert e_0\Vert$	$\sigma$
$1/4$	1.2293e-01		1.2361e-01
$1/8$	1.5203e-02	3.0154	1.9252e-02	2.6827
$1/16$	2.5534e-03	2.5738	3.6999e-03	2.3794
$1/32$	5.8667e-04	2.1217	8.6930e-04	2.0895
$1/64$	1.4530e-04	2.0135	2.1619e-04	2.0075

Equations

$$u_t-\nabla \cdot (A\nabla u)$$

$$u$$

$$u(x,0)$$

$$(Aw,v)$$

$$(Av,v)$$

$$V_h:=$$

$$v_b{|}_e\in P_k(e),e\in \partial \Omega ,k\ge 0\},$$

$$({\nabla }_{w}v,q{)}_T=-(v_0,\nabla \cdot q{)}_T+⟨v_b,q\cdot \mathbf{n}{⟩}_{\partial T},\forall q\in R{T}_k(T).$$

$${\nabla }_w(Q_{h}v)=R_h(\nabla v),\forall v\in H^1(T).$$

$$(\nabla \cdot \mathbf{q},v_0{)}_T=(\nabla \cdot {\Pi }_h\mathbf{q},v_0{)}_T,\forall v_0\in P_k(T).$$

$$\Vert {\Pi }_h(A\nabla u)-A{\nabla }_w(Q_{h}u)\Vert \le C{h}^{k+1}\Vert u{\Vert }_{k+2}.$$

$$a_w(v,\chi ):=(A{\nabla }_{w}v,{\nabla }_w\chi )+J_0(v,\chi ),\forall v,\chi \in V_h,$$

$$J_0(v,\chi ):=\sum _{e\in {\mathcal{E}}_I}|e{|}^{-{\beta }_0}⟨[[v_b]],[[{\chi }_b]]{⟩}_e,{\beta }_0\ge 1.$$

$$((u_0{)}_t,v_0)+a_w(u_h,v_h)=(f,v_0),\forall v_h\in V_h^0.$$

$$\text{interleave}v{\text{interleave}}^2:=a_w(v,v).$$

$$(\bar{\partial }{u}^n,v_0)+a_w(\theta u^n+(1-\theta )u^{n-1},v_h)$$

$$=(\theta f(t^n)+(1-\theta )f(t^{n-1}),v_0),\forall v_h\in V_h^0,$$

$${\int }_{t-\nabla t}^{t+\nabla t}{\int }_K{u}_{t}dxdt+{\int }_{t-\nabla t}^{t+\nabla t}{\int }_{\partial K}\mathbf{q}\cdot \mathbf{n}dsdt={\int }_{t-\nabla t}^{t+\nabla t}{\int }_{K}fdxdt,$$

$${\int }_{t-\nabla t}^{t+\nabla t}{\int }_K(u_h{)}_{t}dxdt$$

$$\Vert v_0\Vert \le C\text{interleave}v\text{interleave}.$$

$$\Vert v_0{\Vert }^2$$

$$=-\sum _{T\in {\mathcal{T}}_h}({\nabla }_{w}v,{\Pi }_h\mathbf{q}{)}_T+\sum _{e\in {\mathcal{E}}_I}⟨[[v_b]],{\Pi }_h\mathbf{q}\cdot {\mathbf{n}}_{\mathbf{e}}{⟩}_e$$

$$\le \Vert {\nabla }_{w}v\Vert \Vert {\Pi }_h\mathbf{q}\Vert +\sum _{e\in {\mathcal{E}}_I}\Vert [[v_b]]{\Vert }_e\Vert {\Pi }_h\mathbf{q}{\Vert }_e$$

$$\displaystyle\leq\|\nabla_{w}v\|\|\Pi_{h}\mathbf{q}\|+\Big{(}\sum_{e\in\mathcal{E}_{I}}\|\llbracket v_{b}\rrbracket\|_{e}\Big{)}h^{-\frac{1}{2}}\|\Pi_{h}\mathbf{q}\|,$$

$$\Vert {\Pi }_h\mathbf{q}\Vert$$

$$\le Ch\Vert \Psi {\Vert }_2+\Vert R_h\nabla \Psi \Vert$$

$$\le C\Vert v_0\Vert .$$

$$\Vert v_0\Vert$$

$$\displaystyle\leq C\Big{(}\|\nabla_{w}v\|+h^{\frac{\beta_{0}(d-1)-1}{2}}\sum_{e\in\mathcal{E}_{I}}|e|^{-\frac{\beta_{0}}{2}}\|\llbracket v_{b}\rrbracket\|_{e}\Big{)}$$

$$\le C\text{interleave}v\text{interleave}.$$

$$\Vert u^n\Vert \le \Vert u^0\Vert +C\underset{t\in [0,\bar{T}]}{\sup }\Vert f(t)\Vert .$$

$$v_h=\theta u^n+(1-\theta )u^{n-1}=(\theta -\frac{1}{2})(u^n-u^{n-1})+\frac{1}{2}(u^n+u^{n-1})$$

$$\frac{1}{2}\Vert u^n{\Vert }^2-\frac{1}{2}\Vert u^{n-1}{\Vert }^2+(\theta -\frac{1}{2})\Vert u^n-u^{n-1}{\Vert }^2+\tau \text{interleave}\theta u^n+(1-\theta )u^{n-1}{\text{interleave}}^2$$

$$=\tau (\theta f(t^n)+(1-\theta )f(t^{n-1}),\theta u^n+(1-\theta )u^{n-1})$$

$$\le \frac{\tau }{4\epsilon }\Vert \theta f(t^n)+(1-\theta )f(t^{n-1}){\Vert }^2+\epsilon \tau \text{interleave}\theta u^n+(1-\theta )u^{n-1}{\text{interleave}}^2,$$

$$\Vert u^n{\Vert }^2\le \Vert u^{n-1}{\Vert }^2+C\tau \Vert \theta f(t^n)+(1-\theta )f(t^{n-1}){\Vert }^2.$$

$$\Vert u^n{\Vert }^2$$

Full text

1 Introduction

In this paper, an extension of weak Galerkin finite element method (WG) in [1] 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 ([1], [2], [3]), double-valued weak functions appeared in [4] has been employed to strengthen flexibility of WG with $(P_{k}(T),P_{k}(e),RT_{k}(T))$ 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 $\Omega\in\mathbb{R}^{d}~{}(d=2,3)$ be an open and bounded polygonal or polyhedral domain. A linear parabolic model is listed as follows

[TABLE]

where the functions $f(x,t)$, $g(x,t)$ and $\varphi(x)$ are known in some specific spaces for well-posedness. The coefficient matrix $A(x)$ is symmetric positive, i.e., there exist two positive constants $\alpha_{1}$ and $\beta_{1}$ such that for each $w,v\in\mathbb{R}^{d}$

[TABLE]

With different approximation spaces for weak gradient operator, WG with $(P_{k}(T),P_{k}(e),RT_{k}(T)),~{}k\geq 0$ element and element $(P_{k+1}(T),P_{k+1}(e),[P_{k}(T)]^{d})$ were developed for the parabolic equations in [5] and [6], respectively, first-order backward Euler full-discrete scheme being investigated. However, there are few publications on second-order fully discrete WG schemes. Based on $(P_{k}(T),P_{k}(e),RT_{k}(T))$ element [7] and $\theta$-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. $2$, the semi-discrete and full-discrete OPWG schemes are introduced and the latter is unconditionally stable. In Sec. $3$, optimal convergence analysis is presented including error estimates in the $L^{2}$ and energy norms. Finally, numerical results demonstrate the efficiency and feasibility of the new method.

Throughout this paper, we denote by $\varepsilon$ an arbitrarily small positive constant, $\|\cdot\|$ the $L^{2}$-norm and $L^{p}(0,\bar{T};V)$ with $p\geq 1$ the spaces with respect to time where $V$ represents Sobolev space (see details in [8] or [9]). Moreover, we use $C$ for a positive constant independent of mesh size $h$ and time step $\tau$.

2 OPWG schemes and stability

Let $\mathcal{T}_{h}$ be a partition of domain $\Omega$ satisfying shape regularity in [10]. For each element $T\in\mathcal{T}_{h}$, $h_{T}$ is its diameter and $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ is the mesh size of $\mathcal{T}_{h}$. Denote by $\mathcal{E}_{I}$ the set of interior edges or flat faces, and $\partial T$ the edges or flat faces of element $T$. Let ${P}_{k}(T)$ be the space of polynomials of degree less than or equal to $k$ in variables. The weak Galerkin finite element space for OPWG is defined as

[TABLE]

in particular, $V_{h}^{0}=\{v\in V_{h}~{}\mbox{and}~{}v_{b}=0~{}\mbox{on}~{}\partial\Omega\}$. For each $v=\{v_{0},v_{b}\}\in V_{h}$, we define a unique local weak gradient $\nabla_{w}v\in RT_{k}(T)$ on each element $T\in\mathcal{T}_{h}$ satisfying

[TABLE]

In addition, we define several local $L^{2}$ projection operators onto space $V_{h}$. For each $T\in\mathcal{T}_{h}$ and $e\in\partial T$, let $Q_{0}$, $Q_{b}$ be the $L^{2}$ projection operators to $P_{k}(T)$ and $P_{k}(e)$, respectively. Denote $Q_{h}v$ by $Q_{h}v=\{Q_{0}v,Q_{b}v\},~{}\forall v\in L^{2}(T)$. Meantime, define $R_{h}$ the $L^{2}$ projection onto $RT_{k}(T)$, and then one can obtain

[TABLE]

Furthermore, we define a div projection $\Pi_{h}$ for $\mathbf{q}\in H(div;\Omega)$ satisfying that $\Pi_{h}\mathbf{q}\in H(div;\Omega)$ and $\Pi_{h}\mathbf{q}\in RT_{k}(T)$ on each element $T$, and (see [1] and [11])

[TABLE]

Then, an approximation property of the projection $\Pi_{h}$ is given.

Lemma 2.1**.**

[1]** For $u\in H^{k+2}(\Omega),~{}k\geq 0$, it holds

[TABLE]

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

[TABLE]

where the penalty term is well defined as

[TABLE]

Let $e\in\mathcal{E}_{I}$ be shared by adjacent elements $T_{1}$ and $T_{2}$, then we define the jump on $e$ by $\llbracket v_{b}\rrbracket=v_{b}|_{T_{1}\cap e}-v_{b}|_{T_{2}\cap e}$.

The semi-discrete OPWG scheme for (1.1)-(1.3) is to seek $u_{h}(t)\in V_{h}$ satisfying the boundary condition $u_{h}(x,t)=Q_{h}g(x,t)$ on $\partial\Omega\times(0,\bar{T}]$ and the initial condition $u_{h}(0)=Q_{h}\varphi$ such that

[TABLE]

Now, we define energy norm as for any $v\in V_{h}$

[TABLE]

The existence and uniqueness of semi-discrete solution of (2.3) are obtained from coercivity and continuity of $a_{w}$ (see Lemma $3.1$ in [4]).

Next, we present full-discrete OPWG schemes. The interval $(0,\bar{T}]$ is divided into subintervals by time step $\tau$ uniformly, i.e. $t^{n}=n*\tau$. With the $\theta$-schemes applied, the full-discrete OPWG schemes are to seek $u^{n}\in V_{h}$ satisfying the boundary condition $u^{n}=Q_{h}g(x,t^{n})$ on $\partial\Omega\times(0,\bar{T}]$ and the initial condition $u^{0}=Q_{h}\varphi$ such that

[TABLE]

where the parameters $\theta$ vary in $[\frac{1}{2},1]$ and difference quotient $\bar{\partial}u^{n}:=\frac{u^{n}-u^{n-1}}{\tau}$. For simplification, we denote $f(t^{n}):=f(x,t^{n})$. Moreover, when $\theta=1$, (2.4) is backward Euler scheme, and Crank-Nicolson (CN) scheme is recovered if $\theta=\frac{1}{2}$.

Let $K\in\Omega$ be a small subdomain. The flux in time interval $(t-\nabla t,t+\nabla t)$ holds

[TABLE]

where $\mathbf{q}=-A\nabla u$ is the flow rate of heat energy. Multiplying a test function $v=\{v_{0},v_{b}=0\}$ such that $v_{0}=1$ in $K$ and $v_{0}=0$ elsewhere in (2.3), we can obtain that

[TABLE]

which implies mass conservation, by taking a numerical flux $\mathbf{q}_{h}\cdot\mathbf{n}=-R_{h}(A\nabla_{w}u_{h})\cdot\mathbf{n}$.

2.1 Stability of the full-discrete scheme

At first, we will give the following Poincar$\acute{\mbox{e}}$-type inequality between the $L^{2}$ norm and the energy norm.

Lemma 2.2**.**

For any $v\in V_{h}^{0}$, it holds

[TABLE]

Proof.

Based on Theorem $2.1$ in [12], we know that the weak solution $\Psi\in H^{1}_{0}(\Omega)$ of elliptic problem $\Delta\Psi=v_{0}$ satisfies $H^{2}$-regularity, i.e. $\|\Psi\|_{2}\leq C\|v_{0}\|$. Denote by $\mathbf{q}=\nabla\Psi\in H(div;\Omega)$ and it is obvious that $\nabla\cdot\mathbf{q}=v_{0}$. With the use of the definitions of $\Pi_{h}$ (2.1) and discrete weak gradient, trace inequality [10] and Cauchy-Schwarz inequality, then we can deduce that

[TABLE]

where $\mathbf{n}_{e}$ is a unit normal on $e$. From Lemma 2.1, it follows

[TABLE]

Combining (2.6) with (2.7) leads to

[TABLE]

∎

Theorem 2.3**.**

Let $u^{n}$ be the numerical solution of (2.4). Assume $g=0$ i.e. the parabolic problem is homogeneous problem and $\|f(t)\|$ is bounded in $[0,\bar{T}]$. Then there exists a positive constant such that

[TABLE]

Proof.

By Cauchy-Schwarz inequality and (2.5), taking

[TABLE]

in (2.4) yields

[TABLE]

where $\varepsilon>0$. Let $\varepsilon=\frac{1}{4}$ in (2.8), then it follows

[TABLE]

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

[TABLE]

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 $v\in H^{k+2}(\Omega)$, we define an elliptic projection $E_{h}v\in V_{h}$ satisfying the following equation

[TABLE]

where $E_{h}v$ is the $L^{2}$ projection of the trace of $v$ on the boundary. Then $E_{h}v$ is the OPWG approximation of the solution of the elliptic problem

[TABLE]

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

Lemma 3.1**.**

[4]** Let $u\in H^{k+2}(\Omega),~{}k\geq 0$, then there exists a positive constant $C$ such that

[TABLE]

3.1 Convergence of the semi-discrete scheme

Denote the error of the semi-discrete scheme (2.3) by $e_{h}:=Q_{h}u-u_{h}$. With the use of Lemma 3.1, error estimates can be derived as follows.

Theorem 3.2**.**

Let $u(t)$ and $u_{h}(t)\in V_{h}$ be the exact solution of (1.1)-(1.3) and the numerical solution of (2.3), respectively. Assume $u\in L^{1}(0,\bar{T};H^{k+2}(\Omega))$, $u_{t}\in L^{1}(0,\bar{T};H^{k+2}(\Omega))$ and $\varphi\in H^{k+2}(\Omega)$ where $k\geq 0$, then there exists a positive constant $C$ such that

[TABLE]

and

[TABLE]

Moreover, optimal convergence orders appear when the penalty parameter satisfies $\beta_{0}(d-1)\geq 2k+3$.

Proof.

It is necessary to decompose $e_{h}$ into two items

[TABLE]

From Lemma 3.1, we just need to estimate $\eta$. On account of the semi-discrete scheme (2.3), the definition of $E_{h}$ and that of $Q_{h}$, we get the following identity for each $v\in V_{h}^{0}$ (see [6])

[TABLE]

Thus, choosing $v=\eta$ in the above identity yields

[TABLE]

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

[TABLE]

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

Moreover, we estimate $\interleave e_{h}(t)\interleave$. By taking $v=\eta_{t}$ in (3.5), one obtains

[TABLE]

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

[TABLE]

and then integrating the inequality on $(0,t)$ leads to

[TABLE]

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

3.2 Convergence of the full-discrete scheme

For each $t^{n}\in(0,\bar{T}]$, we denote the error term of full-discrete schemes by

[TABLE]

Fully discrete error estimates are given in following theorem.

Theorem 3.3**.**

*Let $u$ and $u^{n}$ be the exact solution of (1.1)-(1.3) and the numerical solution of (2.4), respectively. Assume $u\in L^{1}(0,\bar{T};H^{k+2}(\Omega))$, $u_{t}\in L^{1}(0,\bar{T};H^{k+2}(\Omega))$ and $\varphi\in H^{k+2}(\Omega)$ with $k\geq 0$.

When $\frac{1}{2}<\theta\leq 1$, assume $u_{tt}\in L^{1}(0,\bar{T};L^{2}(\Omega))$, then there exists a positive constant $C$ such that*

[TABLE]

and

[TABLE]

where $M1=\int_{0}^{t^{n}}\|u_{tt}\|ds$.

When $\theta=\frac{1}{2}$, assume $u_{ttt}\in L^{1}(0,\bar{T};L^{2}(\Omega))$, then there exists a positive constant $C$ such that

[TABLE]

and

[TABLE]

*where $M2=\int_{0}^{t^{n}}\|u_{ttt}\|ds$.

Here, optimal convergence orders appear when the penalty parameter satisfies $\beta_{0}(d-1)\geq 2k+3$.*

Proof.

Based on Lemma 3.1, it is required to estimate $\eta^{n}$ in the energy norm. For each $v\in V_{h}^{0}$, combining scheme (2.4) with the definition (3.1), we have the following error equation

[TABLE]

By integration by parts, we can deduce the following identities when $\frac{1}{2}<\theta\leq 1$,

[TABLE]

and when $\theta=\frac{1}{2}$,

[TABLE]

(i) In the case $\frac{1}{2}<\theta\leq 1$. Taking $v=\theta\eta^{n}+(1-\theta)\eta^{n-1}$ and substituting (3.11) into (3.10), it holds

[TABLE]

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

[TABLE]

Therefore, by (2.5) and $\frac{1}{2}<\theta\leq 1$, choosing $\varepsilon=\frac{1}{4}$ in (3.13) leads to

[TABLE]

By using $\|\bar{\partial}\rho^{n}\|^{2}\leq\frac{C}{\tau}\int_{t^{n-1}}^{t^{n}}\|(Q_{h}-E_{h})u_{t}\|^{2}ds$ and summing (3.14) from $1$ to $n$, it holds

[TABLE]

The error estimate (3.6) in $L^{2}$-norm follows owing to $\|\eta^{0}\|=\|\rho^{0}\|$ and Lemma 3.1.

Concerning about the error estimate in energy norm, taking $v=\eta^{n}-\eta^{n-1}$ in (3.10) leads to

[TABLE]

With $\theta>\frac{1}{2}$ and appropriate choice $\varepsilon=\frac{1}{2}$, we show that

[TABLE]

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

(ii) In the case $\theta=\frac{1}{2}$, 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 $C\tau^{4}\int_{t^{n-1}}^{t^{n}}\|u_{ttt}\|^{2}ds$. 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 $\Omega=(0,1)\times(0,1)$ and $\bar{T}=1$. We consider the backward Euler $(\theta=1)$ and CN $(\theta=\frac{1}{2})$ schemes on time discretization, respectively. The error estimates are established in time level $t^{n}=\bar{T}$. In the example, we set $A\in\mathbb{R}^{2}$ is an identity matrix, denote the convergence order by $O(h^{\gamma}+\tau^{\sigma})$, and take the penalty parameter $\beta_{0}=2k+3$. The programming is implemented in Matlab while uniform triangular meshes are generated by Gmsh.

Example 1**.**

The exact solution is $u=sin(2\pi(t^{2}+1)+\pi/2)sin(2\pi x+\pi/2)sin(2\pi y+\pi/2)$ as the same as in [5] and the initial condition, the Dirichlet boundary condition and the source function $f$ 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 $h$ while the time steps are taken small enough. When taking $\tau=h^{2}$ in Table 1 for $k=0$, the convergence orders are $O(h)$ and $O(h^{2})$ in the energy norm and $L^{2}$-norm, respectively, which are in agreement with our analysis completely. Moreover, by choosing $\tau=h^{3}$ in Table 2 for $k=1$, the convergence orders are $O(h^{2})$ and $O(h^{3})$ in the energy norm and $L^{2}$-norm, respectively. On the other hand, Table 3 presents convergence orders about time step $\tau$ for $k=1$. When mesh size $h=1/128$ is fixed, the convergence orders on time are $O(\tau)$ in both energy norm and $L^{2}$-norm. Moreover, we illustrate the errors of Table 3 in Fig. 1 with $loglog$ functions. The least squares fitting method is used to get convergence rates.

Next, we apply the full-discrete CN scheme. Table 4 shows that when $\tau=h^{2}$ and $k=0$, the convergence orders are $O(h^{2})$ in the energy and $L^{2}$ norms for CN scheme. It is interesting that surperconvergence results in the energy norm are observed in Table 4. Moreover, with $\tau=h^{3}$ taken in Table 5 for $k=1$, the convergence orders are $O(h^{2})$ and $O(h^{3})$ in the energy and $L^{2}$ norms, respectively. In Table 6, we consider the convergence order for CN scheme on time step $\tau$ while $k=1$. When the mesh size $h=1/128$ is fixed enough fine, the convergence orders are $O(\tau^{2})$ in both energy and $L^{2}$ norms, which are in agreement with our theory. In Fig. 2, the errors in Table 6 are plotted.

Acknowledgements

The author wishes to thank associate professor Lunji Song for his critical reading of the manuscript, helpful discussions and valuable suggestions. The research of the author is supported in part by the Natural Science Foundation of Gansu Province, China (Grant 18JR3RA290).