A projection-based error analysis of HDG methods

Abstract

We introduce a new technique for the error analysis of hybridizable discontinuous Galerkin (HDG) methods. The technique relies on the use of a new projection whose design is inspired by the form of the numerical traces of the methods. This renders the analysis of the projections of the discretization errors simple and concise. By showing that these projections of the errors are bounded in terms of the distance between the solution and its projection, our studies of influence of the stabilization parameter are reduced to local analyses of approximation by the projection. We illustrate the technique on a specific HDG method applied to a model second-order elliptic problem.

1. Introduction

This paper is dedicated to presenting a new technique for the error analysis of an emerging class of numerical methods called hybridizable discontinuous Galerkin (HDG) methods 11. The main idea is to devise a projection that matches the form of the numerical traces of the HDG method. The analysis of the projection of the error then becomes simple, concise and independent of the particular choice of the stabilization parameters.

We do not aim for maximal generality, but rather to convey the advantages of our technique compared to other DG analyses. Hence we only consider one specific HDG method, namely, the LDG-hybridizable (LDG-H) method 11 for the model problem

$$\begin{alignat}{2} c\;{\boldsymbol{q}}+ \mathop{\boldsymbol{\nabla }}u & =0 &&\quad \text{ in }\Omega ,\tag{1.1a}\\ \mathop{\nabla }\cdot \,\, {\boldsymbol{q}}& = f&&\quad \text{ in }\Omega ,\tag{1.1b}\\ u & = g &&\quad \text{ on }\partial \Omega , \tag{1.1c} \end{alignat}$$

where $\Omega$ is a Lipschitz polyhedral domain in $\mathbb{R}^n$ ($n \ge 2$). Here $c:\Omega \mapsto \mathbb{R}^{n\times n}$ is a variable matrix-valued coefficient, which we assume to be symmetric and uniformly positive definite, $f$ is in $L^2(\Omega )$ and $g$ in $H^{1/2}(\partial \Omega ).$

Let us put our result in historical perspective. To do that, we begin by describing the HDG method. Consider a partitioning of the domain $\Omega$ into elements $K$ forming a mesh ${\mathcal{T}_h}$ satisfying the standard finite element conditions 7; the faces of $K$ are going to be denoted by $F$. This method yields a scalar approximation $u_h$ to $u$, a vector approximation ${\boldsymbol{q}}_h$ to ${\boldsymbol{q}}$, and a scalar approximation $\widehat{u}_h$ to the trace of $u$ on element boundaries, in spaces of the form

$$\begin{align} \boldsymbol{V}_h &= \{ {\boldsymbol{v}}: \text{ for all mesh elements }K, {\boldsymbol{v}}|_K \in \boldsymbol{V}(K)\}, {\tag{1.2a}}\\ W_h &= \{ w: \text{ for all mesh elements }K, w|_K \in W(K)\}, {\tag{1.2b}}\\ M_h &= \{ \mu : \text{ for all mesh faces } F, \;\mu |_F \in M(F)\}, {\tag{1.2c}} \end{align}$$

respectively, where $\boldsymbol{V}(K)$, $W(K)$, and $M(F)$ are finite dimensional spaces. The HDG approximations $u_h$ in $W_h$, ${\boldsymbol{q}}_h$ in $\boldsymbol{V}_h$, and the numerical trace $\widehat{u}_h$ in $M_h$, are determined by requiring that

$$\begin{alignat}{2} ( c\, {\boldsymbol{q}}_h,{\boldsymbol{r}})_{{\mathcal{T}_h}} - ( u_h,\mathop{\nabla }\cdot \,{\boldsymbol{r}})_{{\mathcal{T}_h}} + \langle {\widehat{u}_h,{\boldsymbol{r}}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} & = 0, {\tag{1.3a}}\\ -( {\boldsymbol{q}}_h, \mathop{\boldsymbol{\nabla }}w)_{{\mathcal{T}_h}} + \langle { \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}, w} \rangle _{{\partial \mathcal{T}_h}} & = (f, w)_{{\mathcal{T}_h}}, {\tag{1.3b}}\\ \langle {\widehat{u}_h,\mu } \rangle _{\partial \Omega } & = \langle {g, \mu } \rangle _{\partial \Omega }, {\tag{1.3c}}\\ \langle {\widehat{\boldsymbol{q}}_h \cdot {\boldsymbol{n}},\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } & = 0, {\tag{1.3d}} \end{alignat}$$

hold for all ${\boldsymbol{r}}$ in $\boldsymbol{V}_h,$ $w \in W_h,$ and $\mu \in M_h,$ with a specific $\widehat{\boldsymbol{q}}_h$ defined on ${\partial \mathcal{T}_h}=\{\partial K: K\in {\mathcal{T}_h}\}$; see 11. Above and throughout, we use the notation

$$\begin{equation*} (v,w)_{{\mathcal{T}_h}}=\sum _{K \in {\mathcal{T}_h}} (v,w)_K \qquad \text{and}\qquad \langle {v,w} \rangle _{{\partial \mathcal{T}_h}} = \sum _{K \in {\mathcal{T}_h}} \langle { v,w} \rangle _{\partial K}, \end{equation*}$$

where we write $(u,v)_D =\int _D u v \; dx$ whenever $D$ is a domain of $\mathbb{R}^n$, and $\langle {u,v} \rangle _D = \int _D u v \; dx$ whenever $D$ is a domain of $\mathbb{R}^{n-1}$. For vector functions ${\boldsymbol{v}}$ and ${\boldsymbol{w}}$, the notations are similarly defined with the integrand being the dot product ${\boldsymbol{v}}\cdot {\boldsymbol{w}}$. For HDG methods $\widehat{u}_h$ is taken to be an unknown in $M_h$, while the numerical trace $\widehat{\boldsymbol{q}}_h$ is prescribed on ${\partial \mathcal{T}_h}$ in such a way that both ${\boldsymbol{q}}_h$ and $u_h$ can be eliminated from the above equations to give rise to a single equation for $\widehat{u}_h$; see 11. Thus the often made criticism that DG methods have too many unknowns does not apply to HDG methods. Additional advantages of HDG methods include the ability to postprocess to get higher order solutions as we shall see in Section 5.

In contrast, the DG methods of the last century for second-order elliptic problems use only the first two equations 1.3a–1.3b, together with specific prescriptions of both $\widehat{u}_h$ and $\widehat{\boldsymbol{q}}_h$ to define their approximations. For most such methods, the above mentioned elimination is not feasible. Moreover, their analysis is also very different from ours, as it does not require the use of any special projection. This can be seen in 2, where most of the then known DG methods were analyzed in a single unifying framework. All of these methods were shown to converge with order $k+1$ in the scalar variable $u$ and with order $k$ in the flux ${\boldsymbol{q}}$. The order of convergence in $u$ is optimal and that in ${\boldsymbol{q}}$ sub-optimal since the methods use as local spaces $\boldsymbol{V}(K)$ and $W(K)$ the sets $\boldsymbol{\mathcal{P}}_k(K)$ and $\mathcal{P}_k(K)$, respectively, where $\boldsymbol{\mathcal{P}}_k(K):=[\mathcal{P}_k(K)]^n$ and $\mathcal{P}_k(K)$ is the space of polynomials of total degree at most $k$.

Two particular DG methods of the local discontinuous Galerkin (LDG) type, fitting in the above-mentioned unifying framework, deserve special mention as they provide approximations to the flux which converge with better orders of convergence. Their analyses do use special projections. The first is defined in the one-dimensional case and provides an order of convergence for the approximate flux of $k+1$; see 6. The second is an extension of that method to multiple dimensions. It uses Cartesian grids and takes as local spaces $\boldsymbol{V}(K)$ and $W(K)$ the sets $\boldsymbol{\mathcal{Q}}_k(K)$ and $\mathcal{Q}_k(K)$, respectively, where $\boldsymbol{\mathcal{Q}}_k(K):=[\mathcal{Q}_k(K)]^n$ and $\mathcal{Q}_k(K)$ is the space of polynomials of degree at most $k$ in each variable. The order of convergence for the approximate flux can be proven to be $k+1/2$; see 13. In both cases, the special projections used in the analysis were carefully devised to capture the structure of the numerical fluxes $\widehat{u}_h$ and $\widehat{\boldsymbol{q}}_h$.

Another example of a method requiring a special projection to carry out its analysis is the so-called minimal dissipation LDG method. This is a special LDG method whose penalization parameter is identically zero in all the interior faces; as a consequence, it cannot be analyzed as in 2. However, by using a suitably defined projection, the lack of stabilization can be overcome and sharp error estimates can be obtained. The flux can be shown to converge with the sub-optimal but sharp order of $k$; see 8.

Now, let us consider the HDG methods for which

$$\begin{equation} \widehat{\boldsymbol{q}}_h = {\boldsymbol{q}}_h + \tau \, \big ( u_h - \widehat{u}_h \big ) {\boldsymbol{n}}, \qquad \text{ on } {\partial \mathcal{T}_h}, {\tag{1.4}} \end{equation}$$

for some nonnegative penalty function $\tau$ defined on ${\partial \mathcal{T}_h}$ which we assume to be constant on each face of the triangulation. As explained in 11, these methods are called the LDG-hybridizable (LDG-H) methods because the above numerical trace is that of the LDG method applied separately on each mesh element $K$. The well-known hybridized versions of the Raviart-Thomas (RT) method 18 and the Brezzi-Douglas-Marini (BDM) method 5 can be considered to be HDG methods for which $\tau =0$; their analyses are carried out by using the celebrated RT and BDM projections. The analysis of HDG methods for which $\tau$ is not identically equal to zero, has been carried out in 12 also by using special projections. All of the above methods were proven to provide approximations to the flux converging with the optimal order of $k+1$; the first HDG method with this property was introduced in 9.

In this paper, we provide a new approach to the error analysis of HDG methods. It is an alternative to the techniques in 12 for general HDG methods, and those in 9 for the particular method treated therein. We recover all the results of 912 and obtain new superconvergence results for the projection of the scalar variable with our new approach. The novelty of our analysis is the use of new projections $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}, \varPi _{\!\scriptscriptstyle {W}}$ which are fitted to the structure of the numerical traces of the scheme in the sense that

$$\begin{equation*} P_M({\boldsymbol{q}}\cdot {\boldsymbol{n}})=\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}+\tau (\varPi _{\!\scriptscriptstyle {W}}u-P_M u) \quad \quad \text{ on } \partial \mathcal{T}_h, \end{equation*}$$

where $P_M$ is the $L^2$-projection into $M_h$; cf. the expression for $\widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}$ from 1.4.

Projections $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}, \varPi _{\!\scriptscriptstyle {W}}$ tailored to the numerical traces have actually been widely used. Indeed, beginning with the simplest example, notice that the projections used to analyze both the RT and BDM methods capture the structure of its numerical fluxes by satisfying

$$\begin{align*} P_M({\boldsymbol{q}}\cdot {\boldsymbol{n}}) & = \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}\quad \text{ on }{\partial \mathcal{T}_h}, \end{align*}$$

This should be compared with the definition of their numerical trace 11

$$\begin{align*} \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}& = {\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\qquad \text{ on } {\partial \mathcal{T}_h}. \end{align*}$$

Besides the projections used to carry out the analysis of the DG methods in 8, 6 and 13, we have the projections used in 9 to analyze an HDG type method called single-face hybridizable (SFH) method. The name arises due to the fact that for every simplex $K$, the penalty function $\tau$ is nonzero just on a single face, say $F_K$. The projections $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}$ and $\varPi _{\!\scriptscriptstyle {W}}$ used there are such that for each simplex $K\in {\mathcal{T}_h}$,

$$\begin{alignat*}{2} P_M({\boldsymbol{q}}\cdot {\boldsymbol{n}})=&\; \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}&&\quad \text{ on } \partial K\setminus F_K,\\ P_M u=&\; \varPi _{\!\scriptscriptstyle {W}}u &&\quad \text{ on } F_K. \end{alignat*}$$

This is in accordance with the numerical traces of the SFH method given by

$$\begin{alignat*}{2} \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}=&\; {\boldsymbol{q}}_h\cdot {\boldsymbol{n}}&&\quad \text{ on } \partial K\setminus F_K,\\ \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}=&\; {\boldsymbol{q}}_h\cdot {\boldsymbol{n}}+\tau (u_h-\widehat{u}_h) &&\quad \text{ on } F_K, \end{alignat*}$$

The SFH method is also an HDG method of the LDG-H type, so our projection can also be used to analyze the SFH method. It is interesting to point out that the component $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ of the SFH projection introduced in 9 and ours do coincide when $\tau$ is as above, whereas the component $\varPi _{\!\scriptscriptstyle {W}}u$ does not, except when $\nabla \cdot {\boldsymbol{q}}$ is a polynomial of degree $k-1$ on each simplex $K\in \mathcal{T}_h$. The final estimates for the SFH method are the same with both approaches.

The remainder of this paper is organized as follows. In Section 2, we define the new projection $\varPi _h$ and state its key properties. In Section 3, we use an energy argument to obtain an optimal error estimate for the approximate flux ${\boldsymbol{q}}_h$ and its numerical trace $\widehat{\boldsymbol{q}}_h$. In Section 4, we use a duality argument to obtain a superconvergence estimate of the projection of the error in $u_h$ and its numerical trace $\widehat{u}_h$. In Section 5 we discuss previously known element-by-element postprocessing of the flux and the scalar variable. They will result in further approximations ${\boldsymbol{q}}^\star _h$ and $u^\star _h$ with interesting properties. Finally, we conclude in Section 6 by summarizing our results and relating them to the previous work 12. Proof of the properties of the projection have been gathered in the Appendix.

2. The projection

The main ingredient of our error analysis is a new projection $\varPi _h$ into the product space $\boldsymbol{V}_h\times W_h$. In this section, we introduce it and establish its properties. From now on, we are going to use the following local spaces:

$$\begin{equation*} \boldsymbol{V}(K)=\boldsymbol{\mathcal{P}}_k(K), \quad W(K)=\mathcal{P}_k(K) \quad \text{ and } \quad M(F)=\mathcal{P}_k(F). \end{equation*}$$

Also, we consider the stabilization function $\tau :\partial \mathcal{T}_h\to \mathbb{R}$ to be constant on each face.

The projected function is denoted by $\varPi _h ({\boldsymbol{q}}, u)$ or by $(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}, \varPi _{\!\scriptscriptstyle {W}}u)$ where $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ and $\varPi _{\!\scriptscriptstyle {W}}u$ are the components of the projection in $\boldsymbol{V}_h$ and $W_h$, respectively. The values of the projection on any simplex $K$ are fixed by requiring that the components satisfy the equations

$$\begin{align} (\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}, {\boldsymbol{v}})_K &= ( {\boldsymbol{q}}, {\boldsymbol{v}})_K && \quad \text{for all }{\boldsymbol{v}}\in \boldsymbol{\mathcal{P}}_{k-1}(K), {\tag{2.1a}}\\ (\varPi _{\!\scriptscriptstyle {W}}u, w )_K &= ( u, w )_K && \quad \text{for all }w \in \mathcal{P}_{k-1}(K), {\tag{2.1b}}\\ \langle {\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}+ \tau \varPi _{\!\scriptscriptstyle {W}}u,\mu } \rangle _F &= \langle {{\boldsymbol{q}}\cdot {\boldsymbol{n}}+ \tau u,\mu } \rangle _F && \quad \text{for all }\mu \in \mathcal{P}_k(F), {\tag{2.1c}} \end{align}$$

for all faces $F$ of the simplex $K$. If $k=0$, then 2.1a and 2.1b are vacuous and $\varPi _h$ is defined solely by 2.1c. Note that although we denoted the first component of the projection by $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$, it depends not just on ${\boldsymbol{q}}$, but rather on both ${\boldsymbol{q}}$ and $u$, as we see from 2.1. The same is true for $\varPi _{\!\scriptscriptstyle {W}}u$. Hence the notation $(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}, \varPi _{\!\scriptscriptstyle {W}}u)$ for $\varPi _h ({\boldsymbol{q}}, u)$ is somewhat misleading, but its convenience outweighs this disadvantage.

The domain of $\varPi _h$ is a subspace of $L^2(\Omega )^n \times L^2(\Omega )$ on which the right-hand sides of 2.1 are well defined. Indeed, all functions ${\boldsymbol{q}}$ and $u$ that are regular enough for their traces ${\boldsymbol{q}}\cdot {\boldsymbol{n}}$ and $u$ to be in $L^2(\partial K)$ are in the domain of $\varPi _h$, for example, $({\boldsymbol{q}},u)\in \boldsymbol{H}^1({\mathcal{T}_h})\times H^1({\mathcal{T}_h})$, where

$$\begin{equation*} H^1({\mathcal{T}_h}) = \prod _{K\in {\mathcal{T}_h}} H^1(K), \qquad \boldsymbol{H}^1({\mathcal{T}_h})= H^1({\mathcal{T}_h})^n. \end{equation*}$$

That the left-hand sides of 2.1 uniquely determine $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ and $\varPi _{\!\scriptscriptstyle {W}}u$ is proved as part of the next theorem.

To state the theorem, we need to introduce additional notation and conventions. We use $\| \cdot \|_D$ to denote the $L^2(D)$-norm for any $D$. More generally, we denote the norm and seminorm on any Sobolev space $X$ by $\| \cdot \|_X$ and $|\cdot |_X$, respectively. Also, as is usual in finite element analysis, we restrict elements to have shape regularity, i.e., all mesh elements $K$ under consideration satisfy $h_K/\rho _K\le \gamma$, where $h_K=\operatorname {diam}(K)$, $\rho _K$ is the diameter of the largest ball contained in $K$, and $\gamma$ is a fixed constant. We use $C$, with or without subscripts, to denote a generic constant, independent of the elements and functions involved in our inequalities (but it may depend on $\gamma$). The value of $C$ at different occurrences may differ. While we absorb dependencies on $c$ into $C$, the dependencies on $\tau$ will always be explicitly mentioned.

Our first result states that the projection $\varPi _h$ is well defined and has reasonable approximation properties. See Appendix A for a detailed proof.

Theorem 2.1.

Suppose $k\ge 0$, $\tau |_{\partial K}$ is nonnegative and $\tau ^{\max }_K:=\max \tau |_{\partial K} >0$. Then the system 2.1 is uniquely solvable for $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}_h$ and $\varPi _{\!\scriptscriptstyle {W}}u$. Furthermore, there is a constant $C_{\!\scriptscriptstyle {\varPi }}$ independent of $K$ and $\tau$ such that

$$\begin{alignat*}{1} \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\, \|_K\le &\;C\,h^{\ell _{\boldsymbol{q}}+1}_K\,|{\boldsymbol{q}}|_{\boldsymbol{H}^{\ell _{\boldsymbol{q}}+1}(K)} +C\,{h^{\ell _u+1}_K}\,{\tau _K^{*}}\,|u|_{H^{\ell _u+1}(K)}, \\ \| \varPi _{\!\scriptscriptstyle {W}}u - u \|_K \le &\;C\,h^{\ell _u+1}_K\,|u|_{H^{\ell _u+1}(K)} +C\,\frac{h^{\ell _{\boldsymbol{q}}+1}_K}{\tau _K^{\max }}\,|\nabla \cdot {\boldsymbol{q}}|_{H^{\ell _{\boldsymbol{q}}}(K)}, \end{alignat*}$$

for $\ell _u, \ell _{\boldsymbol{q}}$ in $[0,k]$. Here $\tau _K^{*}:=\max \tau |_{\partial K\setminus F^*}$, where $F^*$ is a face of $K$ at which $\tau |_{\partial K}$ is maximum.

Note that if $\tau$ is of unit order, both approximation errors converge with the optimal order of $k+1$, when the functions ${\boldsymbol{q}}$ and $u$ are smooth enough. If $\tau \equiv 1$, or if $\tau$ is bounded above and below uniformly by fixed constants, then a standard Bramble-Hilbert argument proves the order $k+1$ approximation estimates of the theorem. However, to track the dependence of the constants on $\tau$, and to compare with previous works, we need to perform a more careful analysis. This is done in Appendix A in elaborate detail. To make an immediate comparison, note that for the $\tau$ used in the SFH method 9, $\tau _K^{*}=0$ and the approximation properties of $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ are independent of the stabilization function $\tau$.

We conclude this section with a property of the projection that we use critically in the error analysis of the method.

Proposition 2.1 (A weak commutativity property).

For any $w$ in $W_h$ and any $(\boldsymbol{\varPhi },\varPsi )$ in the domain of $\varPi _h$, we have

$$\begin{alignat*}{2} (w ,\mathop{\nabla }\cdot \,\boldsymbol{\varPhi })_K =&\; (w, \mathop{\nabla }\cdot \,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_K + \langle { w, \tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi ) } \rangle _{\partial K}. \end{alignat*}$$

Proof.

For any $w\in W_h$,

$$\begin{equation*} \begin{aligned} (w, \mathop{\nabla }\cdot \,\boldsymbol{\varPhi })_K & = -(\mathop{\boldsymbol{\nabla }}w, \boldsymbol{\varPhi })_K + \langle {w, \boldsymbol{\varPhi }\cdot {\boldsymbol{n}}} \rangle _{\partial K} \\ & = -(\mathop{\boldsymbol{\nabla }}w, \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_K + \langle {w, \boldsymbol{\varPhi }\cdot {\boldsymbol{n}}} \rangle _{\partial K} &&\quad \text{by \xhref[disp-formula]{#texmlid4}{2.1a},} \\ & =-(\mathop{\boldsymbol{\nabla }}w, \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_K + \langle {w,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }\cdot {\boldsymbol{n}}+ \tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi ) } \rangle _{\partial K} &&\quad \text{by \xhref[disp-formula]{#texmlid6}{2.1c},} \\ & = ( w, \mathop{\nabla }\cdot \,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_K + \langle {w,\tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi ) } \rangle _{\partial K}. && & \end{aligned} \end{equation*}$$ ■

3. Flux error estimates by an energy argument

The purpose of this section is to give error estimates under minimal regularity assumptions on the solution. Since the projection $\varPi _h$ is designed to fit the structure of the numerical trace, the equations satisfied by the projection of the errors have a form amenable to simple analysis, as we see next. We will need the $L^2$-orthogonal projection onto $M_h$, which we denote by $P_M$. Note that because $\tau$ is piecewise constant,

$$\begin{equation*} \langle {\tau (P_Mu-u),\mu } \rangle _{{\partial \mathcal{T}_h}}=0 \quad \text{for all }\mu \in M_h, \end{equation*}$$

a fact that we will use repeatedly without explicit mention. The projection of the errors satisfy the following.

Lemma 3.1.

Let $\boldsymbol{\varepsilon }_h^{\;q}= \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}_h$, $\varepsilon _h^u= \varPi _{\!\scriptscriptstyle {W}}u - u_h$, and $\varepsilon _h^{\widehat{u}}= P_{\scriptscriptstyle {M}}u - \widehat{u}_h.$ Then

$$\begin{alignat}{2} ( c\, \boldsymbol{\varepsilon }_h^{\;q},{\boldsymbol{r}})_{{\mathcal{T}_h}} - ( \varepsilon _h^u,\mathop{\nabla }\cdot \,{\boldsymbol{r}})_{{\mathcal{T}_h}} + \langle {\varepsilon _h^{\widehat{u}},{\boldsymbol{r}}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} & = (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}),{\boldsymbol{r}}\,)_{\mathcal{T}_h}, {\tag{3.1a}}\\ -( \boldsymbol{\varepsilon }_h^{\;q}, \mathop{\boldsymbol{\nabla }}w)_{{\mathcal{T}_h}} + \langle { \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}, w} \rangle _{{\partial \mathcal{T}_h}} & = 0, {\tag{3.1b}}\\ \langle {\varepsilon _h^{\widehat{u}},\mu } \rangle _{\partial \Omega } & = 0, {\tag{3.1c}}\\ \langle {\widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}},\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } & = 0, {\tag{3.1d}} \end{alignat}$$for all ${\boldsymbol{r}}\in \boldsymbol{V}_h, \, w\in W_h,$ and $\mu \in M_h$, where$$\begin{alignat}{2} \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}:= \boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}+ \tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}) = P_{\scriptscriptstyle {M}}({\boldsymbol{q}}\cdot {\boldsymbol{n}}) & -\widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\qquad \text{ on }\partial {\mathcal{T}_h}\setminus \partial \Omega . {\tag{3.1e}} \end{alignat}$$

Proof.

Let us begin by noting that the exact solution ${\boldsymbol{q}}$ and $u$ satisfies

$$\begin{alignat*}{2} ( c\, {\boldsymbol{q}},{\boldsymbol{r}})_{{\mathcal{T}_h}} - ( u,\mathop{\nabla }\cdot \,{\boldsymbol{r}})_{{\mathcal{T}_h}} + \langle {u,{\boldsymbol{r}}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} & = 0, \\ -( {\boldsymbol{q}}, \mathop{\boldsymbol{\nabla }}w)_{{\mathcal{T}_h}} + \langle { {\boldsymbol{q}}\cdot {\boldsymbol{n}}, w} \rangle _{{\partial \mathcal{T}_h}} & = (f, w)_{{\mathcal{T}_h}}, \end{alignat*}$$

for all ${\boldsymbol{r}}\in \boldsymbol{V}_h$ and $w \in W_h$. By the definition of $\varPi _h$ and $P_{\scriptscriptstyle {M}}$, the above implies

$$\begin{alignat*}{2} ( c\, \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}},{\boldsymbol{r}})_{{\mathcal{T}_h}} - ( \varPi _{\!\scriptscriptstyle {W}}u,\mathop{\nabla }\cdot \,{\boldsymbol{r}})_{{\mathcal{T}_h}} + \langle {P_{\scriptscriptstyle {M}}u,{\boldsymbol{r}}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} & = (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}),{\boldsymbol{r}})_{\mathcal{T}_h}, \\ -( \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}, \mathop{\boldsymbol{\nabla }}w)_{{\mathcal{T}_h}} + \langle { \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}- \tau (u-\varPi _{\!\scriptscriptstyle {W}}u), w} \rangle _{{\partial \mathcal{T}_h}} & = (f, w)_{{\mathcal{T}_h}}, \end{alignat*}$$

for all ${\boldsymbol{r}}\in \boldsymbol{V}_h$ and $w \in W_h$. Subtracting 1.3a and 1.3b from the above two equations, respectively, we obtain 3.1a and 3.1b. The equation 3.1c follows directly from the boundary condition 1.3c. To prove 3.1d we proceed as follows:

$$\begin{align*} \langle {\mu , \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } & = \langle { (\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}_h )\cdot {\boldsymbol{n}}+ \tau (\varPi _{\!\scriptscriptstyle {W}}u-u_h - P_{\scriptscriptstyle {M}}u + \widehat{u}_h),\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } \\ &= \langle { ({\boldsymbol{q}}- {\boldsymbol{q}}_h )\cdot {\boldsymbol{n}}+ \tau (u-u_h - u + \widehat{u}_h),\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } \\ &= \langle { {\boldsymbol{q}}\cdot {\boldsymbol{n}},\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega } - \langle { \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}},\mu } \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega }, \end{align*}$$

where we have used the definition of $\varPi _h$. Note that this proves the identity 3.1e because $\widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\in M_h$. Since ${\boldsymbol{q}}$ is in $H(\mathrm{div},{\Omega })$ and since $\widehat{\boldsymbol{q}}_h$ satisfies 1.3d, both terms above are zero. This completes the proof.

■

Lemma 3.2.

We have

$$\begin{equation*} (c\,\boldsymbol{\varepsilon }_h^{\;q},\boldsymbol{\varepsilon }_h^{\;q})_{\mathcal{T}_h}+ \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),(\varepsilon _h^u-\varepsilon _h^{\widehat{u}})} \rangle _{\partial \mathcal{T}_h}= (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}), \boldsymbol{\varepsilon }_h^{\;q})_{\mathcal{T}_h}. \end{equation*}$$

Proof.

Taking ${\boldsymbol{r}}:=\boldsymbol{\varepsilon }_h^{\;q}$ in 3.1a, $w=\varepsilon _h^u$ in 3.1b, $\mu =-\widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}$ in 3.1c and $\mu =-\varepsilon _h^{\widehat{u}}$ in 3.1d, and adding the resulting four equations, we get

$$\begin{equation*} (c\,\boldsymbol{\varepsilon }_h^{\;q},\boldsymbol{\varepsilon }_h^{\;q})_{\mathcal{T}_h}+ \Theta _h = \; (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}), \boldsymbol{\varepsilon }_h^{\;q})_{\mathcal{T}_h}, \end{equation*}$$

where

$$\begin{align*} \Theta _h=&-( \varepsilon _h^u, \mathop{\nabla }\cdot \,\boldsymbol{\varepsilon }_h^{\;q})_{{\mathcal{T}_h}} + \langle {\varepsilon _h^{\widehat{u}},\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} -( \boldsymbol{\varepsilon }_h^{\;q}, \mathop{\boldsymbol{\nabla }}\varepsilon _h^u)_{{\mathcal{T}_h}} \\ &+ \langle { \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}, \varepsilon _h^u} \rangle _{{\partial \mathcal{T}_h}} - \langle {\varepsilon _h^{\widehat{u}}, \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}}. \end{align*}$$

But, by integration by parts, we get

$$\begin{alignat*}{1} \Theta _h =&-\langle { \varepsilon _h^u, \boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}} \rangle _{{\mathcal{T}_h}} + \langle {\varepsilon _h^{\widehat{u}},\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} + \langle { \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}, \varepsilon _h^u} \rangle _{{\partial \mathcal{T}_h}} - \langle {\varepsilon _h^{\widehat{u}}, \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} \\ =&-\langle {\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}, \varepsilon _h^u-\varepsilon _h^{\widehat{u}}} \rangle _{{\mathcal{T}_h}} + \langle { \widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}, \varepsilon _h^u-\varepsilon _h^{\widehat{u}}} \rangle _{{\partial \mathcal{T}_h}} \\ =&\; \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),(\varepsilon _h^u-\varepsilon _h^{\widehat{u}})} \rangle _{\partial \mathcal{T}_h}, \end{alignat*}$$

by the definition of $\widehat{\boldsymbol{\varepsilon }}_h$, namely 3.1e. This completes the proof.

■

While Lemma 3.2 identifies an “energy” norm, Lemma 3.1 gives “consistency” relations. These are enough to immediately prove an error estimate for the flux. To state it, we need the norm $\| \cdot \|_h$ defined by $\|\mu \|_h^2 = \sum _{K\in {\mathcal{T}_h}} h_K\| \mu \|_{\partial K}^2$ for any function $\mu \in L^2({\partial \mathcal{T}_h}):=\prod _{K\in {\mathcal{T}_h}} L^2(\partial K)$, and the $c$-weighted $L^2(\Omega )$-norm $\| {\boldsymbol{q}}\|_c:= (c\,{\boldsymbol{q}},{\boldsymbol{q}})_{{\mathcal{T}_h}}^{1/2}.$

Theorem 3.1 (Flux error estimates).

Let ${\boldsymbol{q}}_h,\; u_h,$ and $\widehat{u}_h$ solve the HDG equations 1.3 and let the exact solution ${\boldsymbol{q}}, u$ be in the domain of $\varPi _h$. Then, for $k\ge 0$,

$$\begin{align} \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}_h \|_c & \,\le \, \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\|_c, {\tag{3.2}}\\ \|P_{\scriptscriptstyle {M}}({\boldsymbol{q}}\cdot {\boldsymbol{n}}) - \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\|_h & \,\le \, C_{1,\tau }\,\| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\|_c, {\tag{3.3}} \end{align}$$

where $C_{1,\tau }=C\,\max \{ 1,\left(h_K\,\tau _K^{\max }\right)^{1/2}: K\in {\mathcal{T}_h}\}$.

Proof.

By the Cauchy-Schwarz inequality applied to the identity of Lemma 3.2, we get

$$\begin{equation} (c\,\boldsymbol{\varepsilon }_h^{\;q},\boldsymbol{\varepsilon }_h^{\;q})_{\mathcal{T}_h}+ \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),(\varepsilon _h^u-\varepsilon _h^{\widehat{u}})} \rangle _{\partial \mathcal{T}_h}\le \|\,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}\,\|_c\,\|\boldsymbol{\varepsilon }_h^{\;q}\|_c, {\tag{3.4}} \end{equation}$$

and 3.2 follows since $\tau \ge 0$. For 3.3, we start by using the identity 3.1e to get

$$\begin{alignat*}{1} \| P_{\scriptscriptstyle {M}}({\boldsymbol{q}}\cdot {\boldsymbol{n}}) - \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\|_h & \le \; \| \boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}}\|_h +\| \tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}) \|_h \\ & \le \; C\, \|\boldsymbol{\varepsilon }_h^{\;q}\|_{c} +\max _{K\in {\mathcal{T}_h}} \left(h_K\,\tau ^{\max }_K\right)^{1/2}\,\langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),(\varepsilon _h^u-\varepsilon _h^{\widehat{u}})} \rangle _{\partial \mathcal{T}_h}^{1/2} \end{alignat*}$$

where we used an inverse inequality and the fact that $c^{-1}$ is uniformly bounded. Then 3.3 follows from 3.4.

■

4. Superconvergence of the scalar variableby a duality argument

In the previous section, we established an error estimate for ${\boldsymbol{q}}_h$ that only requires that the solution be as regular as required for the application of the projection. On domains permitting higher regularity estimates, we can perform an analogue of the Aubin-Nitsche duality argument 317 to get higher rates of convergence. In particular, such arguments will give us error estimates for $u_h$ and $\widehat{u}_h$.

We thus begin by introducing the dual problem for any given $\varTheta$ in $L^2(\Omega )$:

$$\begin{align} c\,\boldsymbol{\varPhi }- \mathop{\boldsymbol{\nabla }}\varPsi & = 0 &&\text{ on } \Omega , {\tag{4.1a}}\\ \mathop{\nabla }\cdot \,\boldsymbol{\varPhi }& = \varTheta &&\text{ on } \Omega , {\tag{4.1b}}\\ \varPsi &=0 &&\text{ on } \partial \Omega . {\tag{4.1c}} \end{align}$$

We assume that this boundary value problem admits the regularity estimate

$$\begin{equation} \| \boldsymbol{\varPhi }\,\|_{H^1(\Omega )} + \| \varPsi \|_{H^{2}(\Omega )} \le C_{\!\mathrm{reg}}\|\varTheta \|_\Omega {\tag{4.2}} \end{equation}$$

for all $\varTheta$ in $L^2(\Omega )$. This is well known to hold in several cases, e.g., if $c\equiv 1$ and $\Omega$ is a convex polygon 16. Recall that we have been tacitly assuming that $({\boldsymbol{q}},u)$ is in the domain of $\varPi _h$. By 4.2, $(\boldsymbol{\varPhi },\varPsi )$ is also regular enough to apply $\varPi _h$, so we have the following lemma.

Lemma 4.1 (The duality argument).

For any $\varPsi _h \in W_h$, we have

$$\begin{equation*} (\varepsilon _h^u,\varTheta )_{\mathcal{T}_h}= \; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi })_{\mathcal{T}_h}+ ({\boldsymbol{q}}-\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}},\mathop{\boldsymbol{\nabla }}\varPsi -\mathop{\boldsymbol{\nabla }}\varPsi _h)_{\mathcal{T}_h}. \end{equation*}$$

Consequently,

$$\begin{equation*} \|\,\varepsilon _h^u\,\|_{\mathcal{T}_h}\le C\,H(\varTheta )\,\|\,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}\,\|_{\mathcal{T}_h}, \end{equation*}$$

where

$$\begin{equation*} H(\varTheta ):=\sup _{\varTheta \in L^2(\Omega )\setminus {0}} \frac{\|\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi }\|_{\mathcal{T}_h}}{\|\varTheta \|_{\mathcal{T}_h}} + \sup _{\varTheta \in L^2(\Omega )\setminus {0}}\inf _{\varPsi _h\in W_h} \frac{\|\mathop{\boldsymbol{\nabla }}\varPsi -\mathop{\boldsymbol{\nabla }}\varPsi _h\,\|_{\mathcal{T}_h}}{\|\varTheta \|_{\mathcal{T}_h}}. \end{equation*}$$

Proof.

We have

$$\begin{alignat*}{2} (\varepsilon _h^u,\varTheta )_{\mathcal{T}_h}=&\; (\varepsilon _h^u,\mathop{\nabla }\cdot \,\boldsymbol{\varPhi })_{\mathcal{T}_h}&&\text{ by \xhref[disp-formula]{#texmlid20}{4.1b},} \\ =&\;(\varepsilon _h^u,\mathop{\nabla }\cdot \,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }) + \langle { \varepsilon _h^u, \tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi )} \rangle _{{\partial \mathcal{T}_h}} &&\text{ by Prop. \xhref[statement]{#ltxid10}{2.1},} \\ =&\; (c\,\boldsymbol{\varepsilon }_h^{\;q},\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\varepsilon _h^{\widehat{u}},\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}} \\ & - (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\varepsilon _h^u,\tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi )} \rangle _{\partial \mathcal{T}_h}\quad &&\text{ by \xhref[disp-formula]{#texmlid7}{3.1a},} \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\varepsilon _h^{\widehat{u}},(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi })\cdot {\boldsymbol{n}}} \rangle _{\partial \mathcal{T}_h}\\ & +\langle {\varepsilon _h^u,\tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi )} \rangle _{\partial \mathcal{T}_h}, \end{alignat*}$$

by the continuity of $\boldsymbol{\varPhi }\cdot {\boldsymbol{n}}$ and the fact that $\varepsilon _h^{\widehat{u}}=0$ on $\partial \Omega$ by 3.1c. Then

$$\begin{align*} (\varepsilon _h^u,\varTheta )_{\mathcal{T}_h}=&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\varepsilon _h^u-\varepsilon _h^{\widehat{u}},\tau (\varPi _{\!\scriptscriptstyle {W}}\varPsi - \varPsi )} \rangle _{\partial \mathcal{T}_h}&&\text{ by \xhref[disp-formula]{#texmlid6}{2.1c},} \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}), \varPi _{\!\scriptscriptstyle {W}}\varPsi } \rangle _{\partial \mathcal{T}_h}\\ &-\langle {\varepsilon _h^u-\varepsilon _h^{\widehat{u}},\tau \,P_{\scriptscriptstyle {M}}\varPsi } \rangle _{\partial \mathcal{T}_h}\\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),\varPi _{\!\scriptscriptstyle {W}}\varPsi } \rangle _{\partial \mathcal{T}_h}\\ & -\langle {\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}},P_{\scriptscriptstyle {M}}\varPsi } \rangle _{\partial \mathcal{T}_h}&&\text{ by \xhref[disp-formula]{#texmlid11}{3.1d},} \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ \langle {\tau (\varepsilon _h^u-\varepsilon _h^{\widehat{u}}),\varPi _{\!\scriptscriptstyle {W}}\varPsi } \rangle _{\partial \mathcal{T}_h}\\ & -\langle {\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}},\varPsi } \rangle _{\partial \mathcal{T}_h}. \end{align*}$$

Moreover,

$$\begin{align*} (\varepsilon _h^u,\varTheta )_{\mathcal{T}_h}=&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ (\mathop{\nabla }\cdot \,\boldsymbol{\varepsilon }_h^{\;q},\varPi _{\!\scriptscriptstyle {W}}\varPsi )_{\mathcal{T}_h}-\langle {\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}},\varPsi } \rangle _{\partial \mathcal{T}_h}&&\text{ by \xhref[disp-formula]{#texmlid8}{3.1b},} \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}+ (\mathop{\nabla }\cdot \,\boldsymbol{\varepsilon }_h^{\;q},\varPsi )_{\mathcal{T}_h}-\langle {\boldsymbol{\varepsilon }_h^{\;q}\cdot {\boldsymbol{n}},\varPsi } \rangle _{\partial \mathcal{T}_h}&&\text{ by \xhref[disp-formula]{#texmlid5}{2.1b}}, \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi })_{\mathcal{T}_h}-(\boldsymbol{\varepsilon }_h^{\;q}, \mathop{\boldsymbol{\nabla }}\varPsi )_{\mathcal{T}_h}. \end{align*}$$

To bring this into the needed form, we continue:

$$\begin{align*} (\varepsilon _h^u,\varTheta )_{\mathcal{T}_h}=&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi })_{\mathcal{T}_h}+ (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPhi })_{\mathcal{T}_h}\\ &- (\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}_h, \mathop{\boldsymbol{\nabla }}\varPsi )_{\mathcal{T}_h} \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi })_{\mathcal{T}_h}+ ({\boldsymbol{q}}-\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}},\mathop{\boldsymbol{\nabla }}\varPsi )_{\mathcal{T}_h} && \text{ by \xhref[disp-formula]{#texmlid21}{4.1a}}, \\ =&\; (c\,({\boldsymbol{q}}-{\boldsymbol{q}}_h),\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}\boldsymbol{\varPhi }-\boldsymbol{\varPhi })_{\mathcal{T}_h}+ ({\boldsymbol{q}}-\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}},\mathop{\boldsymbol{\nabla }}\varPsi -\mathop{\boldsymbol{\nabla }}\varPsi _h)_{\mathcal{T}_h} &&\text{ by \xhref[disp-formula]{#texmlid4}{2.1a}}. \end{align*}$$

Finally, the inequality of the lemma follows by applying the Cauchy-Schwarz inequality to the identity and using the first estimate of Theorem 3.1.

■

Theorem 4.1.

Suppose the regularity assumption 4.2 holds. Then

$$\begin{align*} \| \varPi _{\!\scriptscriptstyle {W}}u - u_h \|_{\mathcal{T}_h}& \le C_{2,\tau } \, h^{\min \{k,1\}}\, \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\,\|_{\mathcal{T}_h}&&\text{ for } k\ge 0,\\ \| P_{\scriptscriptstyle {M}}u -\widehat{u}_h\|_h\;\, & \le C_{2,\tau }\, h \, \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\,\|_{\mathcal{T}_h}&&\text{ for } k\ge 1, \end{align*}$$

where $h= \max \{ h_K: K\in {\mathcal{T}_h}\}$ and $C_{2,\tau }=C\max \{ 1, h_K\tau _K^*: K\in {\mathcal{T}_h}\}$.

Proof.

Let us prove the first inequality. By the first estimate of Theorem 2.1 with $\ell _{{\boldsymbol{q}}}$ set to 0 and $\ell _u=\min \{k,1\}$, we get that

$$\begin{align*} H(\varTheta ) &\le \,C\,h^{\min \{k,1\}}\, C_{2,\tau } \,\sup _{\varTheta \in L^2(\Omega )\setminus {0}} \frac{\|\boldsymbol{\varPhi }\|_{H^1(\Omega )} + \| \varPsi \|_{H^{\ell _u+1}(\Omega )}}{\|\varTheta \|_\Omega } \\ &\le \,C\,h^{\min \{k,1\}}\, C_{2,\tau } \,C_{\!\mathrm{reg}}. \end{align*}$$

by the regularity estimate 4.2, and the first estimate follows.

The second estimate follows from the first from the same local argument used in 5 to obtain a similar estimate for the BDM method. Indeed, when $k\ge 1$, we can select a function ${\boldsymbol{r}}\in \boldsymbol{\mathcal{P}}_k(K)$ such that ${\boldsymbol{r}}\cdot {\boldsymbol{n}}= \varepsilon _h^{\widehat{u}}$ on $\partial K$ and $\| {\boldsymbol{r}}\|_K \le Ch_K^{1/2}\| \varepsilon _h^{\widehat{u}}\|_{\partial K}.$ Using $h_K\,{\boldsymbol{r}}$ as the test function in 3.1a, and applying an inverse inequality, we find that

$$\begin{align*} h_K \| \varepsilon _h^{\widehat{u}}\|_{\partial K}^2 &= h_K (c\,(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}),{\boldsymbol{r}})_K + h_K (\varepsilon _h^u,\mathop{\nabla }\cdot \,{\boldsymbol{r}})_K - h_K(c\,\boldsymbol{\varepsilon }_h^{\;q},{\boldsymbol{r}}\,)_K \\ &\le C h_K \| {\boldsymbol{r}}\|_K \big (\| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}\|_K + \|\boldsymbol{\varepsilon }_h^{\;q}\|_K \big ) + C \| {\boldsymbol{r}}\|_K \| \varepsilon _h^u\|_K. \end{align*}$$

Applying the first estimate of the theorem for $\varepsilon _h^u$ and the estimate in Theorem 3.1 for $\boldsymbol{\varepsilon }_h^{\;q}$, we get the second inequality of the theorem. This completes the proof.

■

5. Enhanced accuracy by postprocessing

In this section we describe a few techniques to postprocess the approximate solution and flux.

5.1. Flux postprocessing

We can obtain a postprocessed flux ${\boldsymbol{q}}^\star _h$ with better conservation properties. Although ${\boldsymbol{q}}^\star _h$ converges at the same order as ${\boldsymbol{q}}_h$, it is in $H(\mathrm{div},\Omega )$ and its divergence converges at one higher order than ${\boldsymbol{q}}_h$.

To define ${\boldsymbol{q}}^\star _h$, we use a slight modification of the Raviart-Thomas projection 18, as used in the framework of Darcy flows 4, or for the Navier-Stokes equations 14, or for the postprocessing of superconvergent DG methods for second-order elliptic problems 12. We follow 12. Thus, on each simplex $K \in {\mathcal{T}_h}$, we take ${\boldsymbol{q}}_h^\star :={\boldsymbol{q}}_h+\boldsymbol{\eta }_h$ where $\boldsymbol{\eta }_h$ is the only element of $\boldsymbol{\mathcal{P}}_k(K)+{\boldsymbol{x}}\,\mathcal{P}_k(K)$ satisfying

$$\begin{alignat*}{2} (\boldsymbol{\eta }_h,{\boldsymbol{v}})_K&=0 &&\quad \text{ for all } {\boldsymbol{v}}\in \boldsymbol{\mathcal{P}}_{k-1}(K), \\ \langle \boldsymbol{\eta }_h\cdot {\boldsymbol{n}}, \mu \rangle _F&=\langle (\widehat{\boldsymbol{q}}_h-{\boldsymbol{q}}_h)\cdot {\boldsymbol{n}}, \mu \rangle _F && \quad \text{ for all }\mu \in \mathcal{P}_k(F)\text{ and all faces}\; F \;\text{of}\; K. \end{alignat*}$$

If $L^2(K)$-orthogonal hierarchical basis functions are used, computation of $\boldsymbol{\eta }_h$ reduces to solving, on each element, a linear system of size $n+1$ times the dimension of the space of homogeneous polynomials of degree $k$. It is instructive to compare the above definition to that of the Raviart-Thomas projection, namely on each simplex $K \in {\mathcal{T}_h}$, we set $\boldsymbol{\varPi }_k^{{\scriptscriptstyle {\mathrm{RT}}}}{\boldsymbol{q}}$ as the only element of $\boldsymbol{\mathcal{P}}_k(K)+{\boldsymbol{x}}\,\mathcal{P}_k(K)$ satisfying

$$\begin{alignat*}{2} (\boldsymbol{\varPi }_k^{{\scriptscriptstyle {\mathrm{RT}}}}{\boldsymbol{q}}-{\boldsymbol{q}},{\boldsymbol{v}})_K&=0 && \quad \text{ for all } {\boldsymbol{v}}\in \boldsymbol{\mathcal{P}}_{k-1}(K), \\ \langle (\boldsymbol{\varPi }_k^{{\scriptscriptstyle {\mathrm{RT}}}}{\boldsymbol{q}}-{\boldsymbol{q}})\cdot {\boldsymbol{n}}, \mu \rangle _F&=0, &&\quad \text{ for all }\mu \in \mathcal{P}_k(F)\text{ and all faces}\; F \;\text{of}\; K. \end{alignat*}$$

This comparison yields the following theorem. For a proof, see 12.

Theorem 5.1.

For any $k\ge 0$, we have that ${\boldsymbol{q}}_h^\star \in H(\text{div},\Omega )$. Moreover,

$$\begin{alignat*}{1} \|\,{\boldsymbol{q}}-{\boldsymbol{q}}_h^\star \|_{\Omega } &\;\le \, C \|\,{\boldsymbol{q}}-\boldsymbol{\varPi }_k^{{\scriptscriptstyle {\mathrm{RT}}}}{\boldsymbol{q}}\,\|_{\Omega }+C_{1,\tau }\;\|\,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}\,\|_{{\mathcal{T}_h}}, \\ \|\nabla \cdot ({\boldsymbol{q}}- {\boldsymbol{q}}_h^\star )\,\|_{\Omega } & \;= \inf _{f_h\in W_h} \|\,f-f_h \|_{{\mathcal{T}_h}}. \end{alignat*}$$

5.2. Postprocessing the approximate scalar variable

There are a few well-known ways to postprocess to obtain a new approximation $u^\star _h$ of enhanced accuracy.

As the first postprocessing method, we define $u^\star _{h,1}$ in $\mathcal{P}_{k+1}(K)$ satisfying

$$\begin{align} (c^{-1} \mathop{\boldsymbol{\nabla }}u_{h,1}^\star , \mathop{\boldsymbol{\nabla }}w )_K & = (f,w)_K - \langle {{\boldsymbol{q}}_h\cdot {\boldsymbol{n}},w} \rangle _{\partial K}, &&\text{ for all } w \in \mathcal{P}_{k+1}^0(K), {\tag{5.1a}}\\ m_K(u_h^\star ) &= m_K(u_h), {\tag{5.1b}} \end{align}$$

where $m_K(v)=\operatorname {meas}(K)^{-1} (v,1)_K$ and $\mathcal{P}_{k+1}^0(K) = \{ p \in \mathcal{P}_{k+1}(K): m_K(p)=0\}.$ Clearly, $u_{h,1}^\star$ satisfies a discrete Neumann problem on each element with the computed approximate solution as data. This method was introduced in 151920 in the framework of mixed methods.

A variation of this postprocessing was proposed in 912 for DG methods. It is obtained simply by substituting $\widehat{\boldsymbol{q}}_h$ in place of ${\boldsymbol{q}}_h$ on the right-hand side of 5.1a. This yields our second postprocessing alternative. The solution so obtained is denoted by $u^\star _{h,2}$.

As a third alternative, we define $u^\star _{h,3}$ obtained by following 1920. Let $\mathcal{W}_{k+1}(K)$ denote the $L^2(K)$-orthogonal complement of $\mathcal{P}_{k-1}(K)$ in $\mathcal{P}_{k+1}(K)$. The solution $u_{h,3}^\star$ is of the form $u_h+\eta _h$ where $\eta _h$ is the unique function in $\mathcal{W}_{k+1}(K)$ satisfying

$$\begin{alignat}{2} (\mathop{\boldsymbol{\nabla }}\eta _h , \mathop{\boldsymbol{\nabla }}w )_K =&-(\mathop{\boldsymbol{\nabla }}u_h+c\,{\boldsymbol{q}}_h, \mathop{\boldsymbol{\nabla }}w )_K &&\qquad \text{ for all } w \in \mathcal{W}_{k+1}(K). {\tag{5.2}} \end{alignat}$$

As in the case of the flux postprocessing, if we are using an $L^2(K)$-orthogonal hierarchical basis to find $\eta _h$, we need only invert a symmetric, positive definite matrix whose order is the dimension of $\mathcal{W}_{k+1}(K)$. Note also that to evaluate the right-hand side of 5.2, we need only use $n-1$-dimensional quadratures, as

$$\begin{equation*} -(\mathop{\boldsymbol{\nabla }}u_h+c\,{\boldsymbol{q}}_h, \mathop{\boldsymbol{\nabla }}w )_{\mathcal{T}_h}=\langle {\widehat{u}_h-u_h,\mathop{\boldsymbol{\nabla }}w\cdot {\boldsymbol{n}}} \rangle _{{\partial \mathcal{T}_h}\setminus \partial \Omega }+\langle {g-u_h,\mathop{\boldsymbol{\nabla }}w\cdot {\boldsymbol{n}}} \rangle _{\partial \Omega }. \end{equation*}$$

This follows by setting ${\boldsymbol{r}}=\mathop{\boldsymbol{\nabla }}w$ in 1.3a, the first equation of the HDG method.

All the above postprocessed solutions converge at a higher rate than $u_h$ (whenever $k\ge 1$) as stated in the next theorem. It can be proved using the superconvergence estimate of Theorem 4.1 in a standard way (see 20) so we omit it.

Theorem 5.2.

Under the same assumption as Theorem 4.1, the result $u_{h,i}^\star$ of any of the three above mentioned postprocessings satisfy

$$\begin{equation*} \|\,u-u^\star _{h,i} \,\|_{\mathcal{T}_h}\le (C_{2,\tau }+\delta _{i2}C_{1,\tau })\,h^{\min \{k,1\}}\|\,\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}\,\|_{\mathcal{T}_h}+C\,h^{\ell +2}\,\|\,u\,\|_{H^{\ell +2}({\mathcal{T}_h})}, \end{equation*}$$

for any $k\ge 0$ and $i=1,2,3,$ and any $\ell \in [0,k]$.

6. Concluding remarks

We have presented a technique for error analysis of HDG methods that is remarkable for its brevity, especially in comparison with previous DG analyses 2912. We achieved this through the use of a new projection $\varPi _h$. While brevity and elegance is traditionally achieved in the analysis of mixed methods, like RT and BDM methods, via the use of projections with commutativity properties, in the case of DG methods, commutativity seems not to be of paramount importance. Rather, what seems important is a projection tailored to fit the structure of the DG method, such as our $\varPi _h$.

Because $\varPi _h$ is adapted to the structure of our numerical traces, we found it easy to estimate the projection of the errors. To summarize, we proved,

$$\begin{alignat*}{1} \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}_h \, \|_c & \le \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\|_c, \\ \| P_{\scriptscriptstyle {M}}({\boldsymbol{q}}\cdot {\boldsymbol{n}}) - \widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}\|_h & \le C_{1,\tau }\| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\|_c, \\ \| \varPi _{\!\scriptscriptstyle {W}}u - u_h \|_{\mathcal{T}_h}& \le C_{2,\tau } \, h^{\min \{k,1\}}\, \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\,\|_{\mathcal{T}_h}, \\ \| P_{\scriptscriptstyle {M}}u - \widehat{u}_h \|_h\;\, & \le C_{2,\tau }\, h \, \| \boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}- {\boldsymbol{q}}\,\|_{\mathcal{T}_h}, \end{alignat*}$$

where the first three estimates hold for all $k\ge 0$ and the last for $k\ge 1$. Thus, by the approximation properties of the projection $\varPi _h$ (Theorem 2.1), if the penalty function $\tau$ is such that $\tau _K^{\mathrm{max}}$ is of order one on each $K\in {\mathcal{T}_h}$, we obtain the optimal order of convergence of $k+1$ for the approximate flux and its numerical trace. Of course, by the triangle inequality, the above estimates imply that the error of these variables converges to zero at the optimal order. If $k\ge 1$, the projection of the errors for $u$ and its trace superconverge at order $k+2$. This can be exploited to get locally postprocessed solutions of enhanced accuracy.

To end, note that the above estimates imply that

$$\begin{equation*} \| (\widehat{\boldsymbol{q}}_h-{\boldsymbol{q}}_h)\cdot {\boldsymbol{n}}\|_h = O(h^{k+1}). \end{equation*}$$

According to the main result in 12, such an inequality implies optimal order of convergence for the numerical flux ${\boldsymbol{q}}_h$ and its postprocessing ${\boldsymbol{q}}_h^\star$ for $k\ge 0$. It also implies, for $k\ge 1$, the superconvergence of the orthogonal projection of the error in $u$ into $\mathcal{P}_{k-1}(K)$, and furthermore the superconvergence of $\widehat{u}_h$. This is in perfect agreement with our results.

The extension of our approach to other equations of practical interest appearing in, for example, fluid flow and solid mechanics, can prove to be useful not only to analyze already existing HDG methods but also to devise new ones. This constitutes the objective of ongoing research.