A Hybridizable Discontinuous Galerkin Method for the Helmholtz Equation with High Wave Number

# A Hybridizable Discontinuous Galerkin Method for the Helmholtz Equation with High Wave Number

Huangxin Chen  , Peipei Lu, and Xuejun Xu School of Mathematical Sciences, Xiamen University, Xiamen, 361005, People’s Republic of China. chx@xmu.edu.cn.LSEC, Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, P.O.Box 2719, Beijing, 100190, People’s Republic of China. lupeipei@lsec.cc.ac.cn, xxj@lsec.cc.ac.cn.
###### Abstract

This paper analyzes the error estimates of the hybridizable discontinuous Galerkin (HDG) method for the Helmholtz equation with high wave number in two and three dimensions. The approximation piecewise polynomial spaces we deal with are of order . Through choosing a specific parameter and using the duality argument, it is proved that the HDG method is stable without any mesh constraint for any wave number . By exploiting the stability estimates, the dependence of convergence of the HDG method on and is obtained. Numerical experiments are given to verify the theoretical results.

Key words. Hybridizable discontinuous Galerkin method, Helmholtz equation, high wave number, error estimates

## 1 Introduction

The numerical solutions of Helmholtz problems have been an area of active research for almost half of a century. Because of the well known pollution effect, the standard Galerkin finite element methods can maintain a desired level of accuracy only if the mesh resolution is also appropriately increased. In order to remedy this problem and to obtain more stable and accurate approximation, numerous nonstandard methods have been proposed recently (cf. [21]). One type of methods applies the stabilized discrete variational form to approximate the Helmholtz equation, which includes Galerkin-least-squares finite element methods [10, 19, 25], quasi-stabilized finite element methods [5], absolutely stable discontinuous Galerkin (DG) methods [14, 15, 16] and continuous interior penalty finite element methods (CIP-FEM) [29]. Other approaches include the partition of unity finite element methods [3, 24, 26], the ultra weak variational formulation [9], plane wave DG methods [2, 20], spectral methods [27], generalized Galerkin/finite element methods [7, 23], meshless methods [6], and the geometrical optics approach [13].

Discontinuous Galerkin methods have several attractive features compared with conforming finite element methods. For example, the polynomial degrees can be different from element to element, and they work well on arbitrary meshes. For the Helmholtz equation, the interior penalty discontinuous Galerkin methods (cf. [14, 15]) and the local discontinuous Galerkin methods [16] perform much better than the standard finite element methods, and they are well posed without any mesh constraint. Despite all these advantages, the dimension of the approximation DG space is much larger than the dimension of the corresponding classical conforming space.

The hybridizable discontinuous Galerkin methods were recently introduced to try to address this issue. The HDG methods retain the advantages of the standard DG methods and result in a significant reduced degrees of freedom. New variables on the boundary of elements are introduced such that the solution inside each element can be computed in terms of them. In particular, element by element, volume degrees of freedom can be parameterized by the surface degrees, and the resulting algebraic system is only due to the unknowns on the skeleton of the mesh. For a comprehensive understanding of the HDG methods, we can refer to [11] for a unified framework for second order elliptic problems and to [22] for the implementations.

In [17], the authors give error estimates of the HDG method for the interior Dirichlet problem for the Helmholtz equation, but it is under the condition that and are sufficiently small, where is a constant which is dependent on but is not characterized explicitly, and depend on the parameters defined in the numerical fluxes. Motivated by this work, the primary objective of this paper is to analyze the explicit dependence of convergence of HDG method for the Helmholtz equation on and . In this paper, we consider the Helmholtz equation with Robin boundary condition which is the first order approximation of the radiation condition:

 −△u−κ2u =~fin Ω, (1.1) ∂u∂n+iκu =~gon ∂Ω, (1.2)

where , is a polygonal/polyhedral domain, is known as the wave number, denotes the imaginary unit, and denotes the unit outward normal to .

The main difficulty of analyzing the Helmholtz equation lies in the strong indefiniteness of the problem which makes it hard to establish the stability for the numerical approximation. For the HDG method, we use a duality argument to obtain the stability estimates of the numerical solution. In the analysis, a crucial step lies in the derivation of the dependence of convergence on . We utilize the explicit error estimates of projection operator (see Lemma 4.4) to overcome this problem. Then we obtain that the HDG method for the Helmholtz problem (1.1)-(1.2) attains a unique solution for any , . Furthermore, the stability results not only guarantee the well-posedness of the HDG method but also play a key role in the derivation of the error estimates.

The duality argument can not be directly applied to establish the error estimates. Thus, we first construct an auxiliary problem and show its HDG error estimates by the duality technique. Then, combining the stability estimates, the error estimates of HDG scheme for the original Helmholtz problem (1.1)-(1.2) are deduced. Let and be the HDG approximation to and respectively. We obtain the following results:

(i) The following stability and error estimates hold without any constraint:

 ∥uh∥0,Ω+∥qh∥0,Ω≲(1+κ3h2p2)∥f∥0,Ω+(1+κ32hp)∥g∥0,∂Ω,
 ∥u−uh∥0,Ω≲(κh2p2+κ2h2p2+κ5h4p4)M(~f,~g),
 κ∥q−qh∥0,Ω≲(κhp+κ3h2p2+κ6h4p4)M(~f,~g),

where , and . We use notations and for the inequalities and , where is a positive number independent of the mesh size, polynomial degree and wave number , but the value of which can take on different values in different occurrences.

(ii) Suppose , there hold the following improved results:

 ∥uh∥0,Ω+∥qh∥0,Ω≲∥f∥0,Ω+∥g∥0,∂Ω,
 ∥u−uh∥0,Ω≲(κh2p2+κ2h2p2)M(~f,~g),
 κ∥q−qh∥0,Ω≲(κhp+κ3h2p2)M(~f,~g).

Comparing to the estimates in -IPDG method for Helmholtz problem, we find that the condition for the above improved results weakens the mesh condition which is requested in [15]. For the estimates under the mesh condition , the results in (i) can not be directly applied, but we may still get the following improved estimates.

(iii) Suppose , there hold

 ∥u−uh∥0,Ω≲κ2h2p2M(~f,~g),
 κ∥q−qh∥0,Ω≲(κhp+κ3h2p2)M(~f,~g).

We remark that in this work the local stabilization parameter to determine the numerical flux in the HDG scheme is always selected as (see (6.1)). Our numerical results show that the predicted convergence rates are observed.

The organization of the paper is as follows: We precisely define the HDG method for the Helmholtz equation and give some notations in the next section. Section 3 is dedicated to the characterization of the surface degrees . In section 4, we derive the stability estimates of the HDG method. The error estimates of the auxiliary problem are carried out in section 5 while section 6 states the main results of this paper, i.e., the error estimates of the HDG method for the Helmholtz equation. In the final section, we give some numerical results to confirm our theoretical analysis.

## 2 The hybridizable discontinuous Galerkin method

The HDG scheme is based on a first order formulation of the above Helmholtz equation (1.1)-(1.2) which can be rewritten in mixed form as finding such that

 iκq+∇u =0in Ω, (2.1) iκu+divq =fin Ω, (2.2) −q⋅n+u =gon ∂Ω. (2.3)

Existence and uniqueness of solutions to (2.1)-(2.3) is well known and it is proved in [16] that they satisfy the following regularity result:

 κ−1∥u∥2,Ω+∥u∥1,Ω+∥q∥1,Ω≲M(~f,~g). (2.4)

We consider a subdivision of into a finite element mesh of shape regular triangle in (or tetrahedron in ) and denote the collection of triangles (tetrahedra) by , the collection of edges (faces) by , while the collection of interior edges (faces) by and the collection of element boundaries by . Throughout this paper we use the standard notations and definitions for Sobolev spaces (see, e.g.,Adams[1]).

On each element and each edge (face) , we define the local spaces of polynomials of degree :

 V(T):=(Pp(T))d,W(T):=Pp(T),M(F):=Pp(F),

where denotes the space of polynomials of total degree at most on . The corresponding global finite element spaces are given by

 Vph: ={v∈L2(Ω) | v|T∈V(T) for all T∈Th}, Wph: ={w∈L2(Ω) | w|T∈W(T) for all T∈Th}, Mph: ={μ∈L2(Eh) | μ|F∈M(F) for all F∈Eh},

where and . On these spaces we define the bilinear forms

 (v,w)Th:=∑T∈Th(v,w)T, (v,w)Th:=∑T∈Th(v,w)T, and ⟨v,w⟩∂Th:=∑T∈Th⟨v,w⟩∂T,

with , and .

The hybridizable discontinuous Galerkin method yields finite element approximations which satisfy

 (iκqh,¯¯¯r)Th−(uh,¯¯¯¯¯¯¯¯¯¯¯divr)Th+⟨^uh,¯¯¯¯¯¯¯¯¯¯¯r⋅n⟩∂Th =0, (2.5) (iκuh,¯¯¯¯w)Th−(qh,¯¯¯¯¯¯¯¯¯∇w)Th+⟨^qh⋅n,¯¯¯¯w⟩∂Th =(f,¯¯¯¯w)Th, (2.6) ⟨−^qh⋅n+^uh,¯¯¯μ⟩∂Ω =⟨g,¯¯¯μ⟩∂Ω, (2.7) ⟨^qh⋅n,¯¯¯μ⟩∂Th∖∂Ω =0, (2.8)

for all , , and , where the overbar denotes complex conjugation. The numerical flux is given by

 ^qh=qh+τ(uh−^uh)non ∂Th, (2.9)

where the parameter is the so-called local stabilization parameter which has an important effect on both the stability of the solution and the accuracy of the HDG scheme. We always choose in this paper. The error analysis is based on projection operators which are defined as follows

 Πh:L2(Ω)→VphΠh:L2(Ω)→Wph

for any , they satisfy

 (Πhq,v)T =(q,v)Tfor all v∈V(T), (2.10) (Πhu,w)T =(u,w)Tfor all w∈W(T). (2.11)

We conclude the introduction by setting some notations used throughout this paper. Let the broken space be defined by

 H1(Ωh):={v:v|T∈H1(T), ∀T∈Th},

the seminorm of which is

 |v|21,Ωh:=∑T∈Th|v|21,T.

The trace of functions in belong to . For any , and , if , we set

 {{ϕ}}:=12(ϕ++ϕ−),\llbracketϕ\rrbracket:=ϕ+⋅n++ϕ−⋅n−

and

 {{v}}:=12(v++v−),\llbracketv\rrbracket:=v+⋅n++v−⋅n−.

For , we define

 {{ϕ}}:=ϕ, \llbracketϕ\rrbracket:=ϕ⋅n, {{v}}:=v, \llbracketv\rrbracket:=v⋅n.

## 3 The characterization of ^uh

One of the advantages of hybridizable discontinuous Galerkin methods is the elimination of both and from the equation and obtain a formulation in terms of only. In this section, we show that can be characterized by a simple weak formulation in which none of the other variables appear.

First we define the discrete solutions of the local problems, for each function , satisfies the following formulation

 (iκQλ,¯¯¯r)T−(Uλ,¯¯¯¯¯¯¯¯¯¯¯divr)T =−⟨λ,¯¯¯¯¯¯¯¯¯¯¯r⋅n⟩∂Tfor all r∈V(T), (3.1) (iκUλ,¯¯¯¯w)T−(Qλ,¯¯¯¯¯¯¯¯¯∇w)T+⟨^Qλ⋅n,¯¯¯¯w⟩∂T =0for all w∈W(T), (3.2)

where . For , is defined as follows

 (iκQf,¯¯¯r)T−(Uf,¯¯¯¯¯¯¯¯¯¯¯divr)T =0for all r∈V(T), (3.3) (iκUf,¯¯¯¯w)T−(Qf,¯¯¯¯¯¯¯¯¯∇w)T+⟨^Qf⋅n,¯¯¯¯w⟩∂T =(f,¯¯¯¯w)Tfor all w∈W(T), (3.4)

where . Next we show that the local problem (3.1)-(3.2) is well posed. The uniqueness of (3.3)-(3.4) can be deduced similarly.

###### Lemma 3.1.

There exist a unique solution to the local problem (3.1)-(3.2).

###### Proof.

Since it is a square system, to prove the existence and uniqueness of its solution, it is enough to show that if , we have that , . Taking in (3.1), in (3.2), we get

 iκ(Uλ,¯¯¯¯¯¯Uλ)T−iκ(Qλ,¯¯¯¯¯¯¯Qλ)T+τ⟨Uλ,¯¯¯¯¯¯Uλ⟩∂T=0,

which means . Back to (3.1), we derive that

 −iκQλ=∇Uλ.

Inserting the above expression into (3.2),

 △Uλ=−κ2Uλ

is deduced, which implies that , . ∎

It is worth noting that the solution in (2.5)-(2.8) is exactly correspond to the following relationship

 qh=Q^uh+Qf,uh=U^uh+Uf.

And is the solution of the following formulation

 ah(^uh,μ)=bh(μ)for all μ∈Mph,

where

 ah(λ,μ): =−⟨\llbracket^Qλ\rrbracket,¯¯¯μ⟩∂Th+⟨λ,¯¯¯μ⟩∂Ω, bh(μ): =⟨\llbracket^Qf\rrbracket,¯¯¯μ⟩∂Th+⟨g,¯¯¯μ⟩∂Ω.

## 4 The stability of the hybridizable discontinuous Galerkin method

The goal of this section is to derive stability estimates. We first cite the following lemma which provides some approximation results that will play an important role later. A proof of the lemma can be found in [4, 28].

###### Lemma 4.1.

Let be a standard square or triangle. Then there exists an operator such that for any

 ∥^u−^πp^u∥0,^T ≲p−1|^u|1,^T. (4.1)

Moreover, if ,

 ∥^u−^πp^u∥0,^T ≲p−2|^u|2,^T, (4.2) |^u−^πp^u|1,^T ≲p−1|^u|2,^T. (4.3)

Using the standard scaling technique, we can get the following approximation results.

###### Lemma 4.2.

For any , there exists an operator such that for any there holds

 ∥u−πphu∥0,T ≲hp|u|1,T, (4.4)

Moreover, if ,

 ∥u−πphu∥0,T ≲(hp)2|u|2,T, (4.5) |u−πphu|1,T ≲hp|u|2,T. (4.6)

We also need the following trace inequality and refer to [28] for the proof.

###### Lemma 4.3.

For any and

 ∥v∥0,∂T≲ph−12∥v∥0,T.

Now we derive the following approximation properties of the projection operator which is defined in (2.11). For the sake of simplicity, the proof is restricted to 2-d case.

###### Lemma 4.4.

For any , , the projection operator satisfies

 ∥u−Πhu∥0,T ≲hp|u|1,T, (4.7) ∥u−Πhu∥0,∂T ≲(hp)12|u|1,T. (4.8)

Moreover, if

 ∥u−Πhu∥0,T ≲(hp)2|u|2,T, (4.9) ∥u−Πhu∥0,∂T ≲(hp)32|u|2,T. (4.10)
###### Proof.

An important property of projection operator is that

 ∥u−Πhu∥0,T≤infv∈W(T)∥u−v∥0,T,

hence (4.4) and (4.5) imply (4.7) and (4.9). Let and be the standard triangles with linear mappings and respectively, see Figure 1 for illustration. For any , define and .

For any , we have

 |^v(^x)|2 =^v(^x)⋅¯¯¯¯¯¯¯¯¯¯^v(^x)=^v(^x1,0)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯^v(^x1,0)+∫^x20∂(^v(^x1,η)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯^v(^x1,η))∂η =|^v(^x1,0)|2+2Re∫^x20^v(^x1,η)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,η)∂ηdη.

Hence

 |^v(^x1,0)|2=|^v(^x1,^x2)|2−2Re∫^x20^v(^x1,η)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,η)∂ηdη.

By using integration on we can deduce

 ∫10|^v(^x1,0)|2(1−^x1)d^x1=∫10∫1−^x10|^v(^x1,0)|2d^x2d^x1 =∫^T1|^v(^x1,^x2)|2d^x2d^x1−2Re∫10∫1−^x10∫^x20^v(^x1,η)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,η)∂ηdηd^x2d^x1 =∫^T1|^v(^x1,^x2)|2d^x2d^x1−2Re∫10∫1−^x10∫1−^x1η^v(^x1,η)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,η)∂ηd^x2dηd^x1 =∫^T1|^v(^x1,^x2)|2d^x2d^x1−2Re∫10∫1−^x10(1−^x1−η)^v(^x1,η)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,η)∂ηdηd^x1 =∫^T1|^v(^x1,^x2)|2d^x2d^x1−2Re∫^T1(1−^x1−^x2)^v(^x1,^x2)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂^v(^x1,^x2)∂^x2d^x2d^x1.

Take , where satisfies , and note that

 (u−Πhu,v)T=|T||^T1|(^u−ˆΠhu,^v)^T1=0,

which means . By Lemma 4.2 and scaling technique we get

 ∫^e((^u−^Π^u)(^x1,0))2(1−^x1)d^x1=∫^T1(^u−^Π^u)2d^x2d^x1−2Re∫^T1(^u−^Π^u)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂(^u−^Π^u)∂^x2(1−^x1−^x2)d^x2d^x1=∫^T1(^u−^Π^u)2d^x2d^x1−2Re∫^T1(^u−^Π^u)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂(^u−^πp^u)∂^x2(1−^x1−^x2)d^x2d^x1≤∥^u−^Π^u∥20,^T1+2∥^u−^Π^u∥0,^T1|^u−^πp^u|1,^T1≲h−2∥u−Πhu∥20,T+h−1∥u−Πhu∥0,T|u−πphu|1,T≲h2p3|u|22,T, (4.11)

and

 ∫^e((^u−^Π^u)(^x1,0))2(1−^x1)d^x1=∫^T1(^u−^Π^u)2d^x2d^x1−2Re∫^T1(^u−^Π^u)¯¯¯¯¯¯¯¯¯¯¯∂^u∂^x2(1−^x1−^x2)d^x2d^x1≲∥^u−^Π^u∥20,^T1+∥^u−^Π^u∥0,^T1|^u|1,^T1≲h−2∥u−Πhu∥20,T+h−1∥u−Πhu∥0,T|u|1,T≲p−1|u|21,T. (4.12)

Now we map to and similarly we can derive

 ∫^e((~u−~Π~u)(^x1,0))2^x1d^x1≲h2p3|u|22,T, (4.13)

and

 ∫^e((~u−~Π~u)(^x1,0))2^x1d^x1≲p−1|u|21,T, (4.14)

where satisfies and . Since on , summing up (4.11), (4.13), and (4.12), (4.14) respectively and noting that is not particularly chosen, the lemma is proved. ∎

###### Remark 4.1.

In this paper, we only deal with meshes consisting of triangles or tetrahedra, but we should note that Lemma 4.4 can be extended to the meshes constituted with rectangles or hexahedra. The proof is similar with Lemma 4.4 and can also be found in [8].

###### Lemma 4.5.

Let be the solutions of (2.5)-(2.8). There hold

 τ∥uh−^uh∥20,∂Th ≤∥f∥0,Ω∥uh∥0,Ω+∥g∥20,∂Ω, (4.15) κ∥qh∥