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

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.


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: where k ∈ R, called wave number, is a (large) positive number, D ⊂ Ω 1 ⊂ R d , d = 2, 3, D is known as a scatterer and is assumed to be a bounded Lipschitz domain, Ω 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 ∂Ω = Γ R ∪ Γ D .n Ω denotes the unit outward normal to Ω. i := √ −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 Γ R and to limit the computation of the original scattering problem just on the finite domain Ω.Boundary condition (1.3) implies that the scatterer is sound-soft.We note that the case D = ∅ 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 1norm error bound deteriorates as the wave number k increases under the practical "rule of thumb" mesh constraint kh 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 ∇u h • (x − x Ω ) (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 ≤ C 0 in Section 5, where C 0 is some constant independent of k, h, p, and the penalty parameters.

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 + ia i (a r and a i are real numbers), a := a r − ia i denotes the complex conjugate of a.
We also define 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 B and B A for the inequality A ≤ CB and B ≥ CA.A ≃ B is for the statement A B and B A.
We now give the definition of star-shaped domains.
In this paper, we assume that Ω 1 is a strictly star-shaped domain.Recall that Ω 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 Ω1 as Ω 1 does.This implies that x Ω1 ∈ D ⊂ Ω 1 .Under these assumptions, there hold following stability estimates for problem (1.1)-(1.3).
Theorem 2.1.Suppose Ω 1 ⊂ R d is a strictly star-shaped domain and D ⊂ Ω 1 is a star-shaped domain.Then the solution u to problem (1.1)- (1.3) 3) Ω) for some positive integer q ≥ 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 T h be a family of partitions of the domain Ω := Ω 1 \ D parameterized by h ∈ (0, h 0 ).For any K ∈ T h , we define h K := diam(K).Similarly, for each edge/face e of K ∈ T h , h e := diam(e).We impose the following mild restrictions on the partition T h : (i) The elements of T h satisfy the minimal angle condition, (ii) T h is locally quasi-uniform, that is if two elements K and K ′ are adjacent ( i.e., meas(∂K Where meas(e) stands for (d − 1)-dimensional Lebesgue measure of e.
For any two elements K, K ′ ∈ T h , we call e = ∂K ∩ ∂K ′ an interior edge/face of T h if meas(e) > 0. Note that e could be portion of a side/face of the element K or K ′ in the case of geometrically nonconforming partition.Also, for any element K ∈ T h , we call e = ∂K ∩ ∂Ω a boundary edge/face if meas(e) > 0. Then we define E I h := set of all interior edges/faces of T h , E R h := set of all boundary edges/faces of T h on Γ R , E D h := set of all boundary edges/faces of T h on Γ D , = set of all edges/faces of T h .
We also define the jump [v] of v on an interior edge/face e = ∂K ∩ ∂K ′ as The following convention is adopted in this paper If e ∈ E RD h , set {v}| e = v| e .For every e = ∂K ∩ ∂K ′ ∈ E I h , 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 ′ if the other way around.For every e ∈ E RD h , let n e = n Ω the unit outward normal to ∂Ω.
Let p ≥ 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 ≤ q ≤ p, we define the "energy" space and the sesquilinear form and σ is a real number.γ 0,e , • • • , γ q,e > 0 and β 1,e ≥ 0 are numbers to be specified later.{τ ℓ e } d−1 ℓ=1 denote an orthogonal coordinate frame on the edge/face e ∈ E h , ∂u ∂τ ℓ e := ∇u • τ ℓ e stands for the tangential derivative of u in the direction τ ℓ e , and ∂ j u ∂n j e denotes the jth order normal derivative of u on e.
On the other hand, when σ = 1, a q h (•, •) is non-symmetric.In particular, σ = −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 σ = 1.The penalty constants in , 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 : (2.10) It is easy to check that the sesquilinear form Using the sesquilinear form a q h (•, •) we now introduce the following weak formulation for The above formulation is consistent with the boundary value problem (1.1)-(1.2) because a q h (•, •) is consistent with −∆.For any K ∈ T h , let P p (K) denote the set of all polynomials whose degrees do not exceed p.We define our hp-IPDG approximation space V p h as V p h := We are now ready to define our hp-IPDG methods based on the weak formulation (2.14):For each 0 ≤ q ≤ p, find , 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.
(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 , • • • , 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.

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 q h 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 wellposedness 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 = α•∇u q h (defined element-wise) with α(x) := x − x Ω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 p h (•, •), 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].Remark 3.1.The identity (3.2) can be viewed as a local version of the Rellich identity for the Laplacian ∆ (cf.[19,20]).
Lemma 3.2.For any K ∈ T h and z ∈ P p (K), Therefore, taking real part and imaginary part of the above equation and using (2.12) and (2.13) we get the following lemma.
From (3.5) and (3.6) we can bound u q h 1,h and the jumps in terms of u q h 2 L 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 q h will be derived next.
Step 1: Derivation of a representation identity for Using this v h as a test function in (2.15) and taking the real part of the resulted equation we get It follows from (3.1), (3.4), and (3.10) that (compare with (4.11) of [27]) Using the identity a Using the identity again followed by the Rellich identity (3.2) we get (compare with (4.13) of [27]) Plugging (3.12) and (3.13) into (3.11)gives (compare with (4.15) of [27]) 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.

2(d − 1)
The extra work here is to estimate the terms on the last two lines in (3.14).For an edge/face e ∈ E ID h , let Ω e denote the set of element(s) in T h containing e as one edge/face.By (3.3) we obtain 2 The first term on the line six of (3.14) is bounded as follows (compare with (4.20) of [27]): The penalty term .
We remark that Im L 1 (u q h , v h ) = 0 when p = q = 1.Next we estimate the penalty terms J j (u q h , v h ).Since u q h 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 ∈ T h

.25)
Here we have used the fact that (p + 1)th order derivatives of u q h is zero because u q h is a polynomial of degree at most p. .
If q < p, then, from Lemma 3.2 and the inequality ∂ q ϕ ∂n q e L 2 (∂K) we have .
If q = p, (3.25) and the definition of J p (•, •) immediately imply that (compare with (4.14) of [27]) The estimate for Im J 0 (u q h , v h ) is similar to (3.26), so we get .
We also need the following estimate (compare with (4.22) of [27]) where we have used the inverse inequality and the assumption that D is star-shaped to derive the last inequality.
Therefore, it follows from Lemma 3.3 and (3.9) that where we have used the following inequality, which is a consequence of (3.6), to derive the last inequality.M (f, g) and C sta,p are defined by (3.8) and (3.9), respectively.Hence, C sta,p M (f, g), 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 q h provided that γ 0,e , γ 1,e , • • • γ q,e > 0 and β 1,e ≥ 0. Next we consider the case of quasi-uniform meshes.Note that large penalty parameters (γ j , j ≥ 1) for jumps of normal derivatives may cause a large interpolation error in the norm • 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 sta,q under the constraints of β 1,e ≥ 0 and p γ0,e + q j=1 p 2j−1 γ j,e 1.We have the following consequence of Theorem 3.1.The proof is straightforward and is omitted.
p 4 γ 0 , and that γ 0 ≃ p where It is clear that, in the above theorem, γ 0 ≃ p 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 −ρ ) for some ρ ≥ 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 γ 0,e is sufficiently large to guarantee the well-posedness of numerical solutions, and the low bound for γ 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 γ 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 ≤ 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 1 and k 3 h 2 p −1 1.

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 T h is quasiuniform, that γ j,e ≃ γ j > 0 for j = 0, 1, • • • , q, and that β 1,e ≃ β 1 ≥ 0. Let h = max h e .For simplicity, we also assume that the mesh T h is conforming, that is, 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 where s ≥ 2 is an integer and 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 ∈ E q ∩ H 1 ΓD (Ω) ∩ H 2 loc (Ω), we define its elliptic projection wq h ∈ V p h by (4.2) In other words, wh is an IPDG approximation to the solution w of the following (complex-valued) Poisson problem: for some given functions F and ψ which are determined by w.
Before estimating the projection error, we state the following continuity and coercivity properties for the sesquilinear form a q h (•, •).Since they follow easily from (2.4)-(2.13),so we omit their proofs to save space.Lemma 4.1.For any v ∈ E q and w ∈ E q ∩ H 1 ΓD , the mesh-dependent sesquilinear form a q h (•, •) satisfies In addition, for any 0 < ε < 1, there exists a positive constant c ε independent of k, h, p, and the penalty parameters such that Let u be the solution of problem (1.1)-(1.3)and ũq h be its elliptic projection defined as above.Define η := u − ũq h .Then (4.2) immediately implies the following Galerkin orthogonality: 3) is H max{q+1,2} -regular.Then there hold the following estimates: where λ := 1 + p γ0 and γ 1 = 0 if q = 0. a q h (η h , η h ) + ik η h , η h ΓR = a q h (u − z h , η h ) + ik u − z h , η h ΓR .Take ε = 1 2 in (4.4) and assume without loss of generality that c 1 2 > 1 2 .It follows from (4.4) and (4.8) that

Let ŵ1
h be the continuous linear finite element interpolant of w on T h .From (4.5), Testing the conjugated (4.10) by η and using the above orthogonality we get 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 (i) Let µ = min {p + 1, s} and q < µ.
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 ≤ j ≤ q + 1, Therefore, by the following local trace inequality we have Noting that ûh is continuous, we have from (4.16) and (4.18), 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, where Remark 4.1.The requirement q < s in the lemma is clear since the projection ũq h is not defined for q > s.However, for q < s, u − ũq h 1,h,q can be bounded without using full regularity of u, and such a bound is also useful (see Lemma 5.1).

Error estimates for u−u q
h .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 q h denote the solutions of (1.1)-(1.3)and (2.15), respectively.Assume that u ∈ H s (Ω) with s ≥ q + 1, then (2.14) holds for v h ∈ V p h .Define the error function e h := u − u q h .Subtracting (2.15) from (2.14) yields the following error equation:

Let ũq
h be the elliptic projection of u as defined in the previous subsection.Write e h = η − ξ with η := u − ũq h , ξ := u q h − ũq h .From (4.26) and (4.5) we get The above equation implies that ξ ∈ V p h is the solution of scheme (2.15) with sources terms f = −k 2 η and g ≡ 0. Then an application of Theorem 3.1 and Lemma 4.4 immediately gives Lemma 4.5.ξ = u q h − ũq h satisfies the following estimate: 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 q h denote the solutions of (1.1)-(1.3)and (2.15), respectively.Let µ = min{p + 1, s} and q < µ.Assume assumption (4.1) holds.Then Remark 4.2.q < s is required in the theorem because u − u q h 1,h,q is not defined for q > s.However, we note that the hp-IPDG solution u q h is always well-defined regardless the regularity of underlying PDE solution u.For q < s, u − u q h 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 1 (cf.Theorem 5.1).Theorem 4.2.Under the assumptions of Theorem 3.3 and 4.1, we have u − u q h 1,h,q + k u − u q h L 2 (Ω) (4.31) Proof.The proof is obvious since C err,q ≃ C err,q ≃ 1.They can be improved to optimal order when k 3 h 2 p −1 ≤ 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 u − u q h 1,h,q .(b) Theorem 4.2 shows that u − u q h 1,h,q + k u − u q h L 2 (Ω) → 0 if q = p and p 11 3 −s k s h µ− 4 3 → 0, or if q = 1 < p and p 4−s k s h µ− 3 2 → 0.

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 ≤ 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.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 q h denote the solutions of (1.1)-(1.3)and (2.15), respectively.Suppose u ∈ H max{q+1,2} (Ω) ∩ H 1 ΓD (Ω).Then u − u q h 1,h,q C err,2,q + C sta,q k 3 h C err,2,q M (f, g) k h p , (5.3) u − u q h L 2 (Ω) p C err,2,q 1 k + C sta,q k M (f, g) k 2 h 2 p 2 .(5.4) We are now ready to state our final main theorem of this paper.
From Theorem 3.1 we have (5.9)u q h 1,h,q k C sta,q M (f, g), where C sta,q is defined in (3.9).From Lemma 5.2 we have (5.10)u − u q h 1,h,q C err,2,q + k C sta,q k 2 h C err,2,q k h p M (f, g).
Now it follows from Theorem 2.1 and the triangle inequality that u q h 1,h,q ≤ u 1,h,q + u − u q h 1,h,q = |u| 1,h + u − u q h 1,h,q (5.11) 1 + C err,2,q k h p + k 3 h 2 p C err,2,q k C sta,q M (f, g).
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 Λ j such that u q h 1,h,q ≤ Λ j M (f, g), (5.12) with Λ 0 ≃ k C sta,q , Λ j = C 1 (1 + C err,2,q k h p ) + C 1 C err,2,q k 3 h 2 p Λ j−1 , j = 1, 2, • • • .

,
where x Ω1 denotes the point in the star-shaped domain definition for Ω 1 (see Definition 2.1).