$hp$-discontinuous Galerkin Methods for the Helmholtz Equation with Large Wave Number

TL;DR

This paper introduces stable $hp$-discontinuous Galerkin methods for solving the Helmholtz equation with large wave numbers, providing error estimates and stability analysis that depend explicitly on mesh size, polynomial degree, and wave number.

Contribution

The paper develops novel $hp$-DG methods with complex penalty parameters ensuring stability and derives explicit error and stability estimates for high wave number Helmholtz problems.

Findings

Methods are absolutely stable and well-posed.

Error estimates are sub-optimal without mesh constraints.

Optimal order error estimates are achieved under specific mesh conditions.

Abstract

This paper develops some interior penalty $hp$-discontinuous Galerkin ($hp$-DG) methods for the Helmholtz equation in two and three dimensions. The proposed $hp$-DG methods are defined using a sesquilinear form which is not only mesh-dependent but also degree-dependent. In addition, the sesquilinear form contains penalty terms which not only penalize the jumps of the function values across the element edges but also the jumps of the first order tangential derivatives as well as jumps of all normal derivatives up to order $p$. Furthermore, to ensure the stability, the penalty parameters are taken as complex numbers with positive imaginary parts. It is proved that the proposed $hp$-discontinuous Galerkin methods are absolutely stable (hence, well-posed). For each fixed wave number $k$, sub-optimal order error estimates in the broken $H^1$-norm and the $L^2$-norm are derived without any mesh constraint. The error estimates and the stability estimates are improved to optimal order under the mesh condition $k^3h^2p^{-1}\le C_0$ by utilizing these stability and error estimates and using a stability-error iterative procedure To overcome the difficulty caused by strong indefiniteness of the Helmholtz problems in the stability analysis for numerical solutions, our main ideas for stability analysis are to make use of a local version of the Rellich identity (for the Laplacian) and to mimic the stability analysis for the PDE solutions given in \cite{cummings00,Cummings_Feng06,hetmaniuk07}, which enable us to derive stability estimates and error bounds with explicit dependence on the mesh size $h$, the polynomial degree $p$, the wave number $k$, as well as all the penalty parameters for the numerical solutions.

Equations

$$-\Delta u-k^{2}u$$

$$\frac{\partial u}{\partial n_{\Omega }}+\mathbf{i}ku$$

$$u$$

$$H_{{\Gamma }_D}^1(\Omega ):=\left\{u\in H^1(\Omega );u=0\text{ on }{\Gamma }_D\right\}.$$

$$(x-x_Q)\cdot n_Q\ge c_Q\forall x\in \partial Q.$$

$$\displaystyle\|u\|_{H^{j}(\Omega)}\lesssim\Bigl{(}\frac{1}{k}+k^{j-1}\Bigr{)}\bigl{(}\|f\|_{L^{2}(\Omega)}+\|g\|_{L^{2}(\Gamma_{R})}\bigr{)}$$

$$\Vert u{\Vert }_{H_{\text{mbox}loc}^j(\Omega )}$$

$${\mathcal{E}}_h^I$$

$${\mathcal{E}}_h^R$$

$${\mathcal{E}}_h^D$$

$${\mathcal{E}}_h^{RD}$$

$${\mathcal{E}}_h^{ID}$$

$${\mathcal{E}}_h$$

$$[v]|_{e}:=\left\{\begin{array}[]{ll}v|_{K}-v|_{K^{\prime}},&\quad\mbox{if the global label of $K$ is larger},\\ v|_{K^{\prime}}-v|_{K},&\quad\mbox{if the global label of $K^{\prime}$ is larger}.\end{array}\right.$$

$$\{v\}|_{e}:=\frac{1}{2}\bigl{(}v|_{K}+v|_{K^{\prime}}\bigr{)}\qquad\mbox{if }e=\partial K\cap\partial K^{\prime}.$$

$$E^q:=\prod _{K\in {\mathcal{T}}_h}{H}^{q+1}(K),$$

$$a_{h}^{q}(u,v):=b_{h}(u,v)+{\rm\mathbf{i}}\Bigl{(}L_{1}(u,v)+\sum_{j=0}^{q}J_{j}(u,v)\Bigr{)}\qquad\forall u,v\in E^{q},$$

$$b_h(u,v)$$

$$+\sigma {⟨\left[u\right],\left\{\frac{\partial v}{\partial n_e}\right\}⟩}_e),$$

$$L_1(u,v)$$

$$J_0(u,v)$$

$$J_j(u,v)$$

$$\displaystyle\bigl{|}v\bigr{|}_{1,h}$$

$${\Vert v\Vert }_{1,h,q}$$

$${+\sum _{j=1}^q\sum _{e\in {\mathcal{E}}_h^I}{\gamma }_{j,e}{\left(\frac{h_e}{p}\right)}^{2j-1}{\Vert \left[\frac{{\partial }^{j}v}{\partial n_e^j}\right]\Vert }_{L^2(e)}^2)}^{\frac{1}{2}},$$

$$\Vert |v\Vert {|}_{1,h,q}$$

$$\mathrm{Re}{a}_h^q(v,v)$$

$$\mathrm{Im}{a}_h^q(v,v)$$

$$a_h^q(u,v)-k^2(u,v)+\mathbf{i}k⟨u,v{⟩}_{{\Gamma }_R}=(f,v)+⟨g,v{⟩}_{{\Gamma }_R}$$

$$V_h^p:=\prod _{K\in {\mathcal{T}}_h}{\mathcal{P}}_p(K).$$

$$a_h^q(u_h^q,v_h)-k^2(u_h^q,v_h)+\mathbf{i}k⟨u_h^q,v_h{⟩}_{{\Gamma }_R}=(f,v_h)+⟨g,v_h{⟩}_{{\Gamma }_R}\forall v_h\in V_h^p.$$

$$\displaystyle d\left\|v\right\|_{L^{2}(K)}^{2}+2\operatorname{{Re}}(v,\alpha\cdot\nabla v)_{K}=\int_{\partial K}\alpha\cdot n_{K}\bigl{|}v\bigr{|}^{2},$$

$$\displaystyle(d-2)\left\|\nabla v\right\|_{L^{2}(K)}^{2}+2\operatorname{{Re}}\big{(}\nabla v,\nabla(\alpha\cdot\nabla v)\big{)}_{K}=\int_{\partial K}\alpha\cdot n_{K}\bigl{|}\nabla v\bigr{|}^{2},$$

$${⟨\left\{\frac{\partial v}{\partial n_e}\right\},\left[\alpha \cdot \nabla v\right]⟩}_e-{⟨\alpha \cdot n_e\left\{\nabla v\right\},\left[\nabla v\right]⟩}_e$$

$$=\sum _{\ell =1}^{d-1}{\int }_e\left(\alpha \cdot {\tau }_e^{\ell }\left\{\frac{\partial v}{\partial n_e}\right\}-\alpha \cdot n_e\left\{\frac{\partial v}{\partial {\tau }_e^{\ell }}\right\}\right)\left[\frac{\partial \overline{v}}{\partial {\tau }_e^{\ell }}\right],$$

$${\Vert z\Vert }_{L^2(\partial K)}$$

$${\Vert \nabla z\Vert }_{L^2(K)}$$

$$a_h^q(u_h^q,u_h^q)-k^2{\Vert u_h^q\Vert }_{L^2(\Omega )}^2+\mathbf{i}k{\Vert u_h^q\Vert }_{L^2({\Gamma }_R)}^2=(f,u_h^q)+⟨g,u_h^q{⟩}_{{\Gamma }_R}.$$

$$\displaystyle\bigl{|}u_{h}^{q}\bigr{|}_{1,h}^{2}-2\operatorname{{Re}}\sum_{e\in\mathcal{E}_{h}^{ID}}\left\langle\left\{\frac{\partial u_{h}^{q}}{\partial n_{e}}\right\},\left[u_{h}^{q}\right]\right\rangle_{e}-k^{2}\left\|u_{h}^{q}\right\|_{L^{2}(\Omega)}^{2}$$

$$\displaystyle\hskip 166.2212pt\leq\bigl{|}(f,u_{h}^{q})+\langle g,u_{h}^{q}\rangle_{\Gamma_{R}}\bigr{|},$$

Full text

The work of the first author was partially supported by the NSF grants DMS-0410266 and DMS-0710831. The work of the second author was partially supported by the national basic research program of China under grant 2005CB321701 and by the program for the new century outstanding talents in universities of China. Part of this joint work was carried out during the second author’s recent visit of the University of Tennessee, the author would like to thank the host institution for its hospitality and financial support of the visit.

1. Introduction

This is the second installment in a series (cf. [27]) which devotes to developing and analyzing novel interior penalty discontinuous Galerkin (IPDG) methods for the following Helmholtz problem with large wave number:

[TABLE]

where $k\in\mathbf{R}$, called wave number, is a (large) positive number, $D\subset\Omega_{1}\subset\mathbf{R}^{d},\,d=2,3$, $D$ is known as a scatterer and is assumed to be a bounded Lipschitz domain, $\Omega_{1}$, which is assumed to be a polygonal/polyhedral domain and often taken as a $d$-rectangle in applications, defines the size of the computational domain. Note that $\partial\Omega=\Gamma_{R}\cup\Gamma_{D}$. $n_{\Omega}$ denotes the unit outward normal to $\Omega$. ${\rm\mathbf{i}}:=\sqrt{-1}$ denotes the imaginary unit. Condition (1.2) with $g=0$ is known as the first order absorbing boundary condition (cf. [23]), which is used to minimize the reflection at the boundary $\Gamma_{R}$ and to limit the computation of the original scattering problem just on the finite domain $\Omega$. Boundary condition (1.3) implies that the scatterer is sound-soft. We note that the case $D=\emptyset$ also arises in applications either as a consequence of frequency domain treatment of waves or when time-harmonic solutions of the scalar wave equation are sought (cf. [22]).

In [27] we proposed and analyzed some IPDG methods for problem (1.1)–(1.3) using piecewise linear polynomial trial and test functions. It was proved that the proposed methods are unconditionally (with respect to mesh size $h$) stable and well-posed for all wave number $k>0$. Optimal order error estimates were established showing explicit dependence of the error bounds on $h$, $k$ and all penalty parameters. However, due to the existence of a pollution term, the (broken) $H^{1}$-norm error bound deteriorates as the wave number $k$ increases under the practical “rule of thumb” mesh constraint $kh\lesssim 1$. To improve the accuracy and efficiency of those IPDG methods, it is necessary to use (piecewise) high order polynomial trial and test functions partly because of the rigidity and low approximability of linear functions and partly because of the very oscillatory nature of high frequency waves. However, simply replacing the linear element by high order elements in the IPDG methods of [27] does not reduce the pollution very much, in particular, the theoretical error bounds do not change much because the analysis of [27] indeed strongly depends on the properties of linear functions.

Motivated by the above challenge and observation, the primary goal of this paper is to develop some new $hp$-interior penalty discontinuous Galerkin ($hp$-IPDG) methods which retain the advantages of the IPDG methods of [27] but improve their accuracy and stability by exploiting the efficiency and flexibility of piecewise high order polynomial functions. To the end, our key idea is to construct a sesquilinear form (as a discretization of the Laplacian) which is not only mesh-dependent (or $h$-dependent) but also degree-dependent (or $p$-dependent) by introducing penalty terms which not only penalize the jumps of the function values across the element edges but also the jumps of the first order tangential derivatives as well as jumps of all normal derivatives up to order $p$. In addition, as in [27], to ensure the stability, all penalty parameters are taken as complex numbers with positive imaginary parts. Since the Helmholtz equation with large wave number is non-Hermitian and strongly indefinite, as expected, stability estimates (or a priori estimates) for numerical solutions under practical mesh constraints is a difficult task to carry out regardless which discretization method is used. To overcome the difficulty, as in [27], the cruxes of our analysis are to establish and to make use of a local version of the Rellich identity (for the Laplacian) and to mimic the stability analysis for the PDE solutions given in [19, 20, 33]. The key idea here is to use the special test function $\nabla u_{h}\cdot(x-x_{\Omega})$ (defined element-wise) with $u_{h}$ denoting the $hp$-IPDG solution, such a test function is valid for any DG method. We remark that the same technique was successfully employed by Shen and Wang in [38] to establish the stability and error analysis for the spectral Galerkin approximation of the Helmholtz problem. We also note that although the similar techniques to those in [27, 38] are utilized in this paper to carry out the stability analysis, the analysis of this paper is more involved because the special sesquilinear form of this paper, which contains jumps of high order normal derivatives, is a lot more complicate to deal with, even they are similar conceptually.

Since the Helmholtz equation appears, in one way or another, directly or indirectly, in almost all wave-related problems arisen from many science, engineering, and industry applications, solving the Helmholtz equation, in one form or another, has always been and remains at the center of wave computation. We refer the reader to ([1, 5, 6, 10, 13, 14, 18, 19, 22, 24, 29, 32, 36, 37, 39, 46] and the references therein) for some recent developments on numerical methods, in particular, Galerkin type methods, for the Helmholtz equation. We also refer the reader to [27] for a brief review about some theoretical issues for finite element approximations (and other types of Galerkin approximations) of the Helmholtz equation.

The $hp$-finite element method ($hp$-FEM) is a modern version of the finite element method, capable of achieving exceptionally fast (exponential) convergence. It combines the flexibility of the standard finite element method and the high order accuracy of the spectral method. Consequently, the $hp$-FEM can often attain more accurate results than the standard finite element method does while using less CPU time and resources. The $hp$-FEM has undergone intensive developments both on theory and implementation in the past twenty five years. We refer the reader to the survey paper [7] and two recent monographs [41, 42] for a detailed exposition on the basic theory and advanced topics of the $hp$-FEM.

Discontinuous Galerkin (DG) methods was first proposed in 1970s, they were not popular then because they produce larger algebraic systems than standard finite element methods do. However, due to the emergence of high performance computers and fast solvers since early 1990s, especially, massively parallel computers and parallel solvers such as multilevel and domain decomposition methods, which together with advantages of DG methods has quickly attracted renewed interests in DG methods. They have been heavily developed and tested in the past fifteen years, we refer the reader to [4] and the references therein for a review of recent developments. As is well known now, DG methods have several advantages over other types of numerical methods. For example, the trial and test spaces are easy to construct, they can naturally handle inhomogeneous boundary conditions and curved boundaries; they also allow the use of highly nonuniform and unstructured meshes, and have built-in parallelism which permits coarse-grain parallelism. In addition, the fact that the mass matrices are block diagonal is an attractive feature in the context of time-dependent problems, especially if explicit time discretizations are used. Moreover, as proved in [27], DG methods are also effective and have advantages over finite element methods for the strongly indefinite Helmholtz equation, which has not been well understood before. We refer the reader to [3, 4, 9, 16, 17, 21, 26, 40, 45] and the references therein for a detailed account on DG methods for coercive elliptic and parabolic problems, and to [43, 25, 28, 35, 34] and the references therein for recent developments on $hp$-discontinuous Galerkin ($hp$-DG) methods.

The remainder of this paper is organized as follows. In Section 2, we first introduce notation and gather some preliminaries, and then formulate our $hp$-IPDG methods. Both symmetric and non-symmetric methods are constructed and various possible variants are also discussed. Section 3 devotes to the stability analysis for the $hp$-IPDG methods proposed in Section 2. It is proved that the proposed $hp$-IPDG methods are stable (hence well-posed) without any mesh constraint. In Section 4, using the stability results of Section 3 we prove that for each fixed wave number $k$, sub-optimal order (with respect to $h$ and $p$) error estimates in the broken $H^{1}$-norm and the $L^{2}$-norm are derived without any mesh constraint. Finally, using the stability estimate of Section 3, the error estimates of Section 4 and a stability-error iterative procedure we obtain some much improved (optimal order) stability and error estimates for the $hp$-IPDG solutions under the mesh condition $k^{3}h^{2}p^{-1}\leq C_{0}$ in Section 5, where $C_{0}$ is some constant independent of $k$, $h$, $p$, and the penalty parameters.

2. Formulation of $hp$-interior penalty discontinuous Galerkin methods

2.1. Notation and preliminaries

The space, norm and inner product notation used in this paper all are standard, we refer to [11, 15, 9] for their precise definitions. On the other hand, we note that all functions in this paper are complex-valued, so the familiar terminologies such “symmetric/non-symmetric” and “bilinear” are replaced respectively by terms “Hermitian/non-Hermitian” and “sesquilinear”. For a complex number $a=a_{r}+{\rm\mathbf{i}}a_{i}$ ($a_{r}$ and $a_{i}$ are real numbers), $\overline{a}:=a_{r}-{\rm\mathbf{i}}a_{i}$ denotes the complex conjugate of $a$. $(\cdot,\cdot)_{Q}$ and $\langle\cdot,\cdot\rangle_{\Sigma}$ for $\Sigma\subset\partial Q$ denote the complex $L^{2}$-inner product on $L^{2}(Q)$ and $L^{2}(\Sigma)$ spaces, respectively. $(\cdot,\cdot):=(\cdot,\cdot)_{\Omega}$ and $\langle\cdot,\cdot\rangle:=\langle\cdot,\cdot\rangle_{\partial\Omega}$. We also define

[TABLE]

Throughout the paper, $C$ is used to denote a generic positive constant which is independent of $k$, $h$, $p$, and the penalty parameters. We also use the shorthand notation $A\lesssim B$ and $B\gtrsim A$ for the inequality $A\leq CB$ and $B\geq CA$. $A\simeq B$ is for the statement $A\lesssim B$ and $B\lesssim A$.

We now give the definition of star-shaped domains.

Definition 2.1**.**

$Q\subset\mathbf{R}^{d}$* is said to be a star-shaped domain with respect to $x_{Q}\in Q$ if there exists a nonnegative constant $c_{Q}$ such that*

[TABLE]

$Q\subset\mathbf{R}^{d}$* is said to be strictly star-shaped if $c_{Q}$ is positive.*

In this paper, we assume that $\Omega_{1}$ is a strictly star-shaped domain. Recall that $\Omega_{1}$ is often taken as a $d$-rectangle in practice. We also assume that the scatterer $D$ is a star-shaped domain, without loss of the generality, with respect to the same point $x_{\Omega_{1}}$ as $\Omega_{1}$ does. This implies that $x_{\Omega_{1}}\in D\subset\Omega_{1}$. Under these assumptions, there hold following stability estimates for problem (1.1)–(1.3).

Theorem 2.1**.**

Suppose $\Omega_{1}\subset\mathbf{R}^{d}$ is a strictly star-shaped domain and $D\subset\Omega_{1}$ is a star-shaped domain. Then the solution $u$ to problem (1.1)–(1.3) satisfies

[TABLE]

for $j=0,1$ if $u\in H^{3/2+\varepsilon}(\Omega)$ for some $\varepsilon>0$. (2.2) also holds for $j=2$ if $u\in H^{2}(\Omega)$. Furthermore, there hold

[TABLE]

for $j=3,4,\cdots,q$ if $u\in H^{2}(\Omega)\cap H^{q}_{\mbox{\small loc}}(\Omega)$ for some positive integer $q\geq 3$.

Proof.

Inequality (2.2) for $j=0,1,2$ was proved in [19, 20, 33]. Inequality (2.3) follows from (2.2) and an application of the standard cutoff function technique together with an induction argument. We leave the derivation to the interested reader. ∎

2.2. Formulation of $hp$-IPDG methods

To formulate our $hp$-IPDG methods, we need to introduce some notation, most of them were already appeared in [27]. Let $\mathcal{T}_{h}$ be a family of partitions of the domain $\Omega:=\Omega_{1}\setminus D$ parameterized by $h\in(0,h_{0})$. For any $K\in\mathcal{T}_{h}$, we define $h_{K}:=\mbox{diam}(K)$. Similarly, for each edge/face $e$ of $K\in\mathcal{T}_{h}$, $h_{e}:=\mbox{diam}(e)$. We impose the following mild restrictions on the partition $\mathcal{T}_{h}$:

• (i)

The elements of $\mathcal{T}_{h}$ satisfy the minimal angle condition,

• (ii)

$\mathcal{T}_{h}$ is locally quasi-uniform, that is if two elements $K$ and $K^{\prime}$ are adjacent ( i.e., $\mbox{meas}(\partial K\cap\partial K^{\prime})>0$ ), then $h_{K}\simeq h_{K^{\prime}}$. Where $\mbox{meas}(e)$ stands for $(d-1)$-dimensional Lebesgue measure of $e$.

For convenience, we assume $\mathrm{diam}(\Omega)\simeq 1$, hence $h_{e},h_{K}\lesssim 1.$

For any two elements $K,\,K^{\prime}\in\mathcal{T}_{h}$, we call $e=\partial K\cap\partial K^{\prime}$ an interior edge/face of $\mathcal{T}_{h}$ if $\mbox{meas}(e)>0$. Note that $e$ could be portion of a side/face of the element $K$ or $K^{\prime}$ in the case of geometrically nonconforming partition. Also, for any element $K\in\mathcal{T}_{h}$, we call $e=\partial K\cap\partial\Omega$ a boundary edge/face if $\mbox{meas}(e)>0$. Then we define

[TABLE]

We also define the jump $[v]$ of $v$ on an interior edge/face $e=\partial K\cap\partial K^{\prime}$ as

[TABLE]

If $e\in\mathcal{E}_{h}^{D}$, set $[v]|_{e}=v|_{e}$. The following convention is adopted in this paper

[TABLE]

If $e\in\mathcal{E}_{h}^{RD}$, set $\{v\}|_{e}=v|_{e}$. For every $e=\partial K\cap\partial K^{\prime}\in\mathcal{E}_{h}^{I}$, let $n_{e}$ be the unit outward normal to edge/face $e$ of the element $K$ if the global label of $K$ is bigger and of the element $K^{\prime}$ if the other way around. For every $e\in\mathcal{E}_{h}^{RD}$, let $n_{e}=n_{\Omega}$ the unit outward normal to $\partial\Omega$.

Let $p\geq 1$ be a fixed integer, which will be used to denote the degree of the $hp$-IPDG methods in this paper. For each integer $0\leq q\leq p$, we define the “energy” space

[TABLE]

and the sesquilinear form $a^{q}_{h}(\cdot,\cdot)$ on $E^{q}\times E^{q}$

[TABLE]

[TABLE]

and $\sigma$ is a real number. $\gamma_{0,e},\cdots,\gamma_{q,e}>0$ and $\beta_{1,e}\geq 0$ are numbers to be specified later. $\{\tau^{\ell}_{e}\}_{\ell=1}^{d-1}$ denote an orthogonal coordinate frame on the edge/face $e\in\mathcal{E}_{h}$, $\frac{\partial u}{\partial\tau_{e}^{\ell}}:=\nabla u\cdot\tau_{e}^{\ell}$ stands for the tangential derivative of $u$ in the direction $\tau_{e}^{\ell}$, and $\frac{\partial^{j}u}{\partial n_{e}^{j}}$ denotes the $j$th order normal derivative of $u$ on $e$.

It is easy to check that $(-\Delta u,v)=a_{h}^{q}(u,v)$ for all $u\in H^{q+1}(\Omega)$ and $v\in E^{q}$. Hence, $a_{h}^{q}(\cdot,\cdot)$ is a consistent discretization for $-\Delta$. When $\sigma=1$, $a_{h}^{q}(\cdot,\cdot)$ is symmetric, that is, $a_{h}^{q}(u,v)=a_{h}^{q}(\overline{v},\overline{u})$. On the other hand, when $\sigma\neq 1$, $a_{h}^{q}(\cdot,\cdot)$ is non-symmetric. In particular, $\sigma=-1$ would correspond to the non-symmetric IPDG method studied in [40] for coercive elliptic problems. In this paper, for the ease of presentation, we only consider the case $\sigma=1$. The penalty constants in {\rm\mathbf{i}}\bigl{(}L_{1}(u,v)+J_{0}(u,v)+\cdots+J_{q}(u,v)\bigr{)} are ${\rm\mathbf{i}}\beta_{1,e},{\rm\mathbf{i}}\gamma_{0,e},\cdots,{\rm\mathbf{i}}\gamma_{q,e}$, respectively. So they are pure imaginary numbers with positive imaginary parts. It turns out that if any of them is replaced by a complex number with positive imaginary part, the ideas of the paper still apply. Here we set their real parts to be zero because the terms from real parts do not help much (and do not cause any problem either) in our analysis.

Next, we introduce the following semi-norms on the space $E^{q}$:

[TABLE]

Clearly, $\left\|\cdot\right\|_{1,h,q}$ and $\|{\hskip-0.6pt}|\cdot\|{\hskip-0.6pt}|_{1,h,q}$ are norms on $E^{q}$ if $\partial D\neq\emptyset$ but only semi-norms if $\partial D=\emptyset$.

It is easy to check that the sesquilinear form $a_{h}^{q}(\cdot,\cdot)$ satisfies: For any $v\in E^{q}$

[TABLE]

Using the sesquilinear form $a_{h}^{q}(\cdot,\cdot)$ we now introduce the following weak formulation for (1.1)–(1.2): Find $u\in E^{q}\cap H^{1}_{\Gamma_{D}}(\Omega)\cap H^{2}_{{\rm loc}}(\Omega)$ such that

[TABLE]

for any $v\in E^{q}\cap H^{1}_{\Gamma_{D}}(\Omega)\cap H^{2}_{{\rm loc}}(\Omega)$. The above formulation is consistent with the boundary value problem (1.1)–(1.2) because $a_{h}^{q}(\cdot,\cdot)$ is consistent with $-\Delta$.

For any $K\in\mathcal{T}_{h}$, let $\mathcal{P}_{p}(K)$ denote the set of all polynomials whose degrees do not exceed $p$. We define our $hp$-IPDG approximation space $V_{h}^{p}$ as

[TABLE]

Clearly, $V_{h}^{p}\subset E^{q}\subset L^{2}(\Omega)$. But $V_{h}^{p}\not\subset H^{1}(\Omega)$. We are now ready to define our $hp$-IPDG methods based on the weak formulation (2.14): For each $0\leq q\leq p$, find $u_{h}^{q}\in V_{h}^{p}$ such that

[TABLE]

Remark 2.1**.**

(a) When $p=q=1$, the above method (2.15) is exactly the scheme proposed in [27]. The $L_{1}$ term, which penalizes the jumps of the first order tangential derivatives, plays an important role for getting a better (theoretical) stability estimate in [27]. However, our analysis to be given in the next section suggests that the $L_{1}$ term plays a less pivotal role for high order IPDG methods.

(b) In fact, (2.15) defines $p+1$ different IPDG methods for $q=0,1,\cdots,p$. $q=1$ would correspond to using high order elements in the IPDG formulation proposed in [27].

(c) The idea of penalizing the jumps of normal derivatives (i.e., the $J_{1}$ term above) for second order PDEs was used early by Douglas and Dupont [21] in the context of $C^{0}$ finite element methods, by Baker [9] (with a different weighting, also see [26]) for fourth order PDEs. The idea of using multipenalties $J_{0},J_{1},\cdots,J_{p}$ with positive penalty parameters was first used by Arnold in [3] for coercive elliptic and parabolic PDEs. The use of $L_{1}$ term was first introduced in [27].

In the next two sections, we shall study the stability and error analysis for the $hp$-IPDG method (2.15). Especially, we are interested in knowing how the stability constants and error constants depend on the wave number $k$ (and mesh size $h$ and element degree $p$, of course) and on the penalty parameters, and what are the “optimal” relationship between mesh size $h$ and the wave number $k$.

3. Stability estimates

The goal of this section is to derive stability estimates (or a priori estimates) for schemes (2.15). To the end, momentarily, we assume that the solution $u_{h}^{q}$ to (2.15) exists and will revisit the existence and uniqueness issues later at the end of the section. We like to note that because its strong indefiniteness, unlike in the case of coercive elliptic and parabolic problems (cf. [3, 4, 9, 21, 26, 40, 45]), the well-posedness of scheme (2.15) is difficult to prove under practical mesh constraints.

To derive stability estimates for scheme (2.15), our approach is to mimic the stability analysis for the Helmholtz problem (1.1)–(1.2) given in [19, 20, 33]. The key ingredients of our analysis are to use a special test function $v_{h}=\alpha\cdot\nabla u_{h}^{q}$ (defined element-wise) with $\alpha(x):=x-x_{\Omega_{1}}$ in (2.15) and to use the Rellich identity (cf. [20] and below) on each element. Due to existence of multiple penalty terms in $a_{h}^{p}(\cdot,\cdot)$, which do not appear in [19, 20, 33], the analysis to be given below is much more delicate and complicate than those of [19, 20, 33], although they are similar conceptually. Since most proofs of this section are in the same lines as those of the proofs in Section 4 of [27], we shall omit some details if they are already given in [27], but shall provide them if there are meaningful differences.

We first cite the following lemma which establishes three integral identities and play a crucial role in our analysis. A proof of the lemma can be found in [27, Lemma 4.1].

Lemma 3.1**.**

Let $\alpha(x):=x-x_{\Omega_{1}}$, $v\in E^{1}$, $K,K^{\prime}\in\mathcal{T}_{h}$ and $e\in\mathcal{E}_{h}^{ID}$. Then there hold

[TABLE]

where $x_{\Omega_{1}}$ denotes the point in the star-shaped domain definition for $\Omega_{1}$ (see Definition 2.1).

Remark 3.1**.**

The identity (3.2) can be viewed as a local version of the Rellich identity for the Laplacian $\Delta$ (cf. [19, 20]). Since $V_{h}^{p}\subset E^{1}$, hence, (3.1)–(3.3) hold for any function $v=v_{h}\in V_{h}^{p}$.

We also need the following trace and inverse inequalities (cf. [41, 44, 12]).

Lemma 3.2**.**

For any $K\in\mathcal{T}_{h}$ and $z\in\mathcal{P}_{p}(K),$

[TABLE]

Now, taking $v_{h}=u_{h}^{q}$ in (2.15) yields

[TABLE]

Therefore, taking real part and imaginary part of the above equation and using (2.12) and (2.13) we get the following lemma.

Lemma 3.3**.**

Let $u_{h}^{q}\in V_{h}^{p}$ solve (2.15). Then

[TABLE]

From (3.5) and (3.6) we can bound \bigl{|}u_{h}^{q}\bigr{|}_{1,h} and the jumps in terms of $\left\|u_{h}^{q}\right\|_{L^{2}(\Omega)}^{2}$. In order to get the desired a priori estimates, we need to derive a reverse inequality whose coefficients can be controlled. Such a reverse inequality, which is often difficult to get under practical mesh constraints, and stability estimates for $u_{h}^{q}$ will be derived next.

Theorem 3.1**.**

Let $u_{h}^{q}\in V_{h}^{p}$ solve (2.15) and suppose $\beta_{1,e}\geq 0,\gamma_{0,e},\cdots,\gamma_{q,e}>0$. Then

[TABLE]

where

[TABLE]

Proof.

We divide the proof into three steps.

Step 1: Derivation of a representation identity for $\|u_{h}^{q}\|_{L^{2}(\Omega)}$. Define $v_{h}\in E$ by $v_{h}|_{K}=\alpha\cdot\nabla u_{h}^{q}|_{K}$ for every $K\in\mathcal{T}_{h}$. Since $v_{h}|_{K}$ is a polynomial of degree no more than $p$ on $K$, hence, $v_{h}\in V_{h}^{p}$. Using this $v_{h}$ as a test function in (2.15) and taking the real part of the resulted equation we get

[TABLE]

It follows from (3.1), (3.4), and (3.10) that (compare with (4.11) of [27])

[TABLE]

Using the identity \bigl{|}a\bigr{|}^{2}-\bigl{|}b\bigr{|}^{2}=\operatorname{{Re}}(a+b)(\bar{a}-\bar{b}) we have

[TABLE]

Using the identity again followed by the Rellich identity (3.2) we get (compare with (4.13) of [27])

[TABLE]

Plugging (3.12) and (3.13) into (3.11) gives (compare with (4.15) of [27])

[TABLE]

Step 2: Derivation of a reverse inequality. Our task now is to estimate each term on the right-hand side of (3.14). Since the terms on the first four lines can be bounded in the exactly same way as done in [27], we omit their derivations and only give the final results here for the reader’s convenience.

[TABLE]

The extra work here is to estimate the terms on the last two lines in (3.14). For an edge/face $e\in\mathcal{E}_{h}^{ID}$, let $\Omega_{e}$ denote the set of element(s) in $\mathcal{T}_{h}$ containing $e$ as one edge/face. By (3.3) we obtain

[TABLE]

The first term on the line six of (3.14) is bounded as follows (compare with (4.20) of [27]):

[TABLE]

The penalty term $L_{1}(\cdot,\cdot)$ is estimated as follows. Recall that $v_{h}|_{K}=\alpha\cdot\nabla u_{h}^{q}|_{K}$ with $\alpha=x-x_{\Omega_{1}}$ for each $K\in\mathcal{T}_{h}$. Noting that

[TABLE]

by the definition of $L_{1}(\cdot,\cdot)$ and Lemma 3.2 we get

[TABLE]

We remark that $\operatorname{{Im}}L_{1}(u_{h}^{q},v_{h})=0$ when $p=q=1$.

Next we estimate the penalty terms $J_{j}(u_{h}^{q},v_{h})$. Since $u_{h}^{q}$ and $v_{h}$ are piecewise polynomials of degree $p$ in general and those terms contain jumps of high order normal derivatives, it is quite delicate to control those terms as shown below.

By direct calculations we get that on each edge/face $e$ of $K\in\mathcal{T}_{h}$

[TABLE]

Here we have used the fact that $(p+1)$th order derivatives of $u_{h}^{q}$ is zero because $u_{h}^{q}$ is a polynomial of degree at most $p$.

For $j=1,2,\cdots q-1$, by (3.24) we have

[TABLE]

If $q<p$, then, from Lemma 3.2 and the inequality \left\|\dfrac{\partial^{q}\varphi}{\partial n_{e}^{q}}\right\|_{L^{2}(\partial K)}\lesssim ph_{e}^{-\frac{1}{2}}\bigl{|}\varphi\bigr{|}_{H^{q}(K)}, we have

[TABLE]

If $q=p$, (3.25) and the definition of $J_{p}(\cdot,\cdot)$ immediately imply that (compare with (4.14) of [27])

[TABLE]

The estimate for $\operatorname{{Im}}J_{0}(u_{h}^{q},v_{h})$ is similar to (3.26), so we get

[TABLE]

We also need the following estimate (compare with (4.22) of [27])

[TABLE]

where we have used the inverse inequality and the assumption that $D$ is star-shaped to derive the last inequality.

Step 3: Finishing up. We only prove the case of $q=p$ since the proof for $q<p$ is the same except using (3.27) instead of (3.28). Substituting (3.26), (3.28), (3.29) (with $q=p$) and (3.15)–(3.19) into (3.14), and using (3.30) we obtain

[TABLE]

Therefore, it follows from Lemma 3.3 and (3.9) that

[TABLE]

where we have used the following inequality, which is a consequence of (3.6),

[TABLE]

to derive the last inequality. $M(f,g)$ and $C_{{\rm sta},p}$ are defined by (3.8) and (3.9), respectively. Hence,

[TABLE]

which together with (3.6) gives (3.7). The proof is completed. ∎

As (2.15) can be written as a linear system, an immediate consequence of the above stability estimate is the following well-posedness theorem for (2.15).

Theorem 3.2**.**

For $k>0$ and $h>0$, the $hp$-IPDG method (2.15) has a unique solution $u_{h}^{q}$ provided that $\gamma_{0,e},\gamma_{1,e},\cdots\gamma_{q,e}>0$ and $\beta_{1,e}\geq 0$.

Next we consider the case of quasi-uniform meshes. Note that large penalty parameters ($\gamma_{j},j\geq 1$) for jumps of normal derivatives may cause a large interpolation error in the norm $\left\|\cdot\right\|_{1,h,q}$, and hence, may pollute the error estimates of the IPDG solution (see Section 4). It is interested to minimize the stability constant $C_{{\rm sta},q}$ under the constraints of $\beta_{1,e}\geq 0$ and $\frac{p}{\gamma_{0,e}}+\sum_{j=1}^{q}p^{2j-1}\gamma_{j,e}\lesssim 1.$ We have the following consequence of Theorem 3.1. The proof is straightforward and is omitted.

Theorem 3.3**.**

Let $h=\max h_{e}.$ Suppose the mesh $\mathcal{T}_{h}$ is quasi-uniform, that is $h_{e}\simeq h.$ Suppose $k\gtrsim 1$ and $k\,h\lesssim 1$. Assume that $\gamma_{j,e}\simeq\gamma_{j},j=0,1,\cdots,q$, $\sum_{j=1}^{q}p^{2j-1}\gamma_{j}\lesssim 1$, $0\leq\beta_{1,e}\lesssim\frac{h^{2}}{p^{4}}\gamma_{0}$, and that

[TABLE]

Suppose $D=\emptyset$ if $q\geq 2$ or $q=p=1$. Then

[TABLE]

where

[TABLE]

It is clear that, in the above theorem, $\gamma_{0}\simeq p^{\frac{7}{3}}h^{-\frac{2}{3}}$ if $q=p$ or $q\geq 2$, $\gamma_{0}\simeq p^{\frac{3q+1}{q+1}}h^{-\frac{q}{q+1}}$ otherwise. We conclude this section by several remarks.

Remark 3.2**.**

(a) The $hp$-IPDG method (2.15) is well-posed for all $h,k>0$ provided that all penalty parameters are positive. As a comparison, we recall that the standard finite element method is well-posed only if mesh size $h$ satisfies a constraint $h=O(k^{-\rho})$ for some $\rho\geq 1$, hence, the existence is only guaranteed for very small mesh size $h$ when wave number $k$ is large.

(b). It is well known that [3, 4, 9, 26, 40, 45] symmetric IPDG methods for coercive elliptic and parabolic PDEs often require the penalty parameter $\gamma_{0,e}$ is sufficiently large to guarantee the well-posedness of numerical solutions, and the low bound for $\gamma_{0,e}$ is theoretically hard to determine and is also problem-dependent. However, this is no issue for scheme (2.15), which solves the (indefinite) Helmholtz equation, because they are well-posed for all $\gamma_{0,e}>0$.

(c) In the linear element case of $q=p=1$, a better estimate is obtained in [27] with the help of the penalty term $L_{1}$. Unfortunately the penalty term $L_{1}$ does not help too much for higher order elements.

(d). The stability estimates will be improved greatly when $k^{3}h^{2}p^{-1}\leq C_{0}$, where $C_{0}$ is some constant independent of $k$, $h$, $p$, and the penalty parameters (see Theorem 5.1). The above estimates are only interested when $k\,h\lesssim 1$ and $k^{3}h^{2}p^{-1}\gtrsim 1$.

4. Error estimates

In this section, we derive the error estimates of the solutions of scheme (2.15). This will be done in two steps. First, we introduce elliptic projections of the PDE solution $u$ and derive error estimates for the projections. We note that such a result also has an independent interest. Second, we bound the error between the projections and the IPDG solutions by making use of the stability results obtained in Section 3. In this and the next section, we assume that the mesh $\mathcal{T}_{h}$ is quasi-uniform, that $\gamma_{j,e}\simeq\gamma_{j}>0$ for $j=0,1,\cdots,q$, and that $\beta_{1,e}\simeq\beta_{1}\geq 0$. Let $h=\max h_{e}$. For simplicity, we also assume that the mesh $\mathcal{T}_{h}$ is conforming, that is, $\mathcal{T}_{h}$ contains no hanging nodes, since the parallel results for non-conforming meshes can be derived in a similar way.

In this and the next section the following assumption (cf. Theorem 2.1) is required for some results. Suppose problem (1.1)–(1.3) is $H^{s}$-regular and

[TABLE]

where $s\geq 2$ is an integer and

[TABLE]

Note that the assumption is proved to hold for $s=2$ and $M_{2}(f,g)=M(f,g)$.

4.1. Elliptic projection and its error estimates

For any $w\in E^{q}\cap H^{1}_{\Gamma_{D}}(\Omega)\cap H^{2}_{\mbox{\tiny loc}}(\Omega)$, we define its elliptic projection $\tilde{w}_{h}^{q}\in V_{h}^{p}$ by

[TABLE]

In other words, $\tilde{w}_{h}$ is an IPDG approximation to the solution $w$ of the following (complex-valued) Poisson problem:

[TABLE]

for some given functions $F$ and $\psi$ which are determined by $w$.

Before estimating the projection error, we state the following continuity and coercivity properties for the sesquilinear form $a_{h}^{q}(\cdot,\cdot)$. Since they follow easily from (2.4)–(2.13), so we omit their proofs to save space.

Lemma 4.1**.**

For any $v\in E^{q}$ and $w\in E^{q}\cap H^{1}_{\Gamma_{D}}$, the mesh-dependent sesquilinear form $a_{h}^{q}(\cdot,\cdot)$ satisfies

[TABLE]

In addition, for any $0<\varepsilon<1$, there exists a positive constant $c_{\varepsilon}$ independent of $k$, $h$, $p$, and the penalty parameters such that

[TABLE]

Let $u$ be the solution of problem (1.1)–(1.3) and $\tilde{u}_{h}^{q}$ be its elliptic projection defined as above. Define $\eta:=u-\tilde{u}_{h}^{q}$. Then (4.2) immediately implies the following Galerkin orthogonality:

[TABLE]

Lemma 4.2**.**

Suppose problem (1.1)–(1.3) is $H^{\max\left\{q+1,2\right\}}$-regular. Then there hold the following estimates:

[TABLE]

where $\lambda:=1+\frac{p}{\gamma_{0}}$ and $\gamma_{1}=0$ if $q=0$.

Proof.

For any $z_{h}\in V_{h}^{p}\cap H_{\Gamma_{D}}^{1}(\Omega)$, let $\eta_{h}=\tilde{u}_{h}^{q}-z_{h}$. From $\eta_{h}+\eta=u-z_{h}$ and (4.5), we have

[TABLE]

Take $\varepsilon=\frac{1}{2}$ in (4.4) and assume without loss of generality that $c_{\frac{1}{2}}>\frac{1}{2}$. It follows from (4.4) and (4.8) that

[TABLE]

Therefore,

[TABLE]

which together with $\eta=u-z_{h}-\eta_{h}$ yields (4.6).

To show (4.7), we use the Nitsche’s duality argument (cf. [11, 15]). Consider the following auxiliary problem:

[TABLE]

It can be shown that $w$ satisfies

[TABLE]

Let $\hat{w}_{h}^{1}$ be the continuous linear finite element interpolant of $w$ on $\mathcal{T}_{h}$. From (4.5),

[TABLE]

Testing the conjugated (4.10) by $\eta$ and using the above orthogonality we get

[TABLE]

which together with (4.6) and (4.11) gives (4.7). The proof is completed. ∎

We have the following lemma that gives approximation properties of the space $V_{h}^{p}\cap H_{\Gamma_{D}}^{1}(\Omega)$.

Lemma 4.3**.**

**

1. (i)

Let $\mu=\min\left\{p+1,s\right\}$ and $q<\mu$. Suppose $u\in H^{s}(\Omega)\cap H_{\Gamma_{D}}^{1}(\Omega)$. Then there exists $\hat{u}_{h}\in V_{h}^{p}\cap H_{\Gamma_{D}}^{1}(\Omega)$ such that

[TABLE] 2. (ii)

Suppose $u\in H^{\max\left\{q+1,2\right\}}(\Omega)\cap H_{\Gamma_{D}}^{1}(\Omega)$. Then there exists $\hat{u}_{h}\in V_{h}^{p}\cap H_{\Gamma_{D}}^{1}(\Omega)$ such that (4.12) holds with $s=2$ and

[TABLE]

where $\gamma_{1}=0$ if $q=0.$

Proof.

The following $hp$ approximation properties are well-known for the $hp$ finite element functions (cf. [8, 30, 31]):

• •

There exists $\check{u}_{h}\in V_{h}^{p}$ such that, for $j=0,1,\cdots,s$,

[TABLE]

• •

There exists $\hat{u}_{h}\in V_{h}^{p}\cap H_{\Gamma_{D}}^{1}(\Omega)$ such that

[TABLE]

Here the invisible constants in the above two inequalities depending on $s$ but independent of $k$, $h$, $p$, and the penalty parameters. Then (4.12) follows from (4.16) and the trace inequality.

It follows from the inverse inequalities in Lemma 3.2 that, for $1\leq j\leq q+1$,

[TABLE]

Therefore, by the following local trace inequality

[TABLE]

we have

[TABLE]

Noting that $\hat{u}_{h}$ is continuous, we have from (4.16) and (4.18),

[TABLE]

That is, (4.13) holds.

(4.14) can be proved similarly as above. It is clear that (4.16) and (4.17) hold with $s=2$ and (4.18) holds with $s=2$ and $j=1$, that is,

[TABLE]

Since $u\in H^{q+1}(\Omega)$, we have from Lemma 3.2 and (4.21) that, for $2\leq j\leq q$,

[TABLE]

Now (4.14) follows by combining the equality in (4.19) and (4.20)–(4.23). This completes the proof of the lemma. ∎

By combining Lemma 4.2 and Lemma 4.3(i) we have the following estimates for the projection error.

Lemma 4.4**.**

Let $\mu=\min\left\{p+1,s\right\}$ and $q<\mu$. Suppose problem (1.1)–(1.3) is $H^{s}$-regular and (4.1) holds. Then there hold the following estimates:

[TABLE]

where

[TABLE]

Remark 4.1**.**

The requirement $q<s$ in the lemma is clear since the projection $\tilde{u}_{h}^{q}$ is not defined for $q>s$. However, for $q<s$, $\left\|u-\tilde{u}_{h}^{q}\right\|_{1,h,q}$ can be bounded without using full regularity of $u$, and such a bound is also useful (see Lemma 5.1).

4.2. Error estimates for $u-u_{h}^{q}$

In this subsection we shall derive error estimates for scheme (2.15). This will be done by exploiting the linearity of the Helmholtz equation and making use of the stability estimates derived in Theorem 3.1 and the projection error estimates established in Lemma 4.4.

Let $u$ and $u_{h}^{q}$ denote the solutions of (1.1)–(1.3) and (2.15), respectively. Assume that $u\in H^{s}(\Omega)$ with $s\geq q+1$, then (2.14) holds for $v_{h}\in V_{h}^{p}$. Define the error function $e_{h}:=u-u_{h}^{q}$. Subtracting (2.15) from (2.14) yields the following error equation:

[TABLE]

Let $\tilde{u}_{h}^{q}$ be the elliptic projection of $u$ as defined in the previous subsection. Write $e_{h}=\eta-\xi$ with $\eta:=u-\tilde{u}_{h}^{q},\xi:=u_{h}^{q}-\tilde{u}_{h}^{q}$. From (4.26) and (4.5) we get

[TABLE]

The above equation implies that $\xi\in V_{h}^{p}$ is the solution of scheme (2.15) with sources terms $f=-k^{2}\eta$ and $g\equiv 0$. Then an application of Theorem 3.1 and Lemma 4.4 immediately gives

Lemma 4.5**.**

$\xi=u_{h}^{q}-\tilde{u}_{h}^{q}$* satisfies the following estimate:*

[TABLE]

We are ready to state our error estimate results for scheme (2.15), which follows from Lemma 4.5, Lemma 4.4 and an application of the triangle inequality.

Theorem 4.1**.**

Let $u$ and $u_{h}^{q}$ denote the solutions of (1.1)–(1.3) and (2.15), respectively. Let $\mu=\min\{p+1,s\}$ and $q<\mu$. Assume assumption (4.1) holds. Then

[TABLE]

Remark 4.2**.**

$q<s$* is required in the theorem because $\left\|u-u_{h}^{q}\right\|_{1,h,q}$ is not defined for $q>s$. However, we note that the $hp$-IPDG solution $u_{h}^{q}$ is always well-defined regardless the regularity of underlying PDE solution $u$. For $q<s$, $\left\|u-u_{h}^{q}\right\|_{1,h,q}$ can be bounded without using full regularity of $u$, and such a bound is also useful (see Lemma 5.2).*

By combining Theorem 4.1 and Theorem 3.3 we have the following theorem that gives the best convergence order so far that we can obtain theoretically for the method (2.15) under the mesh condition $k^{3}h^{2}p^{-1}\gtrsim 1$ (cf. Theorem 5.1).

Theorem 4.2**.**

Under the assumptions of Theorem 3.3 and 4.1, we have

[TABLE]

Proof.

The proof is obvious since $C_{{\rm err},q}\simeq\widehat{C}_{{\rm err},q}\simeq 1.$ ∎

Remark 4.3**.**

(a) Estimates (4.29)–(4.31) are so-called preasymptotic error estimates which are suboptimal in $h$ and $k$. They can be improved to optimal order when $k^{3}h^{2}p^{-1}\leq C_{0}$, where $C_{0}$ is some constant independent of $k$, $h$, $p$, and the penalty parameters (see Theorem 5.1). The second term on the right-hand side of (4.29) is called a pollution term for $\left\|u-u_{h}^{q}\right\|_{1,h,q}$.

(b) Theorem 4.2 shows that $\left\|u-u_{h}^{q}\right\|_{1,h,q}+k\,\left\|u-u_{h}^{q}\right\|_{L^{2}(\Omega)}\to 0$ if $q=p$ and $p^{\frac{11}{3}-s}k^{s}h^{\mu-\frac{4}{3}}\to 0$, or if $q=1<p$ and $p^{4-s}k^{s}h^{\mu-\frac{3}{2}}\to 0$.

5. Stability-error iterative improvement

In this section we derive some improved optimal order stability and error estimates for the $hp$-IPDG solution under the mesh condition that $k^{3}h^{2}p^{-1}\leq C_{0}$ by using a stability-error iterative procedure, where $C_{0}$ is some constant independent of $k$, $h$, $p$, and the penalty parameters (see Theorem 5.1).

By combining Lemma 4.2 and Lemma 4.3(ii), we have the following estimates for the projection error when only the $H^{2}$-norm of the solution $u$ is allowed in the error bound.

Lemma 5.1**.**

Suppose problem (1.1)–(1.3) is $H^{\max\left\{q+1,2\right\}}$-regular. Then there hold the following estimates:

[TABLE]

where

[TABLE]

By a similar argument to that used to prove Theorem 4.1, we have the following error bounds which only involves $M(f,g)$.

Lemma 5.2**.**

Let $u$ and $u_{h}^{q}$ denote the solutions of (1.1)–(1.3) and (2.15), respectively. Suppose $u\in H^{\max\left\{q+1,2\right\}}(\Omega)\cap H_{\Gamma_{D}}^{1}(\Omega)$. Then

[TABLE]

We are now ready to state our final main theorem of this paper.

Theorem 5.1**.**

Let $u$ and $u_{h}^{q}$ denote the solutions of (1.1)–(1.3) and (2.15), respectively. Suppose $u\in H^{\max\left\{q+1,2\right\}}(\Omega)\cap H_{\Gamma_{D}}^{1}(\Omega)$. Assume that $k\gtrsim 1$, $k\,h\lesssim 1$, and that $p{\gamma_{0}}^{-1}+\sum_{j=1}^{q}p^{2j-1}\gamma_{j}\lesssim 1$. Then there exists a constant $C_{0}>0$, which is independent of $k$, $h$, $p$, and the penalty parameters, such that if $k^{3}h^{2}p^{-1}\leq C_{0}$, then the following stability estimates hold:

[TABLE]

Moreover, under the assumption (4.1), there hold the following error estimates:

[TABLE]

Proof.

We only prove (5.5) since (5.6) can be proved similarly and (5.7)–(5.8) follow from the improved stability estimates and the argument used in the proof of Theorem 4.1.

From Theorem 3.1 we have

[TABLE]

where $C_{{\rm sta},q}$ is defined in (3.9). From Lemma 5.2 we have

[TABLE]

Now it follows from Theorem 2.1 and the triangle inequality that

[TABLE]

Repeating the above process yields that there exists a constant $C_{1}$ independent of $k$, $h$, $p$, and the penalty parameters, and a sequence of positive numbers $\Lambda_{j}$ such that

[TABLE]

with

[TABLE]

A simple calculation yields that if $C_{1}\,\widehat{C}_{{\rm err},2,q}\,k^{3}h^{2}<\theta\,p$ for some positive constant $\theta<1$ then

[TABLE]

which implies (5.5) by noting that $C_{{\rm err},2,q},\widehat{C}_{{\rm err},2,q}\lesssim 1$ and that C_{{\rm err},2,q}\,k\,h\lesssim\big{(}C_{{\rm err},2,q}^{2}k^{3}h^{2}\big{)}^{\frac{1}{2}}\lesssim\big{(}\widehat{C}_{{\rm err},2,q}^{2}k^{3}h^{2}\big{)}^{\frac{1}{2}}\lesssim p^{\frac{1}{2}}. ∎

Note that the stability estimates in (5.5) and (5.6) are of the same order as the PDE stability estimates given in Theorem 2.1.