By Bernardo Cockburn, Jayadeep Gopalakrishnan, and Francisco-Javier Sayas
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 Reference 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 Reference 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 Reference 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)\}, \cssId{texmlid28}{\tag{1.2a}}\\ W_h &= \{ w: \text{ for all mesh elements }K, w|_K \in W(K)\}, \cssId{texmlid29}{\tag{1.2b}}\\ M_h &= \{ \mu : \text{ for all mesh faces } F, \;\mu |_F \in M(F)\}, \cssId{texmlid30}{\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
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 Reference 11. Above and throughout, we use the notation
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 Reference 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 Equation 1.3a–Equation 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 Reference 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 Reference 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 Reference 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 Reference 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 Reference 8.
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 Reference 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 Reference 18 and the Brezzi-Douglas-Marini (BDM) method Reference 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 Reference 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 Reference 9.
In this paper, we provide a new approach to the error analysis of HDG methods. It is an alternative to the techniques in Reference 12 for general HDG methods, and those in Reference 9 for the particular method treated therein. We recover all the results of Reference 9Reference 12 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
where $P_M$ is the $L^2$-projection into $M_h$; cf. the expression for $\widehat{\boldsymbol{q}}_h\cdot {\boldsymbol{n}}$ from Equation 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
Besides the projections used to carry out the analysis of the DG methods in Reference 8, Reference 6 and Reference 13, we have the projections used in Reference 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}$,
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 Reference 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 Reference 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:
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
for all faces $F$ of the simplex $K$. If $k=0$, then Equation 2.1a and Equation 2.1b are vacuous and $\varPi _h$ is defined solely by Equation 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 Equation 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 Equation 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
That the left-hand sides of Equation 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 Equation 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
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 Reference 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
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,
where we have used the definition of $\varPi _h$. Note that this proves the identity Equation 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 Equation 1.3d, both terms above are zero. This completes the proof.
Taking ${\boldsymbol{r}}:=\boldsymbol{\varepsilon }_h^{\;q}$ in Equation 3.1a, $w=\varepsilon _h^u$ in Equation 3.1b, $\mu =-\widehat{\boldsymbol{\varepsilon }}_h\cdot {\boldsymbol{n}}$ in Equation 3.1c and $\mu =-\varepsilon _h^{\widehat{u}}$ in Equation 3.1d, and adding the resulting four equations, we get
by the definition of $\widehat{\boldsymbol{\varepsilon }}_h$, namely Equation 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 Equation 1.3 and let the exact solution ${\boldsymbol{q}}, u$ be in the domain of $\varPi _h$. Then, for $k\ge 0$,
where we used an inverse inequality and the fact that $c^{-1}$ is uniformly bounded. Then Equation 3.3 follows from Equation 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 Reference 3Reference 17 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 )$:
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 Reference 16. Recall that we have been tacitly assuming that $({\boldsymbol{q}},u)$ is in the domain of $\varPi _h$. By Equation 4.2, $(\boldsymbol{\varPhi },\varPsi )$ is also regular enough to apply $\varPi _h$, so we have the following lemma.
by the continuity of $\boldsymbol{\varPhi }\cdot {\boldsymbol{n}}$ and the fact that $\varepsilon _h^{\widehat{u}}=0$ on $\partial \Omega$ by Equation 3.1c. Then
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
by the regularity estimate Equation 4.2, and the first estimate follows.
The second estimate follows from the first from the same local argument used in Reference 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 Equation 3.1a, and applying an inverse inequality, we find that
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 Reference 18, as used in the framework of Darcy flows Reference 4, or for the Navier-Stokes equations Reference 14, or for the postprocessing of superconvergent DG methods for second-order elliptic problems Reference 12. We follow Reference 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 Reference 12.
Theorem 5.1.
For any $k\ge 0$, we have that ${\boldsymbol{q}}_h^\star \in H(\text{div},\Omega )$. Moreover,
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 Reference 15Reference 19Reference 20 in the framework of mixed methods.
A variation of this postprocessing was proposed in Reference 9Reference 12 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 Equation 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 Reference 19Reference 20. 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). \cssId{texmlid23}{\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 Equation 5.2, we need only use $n-1$-dimensional quadratures, as
This follows by setting ${\boldsymbol{r}}=\mathop{\boldsymbol{\nabla }}w$ in Equation 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 Reference 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
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 Reference 2Reference 9Reference 12. 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,
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.
According to the main result in Reference 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.
We begin by observing that once we prove the approximation estimates of the theorem, then the unisolvency of the equations defining the projection follows as a corollary. This is because the number of equations and unknowns in Equation 2.1, namely
respectively, coincide, i.e. Equation 2.1 is a square linear system. Hence setting $u=0$ and ${\boldsymbol{q}}=0$ in the approximation estimates we find that the projection must vanish.
In view of this, in the remainder of this section, we develop estimates for any$\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ and $\varPi _{\!\scriptscriptstyle {W}}u$ satisfying Equation 2.1 without assuming uniqueness a priori (although it will follow a posteriori). But first, we begin with two auxiliary lemmas.
A.1. Two estimates involving orthogonal polynomials
Let $F$ be any face of a simplex $K$. The trace map
$$\begin{equation*} \mathop{\gamma _F}: \mathcal{P}_k^\perp (K) \longmapsto \mathcal{P}_k(F) \quad \text{ defined by }\quad \mathop{\gamma _F}(p) = p|_F \end{equation*}$$
is a bijection. Moreover,
$$\begin{align*} \| p \|_K & \le C h_K^{1/2} \| p \|_F &&\text{ for all } p\in \mathcal{P}^\perp _k(K). \end{align*}$$Proof.
We first prove that $\mathop{\gamma _F}$ is injective. Suppose $\mathop{\gamma _F}(p)=0$ for some $p\in \mathcal{P}_k^\perp (K)$. Then we can write $p = \lambda _F\,q$, where $\lambda _F$ denotes the barycentric coordinate function of $K$ that vanishes on $F$ and $q$ is some function in $\mathcal{P}_{k-1}(K)$. But since $(p,w)_K=0$ for all $w \in \mathcal{P}_{k-1}(K)$, we then have $(\lambda _F\,q, q )_K =0,$ so $q=0$ and hence $p = 0$. The surjectivity of $\mathop{\gamma _F}$ now follows by counting dimensions and using the injectivity.
Finally, the estimate of the lemma follows from the injectivity and a standard scaling argument.
■Lemma A.2.
Let $\eta$ be a nonnegative function on $\partial K$, constant on each face of $K$, and such that $\eta ^{\max }:=\max \eta > 0$. Let $p\in \mathcal{P}^\perp _k(K)$ satisfy the equation
A.2. Decoupling the projection component $\varPi _{\!\scriptscriptstyle {W}} u$
Now we characterize the second (scalar) component of the projection $\varPi _h ({\boldsymbol{q}}, u)\equiv (\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}, \varPi _{\!\scriptscriptstyle {W}}u)$, namely $\varPi _{\!\scriptscriptstyle {W}}u$, independently of the first, and prove its approximation properties.
Proposition A.1.
On each element $K\in \mathcal{T}_h$, the component $\varPi _{\!\scriptscriptstyle {W}}u$ satisfies
$$\begin{align} (\varPi _{\!\scriptscriptstyle {W}}u, v )_K &= ( u, v )_K && \quad \text{for all }v\in \mathcal{P}_{k-1}(K), \cssId{texmlid24}{\tag{A.1a}}\\ \langle {\tau \varPi _{\!\scriptscriptstyle {W}}u, w } \rangle _{\partial K} &= (\nabla \cdot {\boldsymbol{q}},w)_{K} + \langle {\tau u, w } \rangle _{\partial K} && \quad \text{for all }w \in \mathcal{P}^\perp _k(K). \cssId{texmlid25}{\tag{A.1b}} \end{align}$$Proof.
Proposition A.1 permits comparison with the SFH method. Suppose $\tau$ is selected as in the SFH method and suppose $\nabla \cdot {\boldsymbol{q}}\in \mathcal{P}_{k-1}(K)$. Then the system Equation A.1 becomes
$$\begin{alignat*}{2} (\varPi _{\!\scriptscriptstyle {W}}u, v )_K &= ( u, v)_K && \quad \text{for all }v \in \mathcal{P}_{k-1}(K), \\ \langle {\tau \varPi _{\!\scriptscriptstyle {W}}u, \mu } \rangle _{F_K} &= \langle {\tau u, \mu } \rangle _{F_K} && \quad \text{for all }\mu \in \mathcal{P}_k(F_K). \end{alignat*}$$
Note that to obtain the last equation, we used the surjectivity of $\mathop{\gamma _F}$ (Lemma A.1). Thus, in this case, $\varPi _{\!\scriptscriptstyle {W}}$ coincides with the projection used in the analysis of the SFH method Reference 9 (denoted there by $\mathbb{P}$).
We are now ready to obtain the estimate of $\varPi _{\!\scriptscriptstyle {W}}u-u$ in Theorem 2.1.
Proposition A.2.
Suppose the assumptions on $\tau$ in Theorem 2.1 hold. Then,
for $\ell _u, \ell _{\boldsymbol{q}}$ in $[0,k]$ (and consequently, $\varPi _{\!\scriptscriptstyle {W}}u$ is uniquely determined by Equation A.1).
Proof.
To prove the result, we set $\delta ^u:=\varPi _{\!\scriptscriptstyle {W}}u-u_k$, where $u_k$ is the $L^2(K)$-orthogonal projection of $u$ into $\mathcal{P}_k(K)$, and note that
where $b_{\boldsymbol{q}}(w):=(\nabla \cdot {\boldsymbol{q}},w)_{K}$ and $b_u(w):=\langle {\tau \,(u-u_k), w} \rangle _{\partial K}$. By Lemma A.2 with $\eta :=\tau$,$p:=\delta ^u$ and $b=b_{\boldsymbol{q}}+b_u$, this implies that
where $(\nabla \cdot {\boldsymbol{q}})_{k-1}$ is the $L^2(K)$-projection of $\nabla \cdot {\boldsymbol{q}}$ into $\mathcal{P}_{k-1}(K)$ when $k\ge 1$. If $k=0$, we set $(\nabla \cdot {\boldsymbol{q}})_{-1}\equiv 0$. Hence,
for any $\ell _u$ in $[0,k]$. This completes the proof.
■
A.3. Properties of the flux component $P_{\scriptscriptstyle {V}}\boldsymbol{q}$
To study the flux component, we recall another projection $\boldsymbol{B}_{\!\scriptscriptstyle {V}}$ introduced and studied in Reference 8 and later used in Reference 9Reference 12. Let $F^*$ be a face of $K$ at which $\tau |_{\partial K}$ is a maximum. For any function ${\boldsymbol{q}}$ in the domain of $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}$, the restriction of $\boldsymbol{B}_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ to $K$ is defined to be the unique element of $\boldsymbol{\mathcal{P}}^k(K)$ satisfying
for all faces $F$ of $K$ different from $F^*$. We use it to prove the following lemma. Clearly, the proof of Theorem 2.1 would be complete once we prove the lemma.
Proposition A.3.
Suppose the assumptions on $\tau$ in Theorem 2.1 hold. Then,
A basis for $\mathbb{R}^n$ is furnished by the set of unit normals ${\boldsymbol{n}}_F$ for the $n$ faces $F\ne F^*$ of $K$. Letting {$\widetilde{{\boldsymbol{n}}}_F: F\ne F^*\}$ denote its dual basis, that is, $\widetilde{{\boldsymbol{n}}}_F\cdot {\boldsymbol{n}}_{F'}=\delta _{FF'}$, we can write $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}}=\sum _{F\neq F^*} ((\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}})\cdot {\boldsymbol{n}}_F) \widetilde{{\boldsymbol{n}}}_F$. Hence, it is enough to estimate each of the functions $(\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-{\boldsymbol{q}})\cdot {\boldsymbol{n}}_F$. Moreover, since we have
by the approximation properties of $\boldsymbol{B}_{\!\scriptscriptstyle {V}}$ established in Reference 8, it is enough to estimate the function $\delta ^{\boldsymbol{q}}_F:= (\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}-\boldsymbol{B}_{\!\scriptscriptstyle {V}}{\boldsymbol{q}})\cdot {\boldsymbol{n}}_F$.
where $\eta$ is the characteristic function of $F$ and $b(w):=\langle {\tau \eta \,(u-\varPi _{\!\scriptscriptstyle {W}}u), w} \rangle _{\partial K}$. By Lemma A.2, this implies that
for $\ell _u, \ell _{\boldsymbol{q}}$ in $[0,k]$, by the approximation properties of the $L^2$-projection and the estimates of Proposition A.2. This proves the required estimate for $\delta ^{\boldsymbol{q}}_F$ and completes the proof.
■
We have thus finished the proof of Theorem 2.1. Before closing, let us state one more result that characterizes $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ and compare with the case of the SFH method.
Proposition A.4.
Suppose that the assumptions on $\tau$ of Theorem 2.1 hold. Then, on each simplex $K\in \mathcal{T}_h$,$\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ is the only element of $\boldsymbol{\mathcal{P}}_k(K)$ such that
and all faces $F$ of $K$ except one arbitrarily chosen.
Note that for the SFH method, taking all faces $F$ where $\tau =0$, we immediately see that $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}$ coincides with the projection in the analysis of that method in Reference 9. It was denoted therein by $\boldsymbol{\Pi }_{\scriptscriptstyle V}$ and is nothing but the projection $\boldsymbol{B}_{\!\scriptscriptstyle {V}}$ given by Equation A.2.
Proof.
The component $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ obviously satisfies Equation A.4 as these equations are identical to Equation 2.1a and Equation 2.1c. Since the system Equation A.4 is square, we need only show uniqueness. So, let us set the right-hand side equal to zero. For any face $F$, consider the function $\delta _F:=\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}\cdot {\boldsymbol{n}}_F$. Then by equation Equation A.4a with ${\boldsymbol{v}}:={\boldsymbol{n}}_F\,v$ with $v\in \mathcal{P}_{k-1}(K)$, equation Equation A.4b, and Lemma A.1, we get that $\delta _F=0$. Since this can be done for all faces but one, we readily obtain that $\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}=\boldsymbol{0}$. This completes the proof.
■
Acknowledgements
The authors would like to thank the referees for constructive criticism leading to a better version of this paper.
Suppose $k\ge 0$,$\tau |_{\partial K}$ is nonnegative and $\tau ^{\max }_K:=\max \tau |_{\partial K} >0$. Then the system Equation 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
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.
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
Let ${\boldsymbol{q}}_h,\; u_h,$ and $\widehat{u}_h$ solve the HDG equations Equation 1.3 and let the exact solution ${\boldsymbol{q}}, u$ be in the domain of $\varPi _h$. Then, for $k\ge 0$,
$$\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). \cssId{texmlid23}{\tag{5.2}} \end{alignat}$$
Lemma A.1.
Let $F$ be any face of a simplex $K$. The trace map
$$\begin{equation*} \mathop{\gamma _F}: \mathcal{P}_k^\perp (K) \longmapsto \mathcal{P}_k(F) \quad \text{ defined by }\quad \mathop{\gamma _F}(p) = p|_F \end{equation*}$$
is a bijection. Moreover,
$$\begin{align*} \| p \|_K & \le C h_K^{1/2} \| p \|_F &&\text{ for all } p\in \mathcal{P}^\perp _k(K). \end{align*}$$
Lemma A.2.
Let $\eta$ be a nonnegative function on $\partial K$, constant on each face of $K$, and such that $\eta ^{\max }:=\max \eta > 0$. Let $p\in \mathcal{P}^\perp _k(K)$ satisfy the equation
Suppose that the assumptions on $\tau$ of Theorem 2.1 hold. Then, on each simplex $K\in \mathcal{T}_h$,$\boldsymbol{\varPi }_{\!\scriptscriptstyle {V}}{\boldsymbol{q}}$ is the only element of $\boldsymbol{\mathcal{P}}_k(K)$ such that
and all faces $F$ of $K$ except one arbitrarily chosen.
References
[1]
D. N. Arnold and F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO Modél. Math. Anal. Numér. 19 (1985), 7–32. MR 813687 (87g:65126)
Reference [2]
D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2002), 1749–1779. MR 1885715 (2002k:65183)
Reference [3]
J.-P. Aubin, Behavior of the error of the approximate solutions of boundary value problems for linear elliptic operators by Galerkin’s and finite difference methods, Ann. Scuola Norm. Sup. Pisa (3) 21 (1967), 599–637. MR 0233068 (38:1391)
Reference [4]
P. Bastian and B. Rivière, Superconvergence and $H(\text{div})$ projection for discontinuous Galerkin methods, Internat. J. Numer. Methods Fluids 42 (2003), 1043–1057. MR 1991232 (2004f:65177)
Reference [5]
F. Brezzi, J. Douglas, Jr., and L. D. Marini, Two families of mixed finite elements for second order elliptic problems, Numer. Math., 47 (1985), pp. 217–235. MR 799685 (87g:65133)
Reference [6]
F. Celiker and B. Cockburn, Superconvergence of the numerical traces of discontinuous Galerkin and hybridized mixed methods for convection-diffusion problems in one space dimension, Math. Comp. 76 (2007), 67–96. MR 2261012 (2008e:65225)
Reference [7]
P. Ciarlet, The finite element method for elliptic problems, North-Holland, Amsterdam, 1978. MR 0520174 (58:25001)
Reference [8]
B. Cockburn and B. Dong, An analysis of the minimal dissipation local discontinuous Galerkin method for convection-diffusion problems, J. Sci. Comput. 32 (2007), 233–262. MR 2320571 (2008e:65345)
Reference [9]
B. Cockburn, B. Dong, and J. Guzmán, A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems, Math. Comp. 77 (2008), 1887–1916. MR 2429868 (2009d:65166)
[10]
B. Cockburn, J. Gopalakrishnan, Error analysis of variable degree mixed methods for elliptic problems via hybridization, Math. Comp. 74 (2005), 1653–1677. MR 2164091 (2006e:65215)
Reference [11]
B. Cockburn, J. Gopalakrishnan, and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and conforming Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), 1319–1365. MR 2485455
Reference [12]
B. Cockburn, J. Guzmán, and H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comp. 78 (2009), 1–24. MR 2448694 (2009i:65213)
Reference [13]
B. Cockburn, G. Kanschat, I. Perugia, and D. Schötzau, Superconvergence of the local discontinuous Galerkin method for elliptic problems on Cartesian grids, SIAM J. Numer. Anal. 39 (2001), 264–285. MR 1860725 (2002g:65140)
Reference [14]
B. Cockburn, G. Kanschat, and D. Schötzau, A locally conservative LDG method for the incompressible Navier-Stokes equations, Math. Comp. 74 (2005), 1067–1095. MR 2136994 (2006a:65157)
Reference [15]
L. Gastaldi and R. H. Nochetto, Sharp maximum norm error estimates for general mixed finite element approximations to second order elliptic equations, RAIRO Modél. Math. Anal. Numér. 23 (1989), 103–128. MR 1015921 (91b:65125)
Reference [16]
P. Grisvard, Elliptic Problems in Nonsmooth Domains, no. 24 in Monographs and Studies in Mathematics, Pitman Advanced Publishing Program, Marshfield, Massachusetts, 1985. MR 775683 (86m:35044)
Reference [17]
J. Nitsche, Ein Kriterium für die Quasi-Optimalität des Ritzschen Verfahrens, Numer. Math., 11 (1968), pp. 346–348. MR 0233502 (38:1823)
Reference [18]
P.-A. Raviart and J. M. Thomas, A mixed finite element method for second order elliptic problems, Mathematical Aspects of Finite Element Method (I. Galligani and E. Magenes, eds.), Lecture Notes in Math. 606, Springer-Verlag, New York, 1977, pp. 292–315. MR 0483555 (58:3547)
Reference [19]
R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math 53 (1988), 513–538. MR 954768 (89h:65192)
Reference [20]
R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO Modél. Math. Anal. Numér. 25 (1991), 151–167. MR 1086845 (92a:65303)
The first author was supported in part by the National Science Foundation (Grant DMS-0712955) and by the University of Minnesota Supercomputing Institute.
The second author was supported in part by the National Science Foundation under grants DMS-0713833 and SCREMS-0619080.
The third author was partially supported by MEC/FEDER Project MTM2007–63204, Gobierno de Aragón (Grupo PDIE) and was a Visiting Professor of the School of Mathematics, University of Minnesota, during the development of this work.
Journal Information
Mathematics of Computation, Volume 79, Issue 271, ISSN 1088-6842, published by the American Mathematical Society, Providence, Rhode Island.
Publication History
This article was received on , revised on , and published on .
Copyright Information
Copyright 2010 American Mathematical Society; reverts to public domain 28 years from publication
Note. To explore an equation, focus it (e.g., by clicking on it) and use the arrow keys to navigate its structure. Screenreader users should be advised that enabling speech synthesis will lead to duplicate aural rendering.