On the superconvergence of a hydridizable discontinuous Galerkin method for the Cahn-Hilliard equation

TL;DR

This paper introduces a hybridizable discontinuous Galerkin (HDG) method for the Cahn-Hilliard equation, demonstrating optimal convergence and superconvergence properties, supported by theoretical analysis and numerical validation.

Contribution

The work develops a novel HDG Sobolev inequality and proves superconvergence of scalar variables in the HDG method for nonlinear PDEs.

Findings

Optimal convergence rates in $L^2$ norm for all variables.

Superconvergence of scalar variables in globally coupled degrees of freedom.

Numerical results confirm theoretical convergence rates.

Abstract

We propose a hydridizable discontinuous Galerkin (HDG) method for solving the Cahn-Hilliard equation. The temporal discretization can be based on either the backward Euler method or the convex-splitting method. We show that the fully discrete scheme admits a unique solution, and we establish optimal convergence rates for all variables in the $L^2$ norm for arbitrary polynomial orders. In terms of the globally coupled degrees of freedom, the scalar variables are superconvergent. Another theoretical contribution of this work is a novel HDG Sobolev inequality that is useful for HDG error analysis of nonlinear problems. Numerical results are reported to confirm the theoretical convergence rates.

Tables

Table 1. Table 1: Example   6.1 , k = 0 𝑘 0 k=0 with fully implicit scheme: Errors, observed convergence orders for u 𝑢 u , ϕ italic-ϕ \phi and their fluxes 𝒒 𝒒 \bm{q} and 𝒑 𝒑 \bm{p} .

$h/\sqrt{2}$	1/4	1/8	1/16	1/32	1/64
${\Vert 𝒒-{𝒒}_h\Vert }_{{𝒯}_h}$	8.6745E-04	4.8767E-04	2.5058E-04	1.2614E-04	6.3177E-05
order	-	0.83088	0.96061	0.99025	0.99757
${\Vert 𝒑-{𝒑}_h\Vert }_{{𝒯}_h}$	8.9032E-04	4.9104E-04	2.5102E-04	1.2620E-04	6.3184E-05
order	-	0.85847	0.96803	0.99214	0.99804
${\Vert u-u_h\Vert }_{{𝒯}_h}$	2.5400E-04	6.6287E-05	1.6746E-05	4.1975E-06	1.0501E-06
order	-	1.9380	1.9849	1.9962	1.9990
${\Vert \varphi -{\varphi }_h\Vert }_{{𝒯}_h}$	2.6147E-04	6.7626E-05	1.7040E-05	4.2683E-06	1.0676E-06
order	-	1.9510	1.9887	1.9972	1.9993

Table 2. Table 2: Example   6.1 , k = 1 𝑘 1 k=1 with fully implicit scheme: Errors, observed convergence orders for u 𝑢 u , ϕ italic-ϕ \phi and their fluxes 𝒒 𝒒 \bm{q} and 𝒑 𝒑 \bm{p} .

$h/\sqrt{2}$	1/4	1/8	1/16	1/32	1/64
${\Vert 𝒒-{𝒒}_h\Vert }_{{𝒯}_h}$	1.6623E-04	4.5233E-05	1.1599E-05	2.9202E-06	7.3141E-07
order	-	1.8778	1.9634	1.9899	1.9973
${\Vert 𝒑-{𝒑}_h\Vert }_{{𝒯}_h}$	1.6700E-04	4.5276E-05	1.1602E-05	2.9204E-06	7.3142E-07
order	-	1.8830	1.9644	1.9901	1.9974
${\Vert u-u_h\Vert }_{{𝒯}_h}$	4.8698E-05	6.1714E-06	7.7349E-07	9.6742E-08	1.2094E-08
order	-	2.9802	2.9962	2.9992	2.9998
${\Vert \varphi -{\varphi }_h\Vert }_{{𝒯}_h}$	4.9152E-05	6.1862E-06	7.7391E-07	9.6753E-08	1.2095E-08
order	-	2.9901	2.9988	2.9998	3.0000

Table 3. Table 3: Example   6.1 , k = 0 𝑘 0 k=0 with energy-splitting scheme: Errors, observed convergence orders for u 𝑢 u , ϕ italic-ϕ \phi and their fluxes 𝒒 𝒒 \bm{q} and 𝒑 𝒑 \bm{p} .

$h/\sqrt{2}$	1/4	1/8	1/16	1/32	1/64
${\Vert 𝒒-{𝒒}_h\Vert }_{{𝒯}_h}$	8.6761E-04	4.8768E-04	2.5059E-04	1.2614E-04	6.3177E-05
order	-	0.83111	0.96063	0.99025	0.99757
${\Vert 𝒑-{𝒑}_h\Vert }_{{𝒯}_h}$	8.9460E-04	4.9143E-04	2.5107E-04	1.2620E-04	6.3185E-05
order	-	0.86427	0.96891	0.99232	0.99809
${\Vert u-u_h\Vert }_{{𝒯}_h}$	2.5759E-04	6.7122E-05	1.6952E-05	4.2490E-06	1.0629E-06
order	-	1.9402	1.9853	1.9963	1.9991
${\Vert \varphi -{\varphi }_h\Vert }_{{𝒯}_h}$	2.6295E-04	6.7806E-05	1.7076E-05	4.2768E-06	1.0697E-06
order	-	1.9553	1.9894	1.9974	1.9993

Table 4. Table 4: Example   6.1 , k = 1 𝑘 1 k=1 with energy-splitting scheme: Errors, observed convergence orders for u 𝑢 u , ϕ italic-ϕ \phi and their fluxes 𝒒 𝒒 \bm{q} and 𝒑 𝒑 \bm{p} .

$h/\sqrt{2}$	1/4	1/8	1/16	1/32	1/64
${\Vert 𝒒-{𝒒}_h\Vert }_{{𝒯}_h}$	1.5809E-04	4.3945E-05	1.1415E-05	2.8955E-06	7.2935E-07
order	-	1.8470	1.9448	1.9790	1.9891
${\Vert 𝒑-{𝒑}_h\Vert }_{{𝒯}_h}$	1.5896E-04	4.3991E-05	1.1418E-05	2.8957E-06	7.2940E-07
order	-	1.8534	1.9459	1.9793	1.9891
${\Vert u-u_h\Vert }_{{𝒯}_h}$	4.9741E-05	6.3026E-06	7.9008E-07	9.8850E-08	1.2358E-08
order	-	2.9804	2.9959	2.9987	2.9998
${\Vert \varphi -{\varphi }_h\Vert }_{{𝒯}_h}$	4.9111E-05	6.1809E-06	7.7336E-07	9.6709E-08	1.2090E-08
order	-	2.9902	2.9986	2.9994	2.9998

Equations

$$u_t-\Delta \varphi$$

$$-ϵ\Delta u+{ϵ}^{-1}f(u)$$

$$\nabla u\cdot \mathbf{n}=\nabla \varphi \cdot \mathbf{n}$$

$$u(\cdot ,0)$$

$$\Delta t\sum _{n=1}^N\Vert {\mathbf{p}}^n-{\mathbf{p}}_h^n{\Vert }_{L^2(\Omega )}^2+\Delta t\sum _{n=1}^N\Vert {\mathbf{q}}^n-{\mathbf{q}}_h^n{\Vert }_{L^2(\Omega )}^2\le C(h^{k+1}+\Delta t{)}^2,$$

$$\Delta t\sum _{n=1}^N\Vert {\varphi }^n-{\varphi }_h^n{\Vert }_{L^2(\Omega )}^2+\underset{1\le n\le N}{\max }\Vert u^n-u_h^n{\Vert }_{L^2(\Omega )}^2\le C(h^{k+2}+\Delta t{)}^2.$$

$$(w,v{)}_{{\mathcal{T}}_h}$$

$${⟨\zeta ,\rho ⟩}_{\partial {\mathcal{T}}_h}$$

$${\mathbf{V}}_h$$

$$W_h$$

$${\mathring{W}}_h$$

$$M_h$$

$$\mathbf{p}+\nabla \varphi$$

$$u_t+\nabla \cdot \mathbf{p}$$

$$\mathbf{q}+\nabla u$$

$$ϵ\nabla \cdot \mathbf{q}+{ϵ}^{-1}f(u)$$

$$\mathbf{p}\cdot \mathbf{n}=\mathbf{q}\cdot \mathbf{n}$$

$$u(\cdot ,0)$$

$$N$$

$$({\mathbf{p}}_h^n,{\mathbf{r}}_1{)}_{{\mathcal{T}}_h}-({\varphi }_h^n,\nabla \cdot {\mathbf{r}}_1{)}_{{\mathcal{T}}_h}+⟨{\hat{\varphi }}_h^n,{\mathbf{r}}_1\cdot \mathbf{n}{⟩}_{\partial {\mathcal{T}}_h}=0,$$

$$({\partial }_t^{+}{u}_h^n,w_1{)}_{{\mathcal{T}}_h}-({\mathbf{p}}_h^n,\nabla w_1{)}_{{\mathcal{T}}_h}+⟨{\hat{\mathbf{p}}}_h^n\cdot \mathbf{n},w_1{⟩}_{\partial {\mathcal{T}}_h}=0,$$

$$⟨{\hat{\mathbf{p}}}_h^n\cdot \mathbf{n},{\mu }_1{⟩}_{\partial {\mathcal{T}}_h}=0,$$

$$({\mathbf{q}}_h^n,{\mathbf{r}}_2{)}_{{\mathcal{T}}_h}-(u_h^n,\nabla \cdot {\mathbf{r}}_2{)}_{{\mathcal{T}}_h}+⟨{\hat{u}}_h^n,{\mathbf{r}}_2\cdot \mathbf{n}{⟩}_{\partial {\mathcal{T}}_h}=0,$$

$$-(ϵ{\mathbf{q}}_h^n,\nabla w_2{)}_{{\mathcal{T}}_h}+⟨ϵ{\hat{\mathbf{q}}}_h^n\cdot \mathbf{n},w_2{⟩}_{\partial {\mathcal{T}}_h}+({ϵ}^{-1}{f}^n(u_h^n),w_2{)}_{{\mathcal{T}}_h}=0,$$

$$⟨ϵ{\hat{\mathbf{q}}}_h^n\cdot \mathbf{n},{\mu }_2{⟩}_{\partial {\mathcal{T}}_h}=0,$$

$${\hat{\mathbf{p}}}_h^n\cdot \mathbf{n}={\mathbf{p}}_h^n\cdot \mathbf{n}+h_K^{-1}({\Pi }_k^{\partial }{\varphi }_h^n-{\hat{\varphi }}_h^n),$$

$${\hat{\mathbf{q}}}_h^n\cdot \mathbf{n}={\mathbf{q}}_h^n\cdot \mathbf{n}+h_K^{-1}({\Pi }_k^{\partial }{u}_h^n-{\hat{u}}_h^n),$$

$$⟨{\Pi }_k^{\partial }{u}_h,{\mu }_h{⟩}_E=⟨u_h,{\mu }_h{⟩}_E,\forall {\mu }_h\in {\mathcal{P}}^k(E)\text{and}E\in \partial K.$$

$$\begin{array}{rl}& \mathcal{A}({\mathbf{q}}_h,u_h,{\hat{u}}_h;{\mathbf{r}}_h,w_h,{\mu }_h)\\ & =({\mathbf{q}}_h,{\mathbf{r}}_h{)}_{{\mathcal{T}}_h}-(u_h,\nabla \cdot {\mathbf{r}}_h{)}_{{\mathcal{T}}_h}+⟨{\hat{u}}_h,{\mathbf{r}}_h\cdot \mathbf{n}{⟩}_{\partial {\mathcal{T}}_h}\\ & +(\nabla \cdot {\mathbf{q}}_h,w_h{)}_{{\mathcal{T}}_h}-⟨{\mathbf{q}}_h\cdot \mathbf{n},{\mu }_h{⟩}_{\partial {\mathcal{T}}_h}\\ & +⟨h_K^{-1}({\Pi }_k^{\partial }{u}_h-{\hat{u}}_h),{\Pi }_k^{\partial }{w}_h-{\mu }_h{⟩}_{\partial {\mathcal{T}}_h},\end{array}$$

$$({\partial }_t^{+}{u}_h^n,w_1{)}_{{\mathcal{T}}_h}+\mathcal{A}({\mathbf{p}}_h^n,{\varphi }_h^n,{\hat{\varphi }}_h^n;{\mathbf{r}}_1,w_1,{\mu }_1)$$

$$({ϵ}^{-1}{f}^n(u_h^n),w_2{)}_{{\mathcal{T}}_h}+ϵ\mathcal{A}({\mathbf{q}}_h^n,u_h^n,{\hat{u}}_h^n;{\mathbf{r}}_2,w_2,{\mu }_2)-({\varphi }_h^n,w_2{)}_{{\mathcal{T}}_h}$$

$$(u_h^0,w_3{)}_{{\mathcal{T}}_h}-(u^0,w_3{)}_{{\mathcal{T}}_h}$$

$$\begin{array}{rl}({\mathbf{\Pi }}_k^o\mathbf{q},{\mathbf{r}}_h{)}_K& =(\mathbf{q},{\mathbf{r}}_h{)}_K,\forall {\mathbf{r}}_h\in [{\mathcal{P}}_k(K){]}^d,\\ ({\Pi }_{k+1}^{o}u,w_h{)}_K& =(u,w_h{)}_K,\forall w_h\in {\mathcal{P}}_{k+1}(K),\end{array}$$

$$\Vert \mathbf{q}-{\mathbf{\Pi }}_k^o\mathbf{q}{\Vert }_{{\mathcal{T}}_h}\le C{h}^{k+1}\Vert \mathbf{q}{\Vert }_{H^{k+1}},\Vert u-{\Pi }_{k+1}^{o}u{\Vert }_{{\mathcal{T}}_h}\le C{h}^{k+2}\Vert u{\Vert }_{H^{k+2}},$$

$$\Vert u-{\Pi }_{k+1}^{o}u{\Vert }_{\partial {\mathcal{T}}_h}\le C{h}^{k+\frac{3}{2}}\Vert u{\Vert }_{H^{k+2}},\Vert w{\Vert }_{\partial {\mathcal{T}}_h}\le C{h}^{-\frac{1}{2}}\Vert w{\Vert }_{{\mathcal{T}}_h},\forall w\in V_h,$$

$$\Vert u-{\Pi }_{k+1}^{o}u{\Vert }_{L^{\infty }}\le C{h}^{k+2-d/2}|u{|}_{H^{k+2}}.$$

$$\displaystyle\|v\|_{\mathcal{T}_{h}}^{2}\leq C\left(\|\nabla v\|^{2}_{\mathcal{T}_{h}}+|(v,1)_{\mathcal{T}_{h}}|^{2}+\sum_{E\in\mathcal{E}^{o}_{h}}|E|^{d/(1-d)}\Big{(}\int_{E}[\![v]\!]ds\Big{)}^{2}\right),$$

$$\Vert v_h{\Vert }_{{\mathcal{T}}_h}^2\le C\left(\Vert \nabla v_h{\Vert }_{{\mathcal{T}}_h}^2+|(v_h,1{)}_{{\mathcal{T}}_h}{|}^2+\Vert h_K^{-1/2}({\Pi }_k^{\partial }{v}_h-{\hat{v}}_h){\Vert }_{\partial {\mathcal{T}}_h}^2\right).$$

$$\mathcal{A}({\mathbf{r}}_h,w_h,{\mu }_h,{\mathbf{r}}_h,w_h,{\mu }_h)=\Vert {\mathbf{r}}_h{\Vert }_{{\mathcal{T}}_h}^2+\Vert h_K^{-1/2}({\Pi }_k^{\partial }{w}_h-{\mu }_h){\Vert }_{\partial {\mathcal{T}}_h}^2.$$

$$\mathcal{A}({\mathbf{q}}_h,u_h,{\hat{u}}_h;{\mathbf{p}}_h,-{\varphi }_h,-{\hat{\varphi }}_h)=\mathcal{A}({\mathbf{p}}_h,{\varphi }_h,{\hat{\varphi }}_h;{\mathbf{q}}_h,-u_h,-{\hat{u}}_h).$$

$$\left|\mathcal{A}({\mathbf{q}}_h,u_h,{\hat{u}}_h;{\mathbf{p}}_h,{\varphi }_h,{\hat{\varphi }}_h)\right|\le C\left(\Vert {\mathbf{q}}_h{\Vert }_{{\mathcal{T}}_h}+\Vert \nabla u_h{\Vert }_{{\mathcal{T}}_h}+\Vert h_K^{-1/2}({\Pi }_k^{\partial }{u}_h-{\hat{u}}_h){\Vert }_{\partial {\mathcal{T}}_h}\right)$$

$$\times \left(\Vert {\mathbf{p}}_h{\Vert }_{{\mathcal{T}}_h}+\Vert \nabla {\varphi }_h{\Vert }_{{\mathcal{T}}_h}+\Vert h_K^{-1/2}({\Pi }_k^{\partial }{\varphi }_h-{\hat{\varphi }}_h){\Vert }_{\partial {\mathcal{T}}_h}\right).$$

$$\left|\mathcal{A}({\mathbf{q}}_h,u_h,{\hat{u}}_h;{\mathbf{p}}_h,{\varphi }_h,{\hat{\varphi }}_h)\right|$$

$$+\left|-({\mathbf{q}}_h,\nabla {\varphi }_h)+⟨{\mathbf{q}}_h\cdot \mathbf{n},{\Pi }_k^{\partial }{\varphi }_h-{\hat{\varphi }}_h{⟩}_{\partial {\mathcal{T}}_h}\right|$$

$$+\left|⟨h_K^{-1}({\Pi }_k^{\partial }{u}_h-{\hat{u}}_h),{\Pi }_k^{\partial }{\varphi }_h-{\hat{\varphi }}_h{⟩}_{\partial {\mathcal{T}}_h}\right|.$$

Full text

Keywords— Cahn-Hilliard equation, hybridizable discontinuous Galerkin method, superconvergence

1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) be a polygonal domain with Lipshitz boundary $\partial\Omega$ and $T$ be a positive constant. We consider the following Cahn-Hilliard equation:

[TABLE]

where $f(u)=u^{3}-u$. The Cahn-Hilliard equation is a fourth order, nonlinear parabolic equation which was originally proposed by Cahn and Hilliard [9, 8, 10] as a phenomenological model for phase separation and coarsening in a binary alloy. Since then Cahn-Hilliard-type equations have found applications in a variety of fields, including multiphase flow [46, 1], two-phase flow in porous media [36], tumor growth [64], pattern formation [65], thin films [6] and many others. Owing to its importance, many works have been devoted to the design and analysis of numerical schemes for solving the Cahn-Hilliard equation; see, e.g., finite difference methods [35], mixed and nonconforming finite element methods [29, 28, 30, 5, 34, 26] and Fourier-spectral methods [59, 44, 58].

In recent years, the discontinuous Galerkin (DG) method has become popular for solving the Cahn-Hilliard equation, owing to its flexibility in handling higher order derivatives, high-order accuracy, the property of local conservation which is crucial for applications in porous redmedium flow and transport phenomenon, high parallelizability and ease of achieving $hp$-adaptivity. Applications of DG methods to fourth order elliptic problems have been considered by Babuška and Zlámal in [3], by Baker in [4], and more recently by Mozolevski et al. in a series of works [48, 47, 50, 49, 62]. In [31], Feng and Karakashian design and analyze a DG method of interior penalty type based on the fourth order formulation of the Cahn-Hilliard equation. Optimal error estimates in various energy norms are established: see also [32]. Kay et al. propose and analyze a different DG method [42] that treats the Cahn-Hilliard equation as a system of second order equations allowing a relatively larger penalty term. A fully adaptive version of the interior penalty DG method was recently constructed in [2] for the Cahn-Hilliard equation with a source and optimal $L^{2}$ error bound were derived; see also [33] for solving the advective Cahn-Hilliard equation. The local discontinuous Galerkin (LDG) method has also been proposed for the discretization of the Cahn-Hilliard equation by writing it as a system of four first-order equations. Dong and Shu in [27] analyzed an LDG scheme for general elliptic equations including the linearized Cahn-Hilliard equation and obtained optimal error estimate in $L^{2}$. Recently, an LDG method has been employed for solving a number of Cahn-Hilliard fluid models, cf. [38, 39, 60].

The DG method is however often criticized for the larger amount of degrees of freedom compared to the continuous Galerkin (CG) method . In the seminal work [20] Cockburn et al. propose a hybridizable discontinuous Galerkin (HDG) method for second order elliptic problems. In a nutshell, the HDG method maps the flux and solution into the numerical trace of the solution via a local solver, which are in turn connected by the continuity of fluxes across inter-element boundaries (a transmission condition). Hence the globally coupled degrees of freedom are those numerical traces, resulting in a significant reduction of the number of unknowns in traditional DG methods. Moreover, the HDG methods possess the same favorable properties as classical mixed methods. In particular, HDG methods provide optimal convergence rates for both the gradient and the primal variables of the mixed formulation. This property enables the construction of superconvergent solutions, contrary to other DG methods. These advantages of the HDG methods have made HDG an attractive alternative for solving problems governed by PDEs and PDE control problems, see [24, 21, 25, 23, 11, 55, 54, 18, 56, 61, 13, 40, 37].

Most study currently focuses on establishing optimal and superconvergent rates of HDG methods for second order problems, such as elliptic PDEs [22], convection diffusion equations [16, 17, 52], Stokes equations [25, 21], Oseen equations [11] and Navier-Stokes equations [53, 12]. However, in [19], the authors utilized an HDG method with polynomial degree $k$ for all variables to investigate the biharmonic equations and obtained an optimal convergence rate for the solution and suboptimal convergence rates for the other variables. To the best of our knowledge, there does not exist an HDG work that achieves optimal convergence rates for all variables for a fourth order problem.

In this work, we propose a HDG method for the Cahn-Hilliard equation with Lehrenfeld-Schöberl stabilization function, polynomials of degree $k+1$ for the scalar unknown, and polynomials of degree $k$ for the other unknowns. The HDG framework with reduced stabilization and polynomials of mixed orders was first introduced by Lehrenfeld in [43] where it was alluded that the scheme could be a superconvergent method, i.e., $O(h^{k+2})$ error estimates is expected for the solution variables even though polynomials of degree $k$ are used for the globally coupled unknowns (numerical traces of the solution). Optimal convergence and hence superconvergence was then rigorously established for convection diffusion problems [52], Navier-Stokes equations [53], and more recently for linear elasticity problems [51]. We provide the HDG formulation for the Cahn-Hilliard equation in Section 2 and prove the existence, uniqueness and stability of the HDG method in Section 4. In Section 5, we perform a rigorous error analysis for the HDG method and obtain the following a priori error bounds for the solution $\phi$, $u$ and their fluxes $\bm{p}=-\nabla\phi$ and $\bm{q}=-\nabla u$:

[TABLE]

These convergence rates are further validated by numerical experiments in Section 6. A particular theoretical contribution of this article is the establishment of a novel HDG Sobolev inequality (cf. Theorem 3.14) which is a useful tool in the numerical analysis of nonlinear problems.

2 The HDG formulation

To introduce the fully discrete HDG formulation for the Cahn-Hilliard equation, we first fix some notation. Let $\mathcal{T}_{h}$ be a shape-regular, quasi-uniform triangulation of $\Omega$. Let $\mathcal{E}_{h}$ denote the set of all faces $E$ of all simplexes $K$ of the triangulation $\mathcal{T}_{h}$. Also let $\mathcal{E}_{h}^{o}$ and $\mathcal{E}_{h}^{\partial}$ denote the set of interior faces and boundary faces, respectively. Furthermore, we introduce the discrete inner products

[TABLE]

For any integer $k\geq 0$, let $\mathcal{P}^{k}(K)$ denote the set of polynomials of degree at most $k$ on the element $K$. We introduce the following discontinuous finite element spaces:

[TABLE]

where $L^{2}_{0}(\Omega)$ is the subspace of $L^{2}(\Omega)$ of mean zero functions.

Since the HDG method is based on a mixed formulation, we rewrite the Cahn-Hilliard equation as a first order system by setting $\bm{p}+\nabla\phi=0$ and $\bm{q}+\nabla u=0$ in (1). The mixed formulation of (1) then reads

[TABLE]

Now we introduce the fully discrete HDG formulation of the Cahn-Hilliard equation based on backward Euler method and convex-splitting approach.

[TABLE]

We shall also make use of the standard $L^{2}$ projection, denoted by $P_{M}$, onto $M_{h}$. Note that $\Pi_{k}^{\partial}$ coincides with $P_{M}$ on the space $H^{1}(\Omega)$.

To make the expressions concise, we introduce the operator $\mathcal{A}:[\bm{V}_{h}\times W_{h}\times M_{h}]^{2}\to\mathbb{R}$ by

[TABLE]

for all $(\bm{q}_{h},u_{h},\widehat{u}_{h}),(\bm{r}_{h},w_{h},\mu_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

Then the HDG formulation (3) can be recasted as: for $n=1,2,\cdots,N$, find $(\bm{p}_{h}^{n},\phi_{h}^{n},\widehat{\phi}_{h}^{n})$, $(\bm{q}_{h}^{n},u_{h}^{n},\widehat{u}_{h}^{n})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that

[TABLE]

3 Preliminaries

Throughout, $C$ shall denote a generic constant independent of the mesh parameters $h,\Delta t$. We first recall the standard $L^{2}$ projections $\bm{\Pi}_{k}^{o}:[L^{2}(\Omega)]^{d}\to\bm{V}_{h}$ and $\Pi_{k+1}^{o}:L^{2}(\Omega)\to W_{h}$

[TABLE]

which obey the following classical error estimates (see for instance [15, Lemma 3.3]):

[TABLE]

The same error bounds hold true for the projections of $\bm{p}$ and $\phi$.

We shall also utilize the following version of the piecewise Poincaré-Friedrichs inequality , cf. [7].

Lemma 3.1**.**

Let $v$ be a piecewise $H^{1}$ function with respect to the partition $\mathcal{T}_{h}$. The following Poincaré-Friedrichs inequality holds

[TABLE]

where the generic constant $C$ depends only on the regularity of the partition, and $[\![v]\!]$ denotes the jump of $v$ across a side $E$.

The following HDG Poincaré inequality is then an immediate consequence of Lemma 3.1, the Cauchy-Schwarz inequality and the triangle inequality.

Lemma 3.2** (The HDG Poincaré inequality).**

If $(v_{h},\widehat{v}_{h})\in W_{h}\times M_{h}$, we have

[TABLE]

Now, we glean some basic properties of the operator $\mathcal{A}$. First, the definition of $\mathcal{A}$ in Eq. (4) immediately implies lemmas 3.3 and 3.4.

Lemma 3.3**.**

For any $(\bm{r}_{h},w_{h},\mu_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$, we have

[TABLE]

Lemma 3.4**.**

For all $(\bm{q}_{h},u_{h},\widehat{u}_{h}),(\bm{p}_{h},\phi_{h},\widehat{\phi}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$, the operator $\mathcal{A}$ has the following property

[TABLE]

Next, we show that the operator $\mathcal{A}$ satisfies the following bound.

Lemma 3.5**.**

For all $(\bm{q}_{h},u_{h},\widehat{u}_{h}),(\bm{p}_{h},\phi_{h},\widehat{\phi}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$, we have

[TABLE]

Proof.

By the definition of $\mathcal{A}$ in (4) and integration by parts, one gets

[TABLE]

Then the bound (3.5) follows from the Cauchy-Schwarz inequality and the inverse inequality (7b). This completes the proof. ∎

Furthermore, we establish a crucial lemma that bounds the gradient of the scalar variable in terms of the flux variable and the reduced stabilization.

Lemma 3.6**.**

If $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ satisfies

[TABLE]

then the following inequality holds

[TABLE]

Proof.

By the definition of $\mathcal{A}$ in (4), let $\bm{r}_{h}=\nabla u_{h}$ in (12) and perform integration by parts to get

[TABLE]

Note that $\langle\Pi_{k}^{\partial}u_{h}-u_{h},\nabla u_{h}\cdot\bm{n}\rangle_{\partial\mathcal{T}_{h}}=0$. It follows from the element-wise Cauchy-Schwarz inequality and the inverse inequality (7b) that

[TABLE]

The triangle inequality gives

[TABLE]

The desired inequality (13) now follows by combining the last two inequalities. This completes the proof. ∎

Finally, we show that the operator $\mathcal{A}$ satisfies a version of the discrete LBB condition.

Lemma 3.7** (Discrete LBB Condition of $\mathcal{A}$).**

For all $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, we have

[TABLE]

Proof.

We only give the details of the proof of the inequality (14a), since the argument for (14b) is similar.

First, we note that $\|\bm{p}_{h}\|_{\mathcal{T}_{h}}+\|\nabla\phi_{h}\|_{\mathcal{T}_{h}}+\|h_{K}^{-1/2}(\Pi_{k}^{\partial}\phi_{h}-\widehat{\phi}_{h})\|_{\partial\mathcal{T}_{h}}$ defines a norm in the product space $\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, thanks to the HDG Poincaré inequality (8). Let $\alpha$ be a positive number to be specified later. For any fixed $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, we take $(\bm{p}_{h},\phi_{h},\widehat{\phi}_{h})=(\bm{q}_{h}+\alpha\nabla u_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$ to get

[TABLE]

By choosing $\alpha$ such that $1-C\alpha>0$, one obtains

[TABLE]

On the other hand, $(\bm{p}_{h},\phi_{h},\widehat{\phi}_{h})=(\bm{q}_{h}+\alpha\nabla u_{h},u_{h},\widehat{u}_{h})$ gives

[TABLE]

Then (14a) follows immediately. This completes the proof.

∎

We now introduce the HDG inversion of the Laplace operator equipped with homogeneous Neumann boundary condition.

Definition 3.8**.**

For any $u_{h}\in W_{h}$, we define $(\bm{\Pi}_{\bm{V}}u_{h},\Pi_{W}u_{h},\Pi_{M}u_{h})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$ such that

[TABLE]

for all $(\bm{r}_{h},w_{h},\mu_{h})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$,

Thanks to the discrete LBB condition Lemma 3.7, the inversion (15) in Definition 3.8 is well defined. For all $u_{h}\in\mathring{W}_{h}$, we define the semi-norm

[TABLE]

Then for all $u_{h}\in\mathring{W}_{h}$, by Lemma 3.3 and Definition 3.8, we have

[TABLE]

Next, we show that $\|\cdot\|_{-1,h}$ is a norm on the space $\mathring{W}_{h}$.

Lemma 3.9**.**

$\|\cdot\|_{-1,h}$* defines a norm on the space $\mathring{W}_{h}$.*

Proof.

Thanks to (16), one only needs to show that $\|u_{h}\|_{-1,h}=0$ implies $u_{h}=0$ for $u_{h}\in\mathring{W}_{h}$. It follows readily from (16) that

[TABLE]

Then Definition 3.8 and (4) give that for all $(\bm{r}_{h},w_{h})\in\bm{V}_{h}\times\mathring{W}_{h}$

[TABLE]

This is only possible if $u_{h}=0$. The proof is complete. ∎

For the negative norm $\|\cdot\|_{-1,h}$, the following HDG interpolation inequality holds true.

Lemma 3.10**.**

If $u_{h}\in\mathring{W}_{h}$ and $(w_{h},\mu_{h})\in W_{h}\times M_{h}$, one has

[TABLE]

where $h_{k}$ is the diameter of the element $K$.

Proof.

Let $(w_{h},\mu_{h})\in{W}_{h}\times M_{h}$ and $u_{h}\in\mathring{W}_{h}$. By Definition 3.8 and (4) we have

[TABLE]

By integration by parts, the identity (16) and the inverse inequality (7b) we have

[TABLE]

This concludes our proof. ∎

In addition, by the Definition 3.8, the identity (16) and Lemma 3.6 one can easily establish the following relation.

Lemma 3.11**.**

For any $u_{h}\in\mathring{W}_{h}$ there holds

[TABLE]

For the error analysis of the nonlinear equation we need to establish the discrete HDG Sobolev inequalities for which we will make use of the so-called Oswald interpolation operator [41].

Lemma 3.12** ([41]).**

There exists an interpolation operator, called Oswald interpolation, $\mathcal{I}_{h}^{c}:{W}_{h}\to{W}_{h}\cap H^{1}(\Omega)$, such that for any ${w}_{h}\in{W}_{h}$,

[TABLE]

where $[\![w_{h}]\!]$ denotes the jump of $w_{h}$ across a side $E$

Remark 3.13*.*

From the proof of [41, Page 644, Theorem 2.1], one can also obtain the following estimate:

[TABLE]

where $\varepsilon$ is any fixed real constant.

Now we are ready to prove the HDG Sobolev inequalities.

Theorem 3.14** (Discrete Sobolev inequalities).**

Let $\mu$ be the exponents as in the classical $H^{1}$ Sobolev embedding, i.e., $\mu$ satisfying

[TABLE]

For $w_{h}\in W_{h}$, it holds

[TABLE]

If further $w_{h}\in\mathring{W}_{h}$, then

[TABLE]

Proof.

We only give the proof of the inequality (23), since (24) is a direct consequence of the inequality (23) and the Poincaré-Friedrichs inequality in Lemma 3.1.

First, we note that the case $1\leq\mu<2$ is trivial since $L^{2}$ is embedded in $L^{\mu}$ by Hölder’s inequality. We consider the case $\mu\geq 2$.

By the triangle inequality we have

[TABLE]

Since $\mathcal{I}_{h}^{c}w_{h}\in H^{1}(\Omega)$, by the classical Sobolev embedding, the triangle inequality and Lemma 3.12, we have

[TABLE]

For the term $\|\mathcal{I}_{h}^{c}w_{h}-w_{h}\|_{0,\mu}$, we use the element-wise Sobolev embedding and the discrete Minkowski’s inequality to get

[TABLE]

where the last inequality follows from the fact $\mu\geq 2$ and the inequality $\sum_{i=1}^{n}|a_{i}|^{\mu/2}\leq(\sum_{i=1}^{n}|a_{i}|)^{\mu/2}$. Lemma 3.12 then yields

[TABLE]

The desired inequality (23) now follows from the inequalities (25) and (26). This completes the proof.

∎

The combination of the above theorem and the triangle inequality gives the following HDG Sobolev inequality.

Corollary 3.15** (HDG Sobolev inequality).**

For $(w_{h},\mu_{h})\in W_{h}\times M_{h}$, it holds

[TABLE]

if in addition $w_{h}\in\mathring{W}_{h}$, then

[TABLE]

where $\mu$ satisfying (22).

Proof.

For any $\mu_{h}\in M_{h}$, since $\mu_{h}$ is single-valued, by the triangle inequality, we have

[TABLE]

The HDG Sobolev inequality now follows from the inequality (23).

The second inequality (28) follows from the first inequality (27) and the HDG Poincaré inequality (8). This completes the proof. ∎

4 Well-posedness of the HDG formulation

In this section we establish the well-posedness of the HDG method (5), that is, existence and uniqueness of solutions as well as the energy stability of the solutions. The results differ slightly between the fully implicit discretization (FI) and the convex-spliting method (CS): the CS time marching enjoys unconditionally unique solvability and stability while there is a time-step constraint in the FI scheme for uniqueness and stability. For convenience, we will focus on the analysis of one method and point out the difference of the other.

4.1 Existence and uniqueness

Theorem 4.1**.**

The HDG scheme (5) admits at least one solution.

Proof.

We take $(\bm{r}_{1},w_{1},\mu_{1})=(\bm{0},1,1)$ and $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{0},1,1)$ in (5) to get

[TABLE]

Introduce the space

[TABLE]

with the inner product

[TABLE]

Thanks to Lemma 3.6 and the HDG Sobolev inequality (3.15), the inner product is well-defined on space $X$ with an induced norm $\|\cdot\|_{X}^{2}=(\cdot,\cdot)_{X}$. Now we introduce a nonlinear operator $G:X\to X$ as follows

[TABLE]

where $\alpha=\frac{1}{|\Omega|}(u^{0},1)_{\mathcal{T}_{h}}$. By Lemmas 3.5 and 3.6 and the HDG Sobolev inequality (3.15), the operator $G$ is a continuous operator. Therefore,

[TABLE]

By Lemma 3.5, the Cauchy-Schwarz inequality and noting that

[TABLE]

Hence, we have

[TABLE]

Hence for the fixed $u_{h}^{n-1}$, if $\|(\overline{\bm{p}}_{h}^{n},\overline{\phi}_{h}^{n},\widehat{\overline{\phi}}_{h}^{n},\overline{\bm{q}}_{h}^{n},\overline{u}_{h}^{n},\widehat{\overline{u}}_{h}^{n})\|_{X}$ is large enough, we have

[TABLE]

Then [63, II, Lemma 1.4] implies the existence of $(\overline{\bm{p}}_{h}^{n},\overline{\phi}_{h}^{n},\widehat{\overline{\phi}}_{h}^{n},\overline{\bm{q}}_{h}^{n},\overline{u}_{h}^{n},\widehat{\overline{u}}_{h}^{n})$ such that $G(\overline{\bm{p}}_{h}^{n},\overline{\phi}_{h}^{n},\widehat{\overline{\phi}}_{h}^{n},\overline{\bm{q}}_{h}^{n},\overline{u}_{h}^{n},\widehat{\overline{u}}_{h}^{n})=\bm{0}$.

Now we define

[TABLE]

with

[TABLE]

We proceed to show that $(\bm{p}^{n}_{h},\phi^{n}_{h},\widehat{\phi}^{n}_{h},\bm{q}_{h}^{n},u_{h}^{n},\widehat{u}_{h}^{n})$ is the solution to the HDG scheme (5).

Since

[TABLE]

and $(1,w_{1})_{\mathcal{T}_{h}}=(1,w_{2})_{\mathcal{T}_{h}}=0$, one gets that

[TABLE]

for all $(\bm{r}_{h},w_{1},\mu_{1},\bm{r}_{2},w_{2},\mu_{2})\in[\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}]^{2}$.

Noting that Eqs. (4.1) only hold for $w_{1},w_{2}\in\mathring{W_{h}}$. Next, we prove that they are true for the $w_{1},w_{2}\in W_{h}$. A direct calculation gives

[TABLE]

Likewise,

[TABLE]

Collectively, (4.1), (32), (4.1) implies

[TABLE]

holds for all $(\bm{r}_{h},w_{1},\mu_{1},\bm{r}_{2},w_{2},\mu_{2})\in[\bm{V}_{h}\times W_{h}\times M_{h}]^{2}$. This completes the proof. ∎

Next we show that the solution to the fully discrete scheme with the convex-splitting $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n-1}$ is unique, while the solution is only conditionally unique for the case of the Backward Euler temporal discretization.

Theorem 4.2**.**

The solution to the HDG scheme of (5) with the splitting $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n-1}$ is unique. On the other hand, the solution corresponding to the Backward Euler discretization is unique provided that $\Delta t\leq C\epsilon^{3}$.

Proof.

Given $(\bm{p}_{h}^{n-1},\phi_{h}^{n-1},\widehat{\phi}_{h}^{n-1},\bm{q}_{h}^{n-1},u_{h}^{n-1},\widehat{u}_{h}^{n-1})$, we let $(\bm{p}_{h,1}^{n},\phi_{h,1}^{n},\widehat{\phi}_{h,1}^{n},\bm{q}_{h,1}^{n},\\ u_{h,1}^{n},\widehat{u}_{h,1}^{n})$ and $(\bm{p}_{h,2}^{n},\phi_{htherestof,2}^{n},\widehat{\phi}_{h,2}^{n},\bm{q}_{h,2}^{n},u_{h,2}^{n},\widehat{u}_{h,2}^{n})$ be two solutions of (5). Let

[TABLE]

After inserting the two solutions into (5), subtracting the two equations gives

[TABLE]

Take $(\bm{r}_{1},w_{1},\mu_{1},\bm{r}_{2},w_{2},\mu_{2})=(\Delta t\bm{P}_{h}^{n},\Delta t\Phi_{h}^{n},\Delta t\widehat{\Phi}_{h}^{n},\bm{Q}_{h}^{n},U_{h}^{n},\widehat{U}_{h}^{n})$ to get

[TABLE]

We add the above two equations together to get

[TABLE]

In the case of the convex-splitting $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n-1}$, one has

[TABLE]

Hence all terms on the left hand side of (4.1) are nonnegative. It follows that

[TABLE]

Now, we take $(\bm{r}_{1},w_{1},\mu_{1})=(\bm{0},1,1)$ in (34a) and $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{0},1,1)$ in (34b) to get $(U_{h}^{n},1)_{\mathcal{T}_{h}}=0$, $(\Phi_{h}^{n},1)_{\mathcal{T}_{h}}=0$. By (28), Lemma 3.6 and (36), we have $(\bm{P}_{h}^{n},\Phi_{h}^{n},\widehat{\Phi}_{h}^{n},\bm{Q}_{h}^{n},U_{h}^{n},\widehat{U}_{h}^{n})=0$.

For the case of the Backward Euler method, i.e., In the case of the convex-splitting $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n}$, one has

[TABLE]

In light of Eq. (34a) we estimate $\epsilon^{-1}\|U_{h}^{n}\|_{\mathcal{T}_{h}}^{2}$ using the continuity of $\mathcal{A}$ and Lemma 3.6 as follows

[TABLE]

Hence provided that $\Delta t\leq\frac{\epsilon^{3}}{2C}$, Eq. (4.1) yields

[TABLE]

One then obtains uniqueness of the solution for the Backward Euler scheme. This completes the proof. ∎

4.2 Energy stability

In this subsection, we analyze the stability of the HDG formulation (5), focusing on the fully implicit scheme, i.e., $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n}$. We first recall some useful identities.

Lemma 4.3**.**

Let a,b be two real numbers. Let $\{a_{n}\}_{n=0}^{m}$ and $\{b_{n}\}_{n=0}^{m}$ be two sequences such that $b_{0}=0$. Then the following identities hold

[TABLE]

The first energy identity makes use of the negative norm and takes the following form.

Lemma 4.4** (Discrete energy identy I).**

Let $(\bm{p}_{h}^{n},\phi_{h}^{n},\widehat{\phi}_{h}^{n})$, $(\bm{q}_{h}^{n},u_{h}^{n},\widehat{u}_{h}^{n})$ be the solution of (5). The following energy identity holds for $m=1,\ldots,N$

[TABLE]

Proof.

Due to the mass conservation property of the scheme, it holds $\partial_{t}^{+}u_{h}^{n}\in\mathring{W}_{h}$. We take $(\bm{r}_{1},w_{1},\mu_{1})=-(\bm{\Pi}_{\bm{V}}\partial_{t}^{+}u_{h}^{n},-\Pi_{W}\partial_{t}^{+}u_{h}^{n},-\Pi_{M}\partial_{t}^{+}u_{h}^{n})$ in (5a) to get

[TABLE]

In view of (16), one obtains that

[TABLE]

Next, we take $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{0},\partial_{t}^{+}u_{h}^{n},\partial_{t}^{+}\widehat{u}_{h}^{n})$ in (5b) to get

[TABLE]

Then we apply $\partial_{t}^{+}$ to (5b) and take $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{q}_{h}^{n},0,0)$ so that

[TABLE]

Taking summation of (39), (40) and (41) gives

[TABLE]

Then multiply this equation by $\Delta t$ and use the identities in Lemma 4.3 to get

[TABLE]

The identity (38) follows from summing the above equation from $n=1$ to $n=m$. This completes the proof. ∎

Next, we give the second energy identity which involves the $L^{2}$ norm of $\bm{p}_{h}^{n}$.

Lemma 4.5** (Discrete energy identity II).**

Let $(\bm{p}_{h}^{n},\phi_{h}^{n},\widehat{\phi}_{h}^{n})$, $(\bm{q}_{h}^{n},u_{h}^{n},\widehat{u}_{h}^{n})$ be the solution of (5). Then for $m=1,\ldots N$, we have the following energy identity

[TABLE]

Proof.

We take $(\bm{r}_{1},w_{1},\mu_{1})=(\bm{p}_{h}^{n},\phi_{h}^{n},\widehat{\phi}_{h}^{n})$ in (5a) to get

[TABLE]

Taking the sum of (40), (41) and (43) gives

[TABLE]

The energy identity now follows from applications of the identities in Lemma 4.3 and addition of the resulting equation from $n=1$ to $n=m$. ∎

There is a negative term, i.e., $-\frac{(\Delta t)^{2}}{2\epsilon}\|\partial_{t}^{+}u_{h}^{n}\|^{2}_{\mathcal{T}_{h}}$ in the energy identities I and II in the foregoing Lemmas 4.4 and 4.5. In Theorem 4.6, we bound this term and derive the stability result.

Theorem 4.6**.**

Let $C$ be the product of the constants in (17) and Lemma 3.6. If the time-step satisfies the constraint $\Delta t\leq\frac{2\epsilon^{3}}{C}$, then for $m=1,2,\ldots N$, the following energy bounds hold

[TABLE]

Proof.

Since $\partial_{t}^{+}u_{h}^{n}\in\mathring{W}_{h}$, by the interpolation inequality (17) and Lemma 3.6 we have

[TABLE]

Now, by the energy identity I (38) and the assumption $\Delta t\leq\frac{2\epsilon^{3}}{C}$, we have

[TABLE]

The energy bounds (44) now follow from the energy law above, the energy identity II (42) and Lemma 3.6. ∎

Remark 4.7*.*

In the case of the energy-splitting scheme $f^{n}(u_{h}^{n})=(u_{h}^{n})^{3}-u_{h}^{n-1}$, by the identity (37c), an energy law holds in the spirit of Lemma 4.4 where all terms are associated with the positive sign.

We note that in Theorem 4.6, the energy term $\|\phi_{h}^{m}\|_{\mathcal{T}_{h}}^{2}$ is not contained. Moreover, the HDG Poincaré inequality (28) does not apply to $\phi_{h}^{n}$ since $\phi_{h}^{n}\notin\mathring{W}_{h}$. Hence, we need a refined analysis for this term. In the following we derive further a priori bounds for the solution of the fully implicit HDG scheme (5) with the assumption $\Delta t\leq\frac{2\epsilon^{3}}{C}$.

Lemma 4.8**.**

The assumption is the same as in Theorem 4.6. There holds

[TABLE]

where $C$ may depend on $\epsilon,T$ and the initial condition.

Proof.

Taking $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{0},1,1)$ in Eq. (5b) yields

[TABLE]

where the last inequality follows from the stability bounds in (44). In light of the stability bounds from (44) and the HDG Poincaré inequality (8), one gets that

[TABLE]

We apply $\partial_{t}^{+}$ to (5b) and keep (5a) unchanged to get

[TABLE]

Take $(\bm{r}_{1},w_{1},\mu_{1})=-\epsilon(\partial_{t}^{+}\bm{q}_{h}^{n},-\partial_{t}^{+}u_{h}^{n},-\partial_{t}^{+}\widehat{u}_{h}^{n})$ in (47a) and $(\bm{r}_{2},w_{2},\mu_{2})=(\bm{p}_{h}^{n},\\ -\phi_{h}^{n},-\widehat{\phi}_{h}^{n})$ in (47b) to get

[TABLE]

In light of the symmetry property (10), adding the above two equations together gives

[TABLE]

By the identity (37a) we have

[TABLE]

and hence

[TABLE]

By Hölder’s inequality, we have

[TABLE]

Next, by the HDG Sobolev inequality in Corollary 3.15, we have

[TABLE]

We then use Theorem 4.6 and the Cauchy-Schwarz inequality to get

[TABLE]

Together with (49) and the above inequality, we have

[TABLE]

From the HDG Sobolev embedding inequality Corollary 3.15 again, one gets

[TABLE]

where one uses the estimate (46) and the bounds in Theorem 4.6 in the derivation of the last step. This finishes the proof.

∎

4.3 The uniform estimate of $u_{h}^{n}$

We now bound $u_{h}^{n}$ in the $L^{\infty}$ norm. For this, we introduce the discrete Laplacian operator. For any $u_{h}\in W_{h}$, we define

$(\bm{q}_{h},\Delta_{h}u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that

[TABLE]

for all $(\bm{r}_{h},w_{h},\mu_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$. One can verify that $\Delta_{h}u_{h}$ is well-defined by the linearity of the operator $\mathcal{A}$ and the fact that uniqueness implies existence for a linear square system of finite dimension, cf. (4).

Lemma 4.9**.**

Let $u_{h}^{n}$ be the solution of (5). For all $n=1,2\ldots,N$, we have

[TABLE]

where $C$ depends on $\epsilon,T$ and the initial condition.

Proof.

For $u_{h}=u_{h}^{n}$ in Eq. (50), by the uniqueness of the solution, in light of Eq. (5b), one identifies that $\bm{q}_{h}=\bm{q}_{h}^{n}$ and $\widehat{u}_{h}=\widehat{u}_{h}^{n}$. Now taking $w_{2}=\Delta_{h}u_{h}^{n}$ in (5b), one obtains

[TABLE]

It follows from the Cauchy-Schwarz inequality

[TABLE]

Now, one estimates the term $\|u_{h}\|_{0,6}$ by the HDG Sobolev inequality (27) and the stability bounds in Theorem 4.6 as follows

[TABLE]

One readily obtains the estimate (51) in light of the stability bound(45). This completes the proof. ∎

We now estimate the $L^{\infty}$ norm of functions in $W_{h}$.

Lemma 4.10**.**

For all $w_{h}\in W_{h}$, we have the inequality

[TABLE]

Proof.

Consider the following continuous problem: find $w\in H^{1}(\Omega)$ such that

[TABLE]

Since $\Omega$ is convex, we have the regularity estimate

[TABLE]

where $\overline{w}=\frac{1}{|\Omega|}(w_{h},1)_{\mathcal{T}_{h}}$. The definitions (50) and (52) imply

[TABLE]

for all $(\bm{s}_{h},v_{h},\mu_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$. By the uniqueness of solutions to the elliptic projection (59), in view of Remark 5.4, one uses the HDG elliptic projection result (70) to get

[TABLE]

By the triangle inequality, we have

[TABLE]

The $\{R_{i}\}_{i=1}^{3}$ are estimated as follows

[TABLE]

where the last inequality follows from the elliptic regularity result. Collecting the above estimates, one concludes the proof.

∎

Using Lemma 4.10 and Lemma 4.9 immediately gives the following result.

Lemma 4.11**.**

Let $u_{h}^{n}$ be the solution of (5). For all $n=1,2\ldots,N$, we have

[TABLE]

where $C$ depends on $\epsilon,T$ and the initial condition.

5 Error analysis

In this section, we provide a convergence analysis of the fully implicit HDG method for the Cahn-Hilliard equation. The convex-splitting scheme can be similarly treated. First, we give our main results. Then, we define an HDG elliptic projection as in [14], which is a crucial step to prove the main result. In the end, we provide rigorous error estimation for our fully implicit HDG method.

Throughout, we assume the data and the solution of (1) are smooth enough. As in Section 4, we do not track the dependence on $\epsilon$ and treat as if $\epsilon=O(1)$. The generic constant $C$ may depend on the data of the problem but is independent of $h$ and may change from line to line.

Given $\Theta\in L^{2}_{0}(\Omega)$, let $(\bm{\Psi},\Phi)$ be the solution of the following system

[TABLE]

If $\Omega$ is convex, then we have the following regularity result:

[TABLE]

5.1 The main result

We can now state our main result for the HDG method.

Theorem 5.1**.**

Let $(\bm{p},\phi,\bm{q},u)$ and $(\bm{p}_{h}^{n},\phi_{h}^{n},\bm{q}_{h}^{n},u_{h}^{n})$ be the solutions of (2) and (5), respectively. Assume the solution $(\bm{p},\phi,\bm{q},u)$ attains the maximum regularity for the best approximation results in (7). If $\Delta t\leq C\epsilon^{3}$ for the BE scheme and is arbitrary for the CS scheme, one has the following optimal error estimates

[TABLE]

Furthermore, if the polynomial order $k\geq 1$, one also has the optimal error estimate in the negative norm

[TABLE]

Remark 5.2*.*

To the best of our knowledge, [19] is the only work for fourth order problems using an HDG method with polynomial degree $k$ for all variables. They obtained an optimal convergence rate for the solution and suboptimal convergence rates for the other variables. In contrast, the HDG method proposed in this work deals with a nonlinear fourth order problem and achieves optimal convergence rates for all variables. Moreover, from the view point of degrees of freedom, we obtain the superconvergent rate for the solution.

5.2 The HDG elliptic projection

For all $t\in[0,T]$, we define the HDG elliptic projection: find $(\bm{p}_{Ih},\phi_{Ih},\widehat{\phi}_{Ih}),\\ (\bm{q}_{Ih},u_{Ih},\widehat{u}_{Ih})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that

[TABLE]

for all $(\bm{r}_{1},w_{1},\mu_{1}),(\bm{r}_{2},w_{2},\mu_{2})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$. These projections are well-defined in the sense that there exist unique $(\bm{p}_{Ih},\phi_{Ih},\widehat{\phi}_{Ih}),(\bm{q}_{Ih},u_{Ih},\widehat{u}_{Ih})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that Eqs. (59) hold.

We have the following approximation property for the HDG elliptic projection (59).

Theorem 5.3**.**

Let $(\bm{p},\phi,\bm{q},u)$ and $(\bm{p}_{Ih},\phi_{Ih},\bm{q}_{Ih},u_{Ih})$ be the solution of (2) and (59), respectively. For all integer $s\in[0,k],k\geq 0$ we have

[TABLE]

Remark 5.4*.*

In the proof of Theorem 5.3 below, we only make use of the regularity of $(\bm{p},\phi,\bm{q},u)$ and the fact that $\bm{p}=-\nabla\phi,\quad\bm{q}=-\nabla u$. Hence the approximation properties in Theorem 5.3 are valid for any regular functions with $\bm{p}=-\nabla\phi,\quad\bm{q}=-\nabla u$.

We only give a proof of (60a) and (60b), and we split the proof into three steps. To simplify the notation, we define

[TABLE]

Note that $\varepsilon_{h}^{\widehat{u}}\in M_{h}$ since $\Pi_{k}^{\partial}u=P_{M}u$ for $u\in H^{1}(\Omega)$, where $P_{M}$ is the $L^{2}$ orthogonal projection onto $M_{h}$. The following error estimate of $\Pi_{k}^{\partial}$ ($P_{M}$) is classical

[TABLE]

5.2.1 Step 1: The error equation

Lemma 5.5**.**

For all $(\bm{r}_{2},w_{2},\mu_{2})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, we have

[TABLE]

Proof.

By the definition of $\mathcal{A}$ in (4), we have

[TABLE]

where we used the orthogonality of $\bm{\Pi}_{k}^{o}$, $\Pi_{k+1}^{o}$, $\Pi_{k}^{\partial}$ in the last equality. Since $\bm{q}=-\nabla u$ and $\langle\bm{q}\cdot\bm{n},\mu_{2}\rangle_{\partial\mathcal{T}_{h}}=0$, one gets

[TABLE]

This completes the proof. ∎

Subtracting Eq. (63) and Eq. (59b) gives the error equation.

Lemma 5.6**.**

For all $(\bm{r}_{2},w_{2},\mu_{2})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, we have

[TABLE]

5.2.2 Step 2: An energy argument

Lemma 5.7**.**

Let $(\bm{q},u)$ and $(\bm{q}_{Ih},u_{Ih})$ be the solution of (1) and (59b), respectively. The following error estimate holds for $s\in[0,k]$ and $k\geq 0$.

[TABLE]

In particular, one has

[TABLE]

Proof.

First, the error equation (64) implies that

[TABLE]

Hence Lemma 3.6 gives

[TABLE]

Noting that $\varepsilon_{h}^{u}\in\mathring{W}_{h}$ by the definitions of $\Pi_{k+1}^{o}$ and $u_{Ih}$, we now take $(\bm{r}_{2},w_{2},\mu_{2})=(\varepsilon_{h}^{\bm{q}},\varepsilon_{h}^{u},\varepsilon_{h}^{\widehat{u}})$ in (64). Then by the Cauchy-Schwartz inequality, the triangle inequality, the inequality (68) and inequalities in (7), one get

[TABLE]

The error estimate (65) readily follows. Then the estimate (66) is a consequence of the inequality (68).

Now in light of the definitions of the error functions in (61), one obtains the desired error estimate (67) by the triangle inequality, the $L^{2}$ stability of the projection $\Pi_{k}^{\partial}$, the inequalities in (7) and the fact that $\bm{q}=-\nabla u$. This completes the proof. ∎

5.2.3 Step 3: The error estimate of the scalar variable by the duality argument

Similar to Lemma 5.5 we have the following result.

Lemma 5.8**.**

Let $\Theta$ be in $\mathring{W}_{h}$, and let $(\bm{\Psi},\phi)$ be the solution to the system (55). Then for all $(\bm{r}_{2},w_{2},\mu_{2})\in\bm{V}_{h}\times\mathring{W}_{h}\times M_{h}$, we have the equation

[TABLE]

Lemma 5.9**.**

Let $u$ and $u_{Ih}$ be the solutions of (1) and (59b), respectively. Then for $s\in[0,k]$, $k\geq 0$ we have the error estimates

[TABLE]

Proof.

We take $(\bm{r}_{2},w_{2},\mu_{2})=(\varepsilon_{h}^{\bm{q}},-\varepsilon_{h}^{u},-\varepsilon_{h}^{\widehat{u}})$ and $\Theta=-\varepsilon_{h}^{u}$ in (69) to get

[TABLE]

By Lemma 3.4, the error equation (64), we have

[TABLE]

Since $\Pi_{k}^{\partial}\Phi$ is single valued, one has $\langle\bm{q}\cdot\bm{n},\Pi_{k}^{\partial}\Phi\rangle_{\partial\mathcal{T}_{h}}=\langle\bm{q}\cdot\bm{n},\Phi\rangle_{\partial\mathcal{T}_{h}}=0$. Hence

[TABLE]

It follows that

[TABLE]

where the regularity result (56) with $\Theta=-\varepsilon_{h}^{u}$ is used in the derivation of the last inequality. The rest of the terms can be dealt with similarly by using the $L^{2}$ stability of the projection $\Pi_{k}^{\partial}$, and the inequalities in (7).

The desired error estimate (70) now follows from the error estimates (65) and (66), and the triangle inequality. This completes the proof.

∎

5.3 Error estimate in the negative norm

To establish the approximation properties of the elliptic projection in the negative norm, we introduce a Scott-Zhang type (cf. [57]) interpolation operator \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+2+d} in Section 7. For all $(u_{h},\widehat{u}_{h})\in L^{2}(\mathcal{T}_{h})\times L^{2}(\partial\mathcal{T}_{h})$, \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+2+d}(u_{h},\widehat{u}_{h})|_{K}\in\mathcal{P}^{\color[rgb]{0,0,0}k+2+d}(K) and satisfies

[TABLE]

Theorem 5.10**.**

Let $(\phi,u)$ and $(\phi_{Ih},u_{Ih})$ be the solution of (1) and (59a)-(59b), respectively. Then if $k\geq 1$, we have the following error estimates

[TABLE]

Proof.

We only give a proof for (72a), the proof for (72b) is similar.

Let $\xi_{h}^{u}=\Pi_{k+1}^{o}u-u_{Ih}$, by Definition 3.8 and (71a) one gets

[TABLE]

We take $(\bm{r}_{2},w_{2},\mu_{2})=(\xi_{h}^{\bm{q}},-\xi_{h}^{u},-\xi_{h}^{\widehat{u}})$ and \Theta=-\mathcal{I}_{h}^{\color[rgb]{0,0,0}k+2+d}(\Pi_{W}\xi_{h}^{u},\Pi_{M}\xi_{h}^{u}) in (69), and use (73) to get

[TABLE]

Hence

[TABLE]

Since $\Pi_{W}\xi_{h}^{u}\in\mathring{W}_{h}$, and by the $H^{1}$ stability of the interpolation operator (120) one has

[TABLE]

which combining (74), then implies

[TABLE]

This completes the proof. ∎

In a similar fashion as Lemma 4.11 one can establish the stability bound of $u_{Ih}$ in the uniform norm.

Lemma 5.11**.**

Let $u_{Ih}$ be the solution to the elliptic projection (59b). Assume $u\in L^{\infty}(0,T;H^{2}(\Omega))$. Then one has

[TABLE]

where $C$ depends on $\|u\|_{L^{\infty}(H^{2})}$.

5.4 Proof of Theorem 5.1

To simplify notation, we define

[TABLE]

Lemma 5.12**.**

For all $(\bm{r}_{1},w_{1},\mu_{1}),(\bm{r}_{2},w_{2},\mu_{2})\in\bm{V}_{h}\times W_{h}\times M_{h}$, we have the following error equations

[TABLE]

Proof.

We use the definition of $\mathcal{A}$ in (4) to get

[TABLE]

Next, we have

[TABLE]

This completes the proof. ∎

We start the error analysis in the negative norm. We have

Lemma 5.13** (Error estimates in the $-1$ norm).**

The following error estimates hold

[TABLE]

Proof.

Taking $(\bm{r}_{1},w_{1},\mu_{1})=(-\bm{\Pi_{V}}e_{h}^{u^{n}},\Pi_{W}e_{h}^{u^{n}},\Pi_{M}e_{h}^{u^{n}})$ in (76a) gives

[TABLE]

Lemma 3.4 and Definition 3.8 imply

[TABLE]

Hence

[TABLE]

Next, one takes $(\bm{r}_{2},w_{2},\mu_{2})=(e_{h}^{\bm{q}^{n}},e_{h}^{u^{n}},e_{h}^{\widehat{u}^{n}})$ in (76b) to get

[TABLE]

Combining Eqs. (82) and (83) one obtains

[TABLE]

We first calculate the first term on the left hand side of Eq. (84). Utilizing Definition 3.8 and (4), one has

[TABLE]

On the other hand,

[TABLE]

Hence taking $\bm{r}_{h}=\bm{\Pi_{V}}\partial_{t}^{+}e_{h}^{u^{n}}$, $w_{h}=\mu_{h}=0$ in (86) yields

[TABLE]

Eq. (5.4) now becomes

[TABLE]

where one has utilized Eq. (16). Substituting (87) into the error equation (84), one has

[TABLE]

The three terms on the right-hand side of Eq. (5.4) are estimated as follows. By Equation 17 and Lemma 3.6 one has

[TABLE]

with $\theta$ an arbitrary positive constant. Likewise, Equation 17 and Lemma 3.11 implies

[TABLE]

where $\kappa$ is another free parameter. For the nonlinear term one has

[TABLE]

It follows

[TABLE]

where $\xi_{h}^{n}:=(u_{h}^{n})^{2}+u_{h}^{n}u_{Ih}^{n}+(u_{Ih}^{n})^{2}\geq 0$, and one has utilized the element-wise duality of $H^{1}$, the uniform bound in Lemma 5.11, and the stability bounds in Lemma 4.8.

Taking $\theta=\frac{\epsilon}{4}$ and substituting the inequalities (5.4)–(5.4) back into (5.4), then multiplying the resulting equation by $\Delta t$ and taking summation from $n=1$ to $n=m$ gives

[TABLE]

Since $\kappa$ is an arbitrary positive number, one can choose the maximum of $\kappa$ such that $\kappa\Delta t\leq\frac{1}{4}$. An application of Gronwall’s inequality then gives the error estimate (79). This completes the proof. ∎

Next we derive the error estimates of the scalar variables in the $L^{2}$ norm.

Lemma 5.14**.**

The scalar variables satisfy the following error bounds

[TABLE]

Proof.

First, we take $(\bm{r}_{1},w_{1},\mu_{1})=(-e_{h}^{\bm{q}^{n}},e_{h}^{u^{n}},e_{h}^{\widehat{u}^{n}})$ in (76a) to get

[TABLE]

Next, we take $(\bm{r}_{2},w_{2},\mu_{2})=(e_{h}^{\bm{p}^{n}},-e_{h}^{\phi^{n}},-e_{h}^{\widehat{\phi}^{n}})$ in (76b) to get

[TABLE]

We multiply (94) by $\epsilon$ and then add (95) to get

[TABLE]

By the $\L^{\infty}$ stability bound of $u_{h}^{n}$ in Lemma 4.11, we have

[TABLE]

By Lemma 4.3 it follows that

[TABLE]

Apply Cauchy-Schwartz inequality to (96), and then add the resulting equation from $n=1$ to $n=m$ to get

[TABLE]

where one has applied the approximation properties of the elliptic projection in Theorem 5.3. The error estimate in the negative norm (79) implies that

[TABLE]

Hence one has

[TABLE]

This establishes the optimal error estimates of $u$ and $\phi$ in the $L^{2}$ norm.

∎

Finally one performs the error analysis of the flux variables. One has

Lemma 5.15**.**

The scalar variables satisfy the following error bounds

[TABLE]

Proof.

First we take $(\bm{r}_{1},w_{1},\mu_{1})=(e_{h}^{\bm{p}^{n}},e_{h}^{\phi^{n}},e_{h}^{\widehat{\phi}^{n}})$ in Eq. (76a) to get

[TABLE]

Applying $\partial_{t}^{+}$ to Eq. (76b) and then setting $(\bm{r}_{2},w_{2},\mu_{2})=(e_{h}^{\bm{q}^{n}},0,0)$ in the resulting equation gives

[TABLE]

One now takes $(\bm{r}_{2},w_{2},\mu_{2})=(0,\partial_{t}^{+}e_{h}^{u^{n}},\partial_{t}^{+}e_{h}^{\widehat{u}^{n}})$ in Eq. (76b), and take summation of the result with Eqs. (100)–(101) to obtain

[TABLE]

For $I_{1}$, by the Cauchy-Schwartz inequality, the HDG Sobolev inequality Corollary 3.15, Lemma 3.6 and the approximation properties in Theorem 5.3, one obtains

[TABLE]

with $\theta$ a positive number to be chosen later.

Denote by $P_{W}$ the $L^{2}$ projection operator onto $W_{h}$. By Eq. (76a) with $(\bm{r}_{1},\mu_{1})=(0,0)$, one has

[TABLE]

where one has utilizes the continuity of the operator $\mathcal{A}$ in Lemma 3.5, the inverse inequality, the $L^{2}$ stability of the projections $\Pi_{k}^{\partial}$ and $P_{W}$.

The term $I_{3}$ is estimated similarly as the term $I_{2}$ as follows.

[TABLE]

where the uniform bound of $u_{Ih}^{n}$ in Lemma 5.11 has been applied here.

For $I_{4}$ one has

[TABLE]

We estimate $I_{5}$ following the approach in [45]. Choosing $\theta=\frac{1}{6}$, substituting the inequalities (103)–(5.4) back to the Eq. (5.4), multiplying the result by $\Delta t$ and then taking summation from $n=1$ to $n=m$ implies

[TABLE]

where the last step follows from the $L^{2}$ error estimate in (98).

By the identity (37d) the last term in (5.4) can be written as

[TABLE]

with $\xi^{n}=(u_{Ih}^{n})^{2}+u_{Ih}^{n}u_{h}^{n}+(u_{h}^{n})^{2}$. Noting that $\xi^{n}\geq 0$, the last two terms in (5.4) are non-positive. One has

[TABLE]

Hence by the uniform stability of $u_{h}^{n}$ and $u_{Ih}^{n}$ in Lemma 4.11 and Lemma 5.11, one obtains

[TABLE]

where one has applied the HDG Sobolev inequality (28) and Lemma 3.6 in deriving the last inequality.

With (5.4) and (5.4), the inequality (5.4) becomes

[TABLE]

By Lemma 4.8 and Theorem 5.3, one has

[TABLE]

An application of the Gronwall’s inequality gives

[TABLE]

This completes the proof. ∎

The combination of (60a), (60d), (60e) and the triangle inequality finishes the proof of Theorem 5.1.

6 Numerical results

We consider two examples on unit square domains in $\mathbb{R}^{2}$. In the first example we have an explicit solution of the system (1); in the second example an explicit form for the exact solution is not known.

Example 6.1**.**

The problem data $u^{0}$ and the artificial $f$ are chosen so that the exact solution of the system (1) is given by

[TABLE]

We report the errors at the final time $T=1$ for polynomial degrees $k=0$ and $k=1$ in Tables 1 and 2 for the fully implicit scheme and Tables 3 and 4 for the energy-splitting scheme. The observed convergence rates match the theory, where $\Delta t=h^{k+1}$.

7 Appendix

Definition 7.1**.**

[TABLE]

(1) For every vertex $A_{i}$ on mesh $\mathcal{T}_{h}$, let $N_{A_{i}}$ be the number of elements adjoint at $A_{i}$, and $\mathcal{K}_{A_{i}}$ denote all these elements, then \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d} at $A_{i}$ is defined as

[TABLE]

we note that $N_{A_{i}}$ is a fixed finite number since $\mathcal{T}_{h}$ is shape-regular.

(2) If, in addition, for $d=3$, for every edge $E$ of element $K$, there are $k+d$ interior Lagrange points on edge $E$, for any of these points $B_{i}$, let $N_{B_{i}}$ be the number of elements adjoint at $B_{i}$, and $\mathcal{K}_{B_{i}}$ denote all these elements, then \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d} at $B_{i}$ is defined as

[TABLE]

Again, $N_{B_{i}}$ is finite since $\mathcal{T}_{h}$ is shape-regular.

(3) Since \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h})\in\mathcal{P}^{k+1+d}(K), there are $\binom{k+d}{d-1}$ Lagrange points on every face $F$ of $K$, the value of \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) on these points are determined by

[TABLE]

holds for all face $F$ of $K$.

(4) Since \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h})\in\mathcal{P}^{k+1+d}(K), there are $\binom{k+d}{d}$ Lagrange points in every element $K$, the value of \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) on these points are determined by

[TABLE]

It is easy to check that the degrees of freedom of $\mathcal{P}^{k+1+d}(K)$ is $\binom{k+1+2d}{d}$, the constrains for (111a), (111b), (111c) and (111d) are $d+1$, $(d-1)d\binom{k+d}{1}$, $(d+1)\binom{k+d}{d-1}$, and $\binom{k+d}{d}$. For $d=3$, there holds

[TABLE]

and for $d=2$, it holds

[TABLE]

Then the definition of \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) is a square system, therefore, the uniqueness and the existence of \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) are equivalence. In addition, it is obviously that when $u_{h}=\widehat{u}_{h}=0$ we have \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h})=0, then the operator \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d} is well-defined. We define \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+1+d}|_{K}=\mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}, if for all $(u_{h},\widehat{u}_{h})\in W_{h}\times M_{h}$, we have \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) is unique defined at every face of $\mathcal{T}_{h}$ due to (111a), (111b), and (111c). Then \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h})\in H^{1}(\Omega).

Lemma 7.2**.**

For all $(u_{h},\widehat{u}_{h})\in L^{2}(\mathcal{T}_{h})\times L^{2}(\partial\mathcal{T}_{h})$, we have the stability:

[TABLE]

where $S(K)$ is the set of all the simplex $K^{\star}\in\mathcal{T}_{h}$ such that $K^{\star}$ and $K$ has at least one common node, and $\partial S(K)$ is the set of all the faces of those simplex.

Proof.

To simplify the proof, we only give a proof for $d=3$, the proof of $d=2$ is similar. According to (111), we divide the Lagrange points on $K$ of degree $k+1+d$ into 4 parts, and the corresponding Lagrange basis denoted as $\{\phi_{1,j}\}_{j=1}^{N_{1}}$, $\{\phi_{2,j}\}_{j=1}^{N_{2}}$, $\{\phi_{3,j}\}_{j=1}^{N_{3}}$, and $\{\phi_{4,j}\}_{j=1}^{N_{4}}$, which are determined by (1)(2)(3)(4) in (111), respectively. It is known that $\phi_{4,j}|_{\partial K}=0$, since the corresponding Lagrange points are inside $K$. We also denote the dual basis of $\{\phi_{3,j}\}_{j=1}^{N_{3}}$ and $\{\phi_{4,j}\}_{j=1}^{N_{4}}$ as $\{\psi_{3,j}\}_{j=1}^{N_{3}}$ and $\{\psi_{4,j}\}_{j=1}^{N_{4}}$, respectively, such that

[TABLE]

where $\delta_{j,\ell}$ is the Kronecker delta.

A result in show that

[TABLE]

We can write \mathcal{I}_{K}^{\color[rgb]{0,0,0}k+1+d}(u_{h},\widehat{u}_{h}) as

[TABLE]

where

[TABLE]

By a scaling argument, one can get

[TABLE]

Again, by a scaling argument, for the Lagrange point $A_{i}$ on a face $F\subset\partial K$, and $A_{i}$ is also the vertex of $\mathcal{T}_{h}$, one can get

[TABLE]

and the similar for $a_{2,i}$, for the Lagrange point $B_{i}$ on an edge $E\subset\partial F$, $F\subset\partial K$:

[TABLE]

We use (114c), (116), (117), (112), (115a), (116), (117) and a scaling argument, to get

[TABLE]

We use (114d), (116), (117), (7), (112), (115a), (116), (117) and a scaling argument, to get

[TABLE]

Then desired result is followed by (113), (116), (117), (7), (7) and (115a) with $p=2$. ∎

Finally, we have the following estimation

Lemma 7.3**.**

For all $(u_{h},\widehat{u}_{h})\in W_{h}\times M_{h}$, we have

[TABLE]

Proof of Lemma 7.3.

Since \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+2+d}(u_{h},u_{h})=u_{h} for every $u_{h}\in V_{h}$, then (120a) follows from Lemma 7.2 and the fact that \mathcal{I}_{h}^{\color[rgb]{0,0,0}k+2+d} is linear. Then (120b) follows by an application of the inverse inequality. This completes the proof. ∎