Error analysis for a hybridizable discontinuous Galerkin method for the ...
Error analysis for a hybridizable discontinuous Galerkin method for the Helmholtz equation Roland Griesmaier∗
Peter Monk†
Department of Mathematical Sciences, University of Delaware, Newark, Delaware 19716, U.S.A.
Abstract Finite element methods for acoustic wave propagation problems at higher frequency result in very large matrices due to the need to resolve the wave. This problem is made worse by discontinuous Galerkin methods that typically have more degrees of freedom than similar conforming methods. However hybridizable discontinuous Galerkin methods offer an attractive alternative because degrees of freedom in each triangle can be cheaply removed from the global computation and the method reduces to solving only for degrees of freedom on the skeleton of the mesh. In this paper we derive new error estimates for a hybridizable discontinuous Galerkin scheme applied to the Helmholtz equation. We also provide extensive numerical results that probe the optimality of these results. An interesting observation is that, after eliminating the internal element degrees of freedom, the condition number of the condensed hybridized system is seen to be almost independent of the wave number.
Keywords
discontinuous Galerkin methods, hybridization, error analysis, Helmholtz equation.
AMS subject classification: Short title:
1
65N15, 65N30, 35J05.
HDG for the Helmholtz equation
Introduction
Discontinuous Galerkin (DG) methods have several attractive features when applied to the Helmholtz equation. In particular the polynomial degree can be changed from element to element to compensate for different element sizes in order to control dispersion, and some types of non-standard meshes can be handled. However, they possess an important disadvantage: the dimension of the standard DG piecewise polynomial space is much larger than the dimension of the corresponding conforming space. This increase in dimension is particularly troublesome for wave propagation problems because as the wave-number κ increases (i.e., as the wave length decreases), the dimension of the space increases faster than O(κn ) in n = 2, 3 dimensions and hence becomes very large for high-frequency problems. In addition there is currently, to our knowledge, no optimal solver for the discrete linear system resulting from a finite element discretization of the Helmholtz equation. Since many practical problems involve high-frequency sound (for example simulations of ultra-sound propagation), DG methods appear less practical than conforming elements. Recently, however, hybridizable DG (HDG) methods have been developed that are the focus of this paper (see Cockburn, Gopalakrishnan, and Lazarov [4], Cockburn, Dong and Guzman [3]). HDG methods are DG ∗ e-mail: † e-mail:
[email protected] [email protected] (corresponding author).
1
methods in which it is possible to introduce new variables on the boundary of elements such that the solution inside each element can be computed in terms of these new variables. In particular, element by element, volume degrees of freedom can be cheaply parameterized by the surface degrees. The linear system can then be reduced, at least logically, to solving only for the unknowns on the skeleton of the mesh (i.e. the edges or faces of the mesh depending on the dimension of the space for the Helmholtz equation). Thus the size of the linear system is considerably reduced compared to a standard DG scheme, and objections based on the bloated dimension of DG methods are no longer relevant. This paper is motivated by the analysis of HDG methods via projections due to Cockburn, Gopalakrishnan, and Sayas in [6]. That analysis, applied to the HDG solution of the Poisson problem, shows that the scalar variable in the discrete problem possesses super-convergence to an appropriate projection of the exact solution in the L2 norm. This is attractive from the point of view of the Helmholtz equation. The standard approach to error analysis of the H 1 -conforming finite element approximation of the Helmholtz equation of Schatz [10] uses duality and L2 estimates for the scalar variable. This is then combined with a Gårding inequality to show that, as the mesh size h approaches zero, the conforming finite element approximation approaches the elliptic projection of the solution. The analysis in [6] suggests that a modification of Schatz’ analysis can be used, and we shall indeed show this is the case. Besides proving error estimates, we also provide some numerical results when the method is applied to the interior Dirichlet problem for the Helmholtz equation (away from a Dirichlet eigenvalue, also called a resonance of the domain). Our numerical results show that the predicted asymptotic convergence rates are observed at practical mesh sizes, but suggest that the regularity required of the solution by our theory is too high. We also test the method for non-smooth solutions and on a non-convex domain. The numerical results show that it can be used even when our analysis is not applicable, with the expected rates of convergence. This suggests that optimal convergence can be proved even for non-convex domain as we might expect. An important observation is that, after eliminating the internal element degrees of freedom, the condensed hybridized system for the boundary degrees of freedom has a condition number that is roughly constant as the wave number increases provided we stay away from resonances (until the mesh becomes too coarse to resolve the wave). It may be that hybridized methods offer an approach that will allow the development of wave number robust linear system solvers. The major contributions of our paper are the first h-version analysis of an HDG scheme for the Helmholtz equation which involves also complex valued functions appropriate to this problem. In addition our numerical results point to specific issues related to the optimality of this analysis. Finally our observation concerning the wave number dependence of the condition number of the condensed hybridized system implies that the condensed hybridized system displays better properties with respect to wave number than the global system. Our analysis, like that in [6], is devoted to the h-version of the method. The most desirable extension to our analysis would include the dependence of the estimates on the degree p of the polynomials on each element. This will be the subject of future work. In addition, because of the use of a specific estimate for the dual problem, we need to assume that the domain is convex. Obviously this restriction needs to be lifted in the future. The plan of the paper is as follows. In the next section (Section 2) we state precisely the problem to be solved and define the HDG methods. We also collect some results from [6] used in our paper. In Section 3 we derive stability (a Gårding-type equality), consistency and duality estimates for the Helmholtz problem. The error estimates for the HDG scheme are obtained in Section 4. In Section 5 we present our numerical results, and finally we draw some conclusions.
2
The hybridizable discontinuous Galerkin method
We start by describing in detail the model problem we shall analyze: the interior Dirichlet problem for the Helmholtz equation. We denote by Ω ⊂ Rn , n = 2, 3, a Lipschitz polyhedral domain with boundary ∂Ω and
2
seek the complex valued solution u of the following problem ∆u + κ2 u = f˜ u=g
in Ω,
(1a)
on ∂Ω,
(1b)
with wave number κ > 0, source term f˜ ∈ L2 (Ω), and boundary values g ∈ H 1/2 (∂Ω). The HDG scheme is based on a first order formulation of the above problems so we introduce q := (i /κ)∇u and f := (i f˜)/κ. Then the Helmholtz equation (1) can be rewritten in mixed form as finding (u, q) such that i κq + ∇u = 0
in Ω,
(2a)
i κu + div q = f
in Ω,
(2b)
on ∂Ω.
(2c)
u=g
Assuming that κ2 is not a Dirichlet eigenvalue for Ω, existence and uniqueness of solutions u ∈ H 1 (Ω) to (1) (or equivalently of solutions (q, u) ∈ H (div, Ω) × H 1 (Ω) to (2)) is well known and these solutions satisfy standard elliptic regularity results, for example kqkH (div,Ω) + kukH 1 (Ω) ≤ C(kf kL2 (Ω) + kgkH 1/2 (∂Ω) )
(3)
where C depends on κ and Ω. Here and throughout C denotes a generic constant, the value of which might be different at different occurrences. The development of the HDG scheme starts by meshing the domain. We consider a subdivision of Ω into a finite element mesh of regular tetrahedra K with faces F in R3 (or triangles K and edges F in R2 ) and denote the resulting collection of triangles/tetrahedra by Th , the collection of edges/faces by Eh , and the collection of element boundaries by ∂Th := {∂K | K ∈ Th }. Regularity here means that the elements are non-degenerate [2]. On each element K and each edge/face F we consider the local spaces of polynomials of degree smaller or equal to k ≥ 0, V (K) := P k (K) := (Pk (K))n ,
W (K) := Pk (K),
and M (F ) := Pk (F ),
where Pk (S) denotes the space of polynomials of total degree at most k in m variables on S where S ⊂ Rm . The corresponding global finite element spaces are given by V h := v ∈ L2 (Ω) v|K ∈ V (K) for all K ∈ Th , Wh := w ∈ L2 (Ω) w|K ∈ W (K) for all K ∈ Th , Mh := µ ∈ L2 (Eh ) µ|F ∈ M (F ) for all F ∈ Eh }, Q where L2 (Ω) := (L2 (Ω))n and L2 (Eh ) := F ∈Eh L2 (F ). Similarly, we will write H s (Ω) := (H s (Ω))n for any s > 0. On these spaces we use the bilinear forms X X X (v, w)Th := (v, w)K , (v, w)Th := (v, w)K , and hv, wi∂Th := hv, wi∂K , K∈Th
R
K∈Th
K∈Th
R
R
with (v, w)K := K v · w dx, (v, w)K := K vw dx, and hv, wi∂K := ∂K vw ds. The hybridizable discontinuous Galerkin method considered in this work belongs to the class of hybridizable local discontinuous Galerkin (LDG-H) methods (cf., e.g., [4]) and yields approximations uh ∈ Wh of u, q h ∈ V h of q, and a numerical trace u bh ∈ Mh approximating u on Eh , which satisfy uh , r · ni∂Th = 0, (i κq h , r)Th − (uh , div r)Th + hb
(4a)
q h · n, wi∂Th = (f, w)Th , (i κuh , w)Th − (q h , ∇w)Th + hb
(4b)
hb uh , µi∂Ω = hg, µi∂Ω , hb q h · n, µi∂Th \∂Ω = 0,
(4c) (4d)
3
for all r ∈ V h , w ∈ Wh , and µ ∈ Mh where the overbar denotes complex conjugation. Here the numerical bh is given by flux q bh = qh + τ (uh − u q bh )n on ∂Th (4e) for some non-negative stabilization function τ : ∂Th → R that is assumed to be constant on each edge/face F ∈ Eh , and n is the outward normal to K ∈ Th . From (2a) we observe that u and q are dimensionally equivalent, and we will show in the next section that this definition yields a consistent scheme. The error analysis for this method presented in the following two sections consists of a variation and extension of arguments used in [6] to analyze a corresponding LDG-H method for a second-order uniformly elliptic problem. It is based on a projection operator Πh inspired by the particular form of the numerical bh from (4e). This projection is defined by flux q Πh : L2 (Ω) × L2 (Ω) ⊃ dom(Πh ) → V h × Wh ,
Πh (q, u) := (ΠV q, ΠW u),
where for any K ∈ Th and all edges/faces F of K the functions ΠV q and ΠW u satisfy (ΠV q, v)K = (q, v)K
for all v ∈ P k−1 (K),
(5a)
(ΠW u, w)K = (u, w)K
for all w ∈ Pk−1 (K),
(5b)
for all µ ∈ Pk (F ).
(5c)
hΠV q · n + τ ΠW u, µiF = hq · n + τ u, µiF
The domain dom(Πh ) of Πh is such that the right hand sides of (5) are well defined, and we note that by (3) solutions to (2) belong to dom(Πh ). max := max τ |∂K > 0, it has been shown in Using the assumed regularity of the mesh, and denoting by τK [6, Thm. 2.1] that the system (5) is uniquely solvable and that its solution (ΠV q, ΠW u) has the following approximation properties: There exists a constant C independent of K, κ and τ such that l +1
kΠV q − qkL2 (K) ≤ ChKq
∗ |q|H lq +1 (K) + ChlKu +1 τK |u|H lu +1 (K) ,
(6a)
l +1
kΠW u − ukL2 (K) ≤ ChlKu +1 |u|H lu +1 (K) + C
hKq max | div q|H lq (K) , τK
(6b)
∗ := max τ |∂K\F ∗ with F ∗ being an edge/face of K at which τ |∂K is maximum. for lu , lq ∈ [0, k], where τK We will also use the following weak commutativity property which is proved in Proposition 2.1 of [6]. For any w ∈ Wh and any (Φ, Ψ) ∈ dom(Πh ),
(w, div Φ)K = (w, div ΠV Φ)K + hw, τ (ΠW Ψ − Ψ)i∂K .
3
(7)
Consistency, Gårding-type identity, and duality
This section is devoted to analyzing the key properties of the scheme that allow us to prove the error estimates in the next section. The arguments in this section rest heavily on those in [6] suitably modified for the non-coercive Helmholtz equation and complex valued solution functions. We denote by PM the L2 -orthogonal projection onto Mh and note that since τ is piecewise constant,
τ (PM u − u), µ ∂T = 0 for all µ ∈ Mh . (8) h
Lemma 3.1 (Consistency). Introducing the projected errors εqh := ΠV q − q h , εuh := ΠW u − uh , and u εb bh , we have h := PM u − u u (i κεqh , r)Th − (εuh , div r)Th + hεb h , r · ni∂Th = (i κ(ΠV q − q), r)Th ,
(i κεuh , w)Th
−
(εqh , ∇w)Th
+ hb εh · n, wi∂Th = (i κ(ΠW u − u), w)Th . u hεb h , µi∂Ω = 0, hb εh · n, µi∂Th \∂Ω = 0, 4
(9a) (9b) (9c) (9d)
for all r ∈ V h , w ∈ Wh , and µ ∈ Mh . Here, u b bh · n εh · n := εqh · n + τ (εuh − εb h ) = PM (q · n) − q
on ∂Th \ ∂Ω.
(9e)
Proof. Recalling that the exact solution (q, u) of (2) satisfies (i κq, r)Th − (u, div r)Th + hu, r · ni∂Th = 0, (i κu, w)Th − (q, ∇w)Th + hq · n, wi∂Th = (f, w)Th for all r ∈ V h and w ∈ Wh , the definitions of Πh and PM yield that (i κΠV q, r)Th − (ΠW u, div r)Th + hPM u, r · ni∂Th = (i κ(ΠV q − q), r)Th (i κΠW u, w)Th − (ΠV q, ∇w)Th + hΠV q · n + τ (ΠW u − u), wi∂Th = (f, w)Th + (i κ(ΠW u − u), w)Th for all r ∈ V h and w ∈ Wh . Now subtracting (4a) from the first equation we obtain (i κ(ΠV q − q h ), r)Th − (ΠW u − uh , div r)Th + hPM u − u ˆh , r · ni∂Th = (i κ(ΠV q − q), r)Th , which yields (9a). Similarly, subtracting (4b) from the second equation above shows that bh ) · n + τ (ΠW u − u), wi∂Th = (i κ(ΠW u − u), w)Th . (i κ(ΠW u − uh ), w)Th − (ΠV q − q h , ∇w)Th + h(ΠV q − q bh in (4e), we find that Recalling the definition of b εh in (9e) and of the numerical flux q b εh · n = (ΠV q − q h ) · n + τ (ΠW u − uh − PM u + u bh ) bh ) · n + τ (uh − u = (ΠV q − q bh ) + τ (ΠW u − u) − τ (uh + PM u − u bh − u), which together with (8) yields (9b). To show (9c) we observe that (4c), (8), and (2c) imply u hεb h , µi∂Ω = hPM u − g, µi∂Ω = hu − g, µi∂Ω = 0.
Finally, using the definition of b εh · n in (9e), (5c), (8), and (4e) we obtain hb εh · n, µi∂Th \∂Ω = h(ΠV q − q h ) · n + τ (ΠW u − uh − PM u + u bh ), µi∂Th \∂Ω = h(q − q h ) · n + τ (u − uh − u + u bh ), µi∂Th \∂Ω q h · n, µi∂Th \∂Ω . = hq · n, µi∂Th \∂Ω − hb bh · n ∈ Mh . Furthermore, since q ∈ H (div, Ω) and q bh satisfies This proves (9e) because b εh · n ∈ Mh and q (4d), both terms on the right hand side vanish, which implies (9d). The following Gårding-type identity is a little more complex than the coercivity bound for the Poisson problem. In addition compared to conforming methods we see the presence of coupling terms on the boundary. From this estimate we see that choosing τ to be real valued gives well controlled terms on element boundaries. Lemma 3.2 (Gårding-type identity). The projected errors satisfy
q u u u b u i κ (εqh , εqh )Th − (εuh , εuh )Th + τ (εuh − εb h ), εh − εh ∂T = i κ(ΠV q − q), εh T + i κ(ΠW u − u), εh T . h
εqh
h
h
εuh
Proof. If we choose r = in (9a), w = in the complex conjugate of (9b), µ = −b εh · n in (9c), and b u µ = −εh in the complex conjugate of (9d) and add these equations we obtain that q u q u u u (i κεqh , εqh )Th − (εuh , div εqh )Th + hεb εh · n, εuh i∂Th h , εh · ni∂Th + (i κεh , εh )Th − (εh , ∇εh )Th + hb q u u u − hεb εh · ni∂Ω − hb εh · n, εb h, b h i∂Th \∂Ω = (i κ(ΠV q − q), εh )Th + (i κ(ΠW u − u), εh )Th .
5
Therefore, (i κεqh , εqh )Th + (i κεuh , εuh )Th + Rh = (i κ(ΠV q − q), εqh )Th + (i κ(ΠW u − u), εuh )Th , with q u q u u εh · n, εuh i∂Th − hεb εh · ni∂Th . Rh = −(εuh , div εqh )Th + hεb h , εh · ni∂Th − (εh , ∇εh )Th + hb h, b
Integrating by parts and applying (9e) yields q u q u u Rh = (∇εuh , εqh )Th − hεuh , εqh · ni∂Th + hεb εh · n, εuh − εb h , εh · ni∂Th − (εh , ∇εh )Th + hb h i∂Th q u u i + hb ε · n, εu − εb i = −hε · n, εu − εb h
h
h ∂Th
h
h
h ∂Th
u u ), εu − εb = hτ (εuh − εb h h i∂Th . h
To obtain error estimates for q h and uh in Section 4 we use a duality argument. To this end, given Θ ∈ L2 (Ω), we introduce the dual problem − i κΦ + ∇Ψ = 0
in Ω,
(10a)
− i κΨ + div Φ = Θ
in Ω,
(10b)
on ∂Ω
(10c)
Ψ=0
and note that (10) is uniquely solvable since we assumed that κ is not a Dirichlet eigenvalue for Ω. Lemma 3.3 (Duality). Given Θ ∈ L2 (Ω) we denote by (Φ, Ψ) the solution to (10). Then, for any Ψh ∈ Wh and any ψh ∈ Wh such that ψh |K ∈ Pk−1 (K) for all K ∈ Th , (εuh , Θ)Th = i κ(q − q h ), ΠV Φ − Φ T + ΠV q − q, ∇Ψ − ∇Ψh T h (11) h + i κ(ΠW u − u), ΠW Ψ − ψh T − i κεuh , ΠW Ψ − Ψ T . h
h
Proof. Using (10b), and the weak commutativity property (7), we obtain (εuh , Θ)Th = (εuh , − i κΨ + div Φ)Th = (εuh , − i κΨ)Th + (εuh , div ΠV Φ)Th + hεuh , τ (ΠW Ψ − Ψ)i∂Th Hence using and (9a) this equality can be rewritten as u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κεqh , ΠV Φ)Th + hεb h , ΠV Φ · ni∂Th − (i κ(ΠV q − q), ΠV Φ)Th
+ hεuh , τ (ΠW Ψ − Ψ)i∂Th u u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hεb h , (ΠV Φ − Φ) · ni∂Th + hεh , τ (ΠW Ψ − Ψ)i∂Th , u where we used in the last step that Φ · n is continuous across all edges/faces and that εb h = 0 on ∂Ω by (9d). Then, using (5c) and (8) this becomes u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hεuh − εb h , τ (ΠW Ψ − Ψ)i∂Th u u b u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th − hτ (εh − εh ), PM Ψi∂Th .
Using the discrete error equations (9e), and (9d), the boundary condition (10c), and the properties of PM in (8), we obtain q u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th + hεh · n, PM Ψi∂Th q u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th + hεh · n, Ψi∂Th .
6
Applying (9e), (9b), and integration by parts, we find that εh · n − εqh · n, ΠW Ψi∂Th + hεqh · n, Ψi∂Th (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hb = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (i κεuh , ΠW Ψ)Th + (εqh , ∇ΠW Ψ)Th − hεqh · n, ΠW Ψi∂Th + hεqh · n, Ψi∂Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (i κεuh , ΠW Ψ)Th − (div εqh , ΠW Ψ)Th + hεqh · n, Ψi∂Th . Next using (5b), and integration by parts once more, we obtain (εuh , Θ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (div εqh , Ψ)Th + hεqh · n, Ψi∂Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (εqh , ∇Ψ)Th . Finally, we use (10a), (5a), and (5b) to obtain (εuh , Θ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th + (i κ(q − q h ), Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (ΠV q − q h , ∇Ψ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th − (q − q h , ∇Ψ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (ΠV q − q h , ∇Ψ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th + (ΠV q − q, ∇Ψ − ∇Ψh )Th + (i κ(ΠW u − u), ΠW Ψ − ψh )Th which is the desired identity.
4
The error estimates
Before we state the first error estimates in Theorem 4.1 below,Qwe introduce the mesh parameter h := maxK∈Th hK and the mesh dependent norm k · kh on L2 (∂Th ) := K∈Th L2 (∂K) by X kµk2h := hK kµk2∂K for any µ ∈ L2 (∂Th ). K∈Th
Theorem 4.1 (Flux error estimates). Assume k ≥ 1. Given f ∈ L2 (Ω) and g ∈ H 1/2 (∂Ω) we denote by q and u the exact solution of (2), and by q h , uh , and u bh the solution of (4). Furthermore, for any Θ ∈ L2 (Ω) let (Φ, Ψ) be the solution of the dual problem (10) and suppose that kΦkH 1 (Ω) + kΨkH 2 (Ω) ≤ Cκ kΘkL2 (Ω) , where the constant Cκ is independent of Θ but depends on the wave number κ. If 1 ∗ Cκ κhMτmin := Cκ κh max 1, and Cκ κhMτmax := Cκ κh max(1, max hK τK ) max K∈Th minK∈Th τK
(12)
(13)
are sufficiently small, then kΠV q − q h kL2 (Ω) ≤ C (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω)
(14)
and bh · nkh kPM (q · n) − q
κ max 1/2 ≤ C max 1, max (hK τK ) (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) , (15) c K∈Th
where c and C are independent of κ and h for hκ small enough. 7
Note that the regularity assumption (12) on the solution of the dual problem used in Theorem 4.1 is, e.g., well known to hold if Ω is convex (see Brenner and Scott [1, p 139]). Proof. Substituting Θ = εuh into (10) we find from (11) that for any Ψh ∈ Wh and any ψh ∈ Wh such that ψh |K ∈ Pk−1 (K) for all K ∈ Th , (εuh , εuh )Th = (i κ(q − q h ), ΠV Φ − Φ)Th + (ΠV q − q, ∇Ψ − ∇Ψh )Th + (i κ(ΠW u − u), ΠW Ψ − ψh )Th − (i κεuh , ΠW Ψ − Ψ)Th . So, using the Cauchy-Schwartz inequality we obtain that kεuh k2L2 (Ω) ≤ C κkq − q h kL2 (Ω) kΠV Φ − ΦkL2 (Ω) + kΠV q − qkL2 (Ω) k∇Ψ − ∇Ψh kL2 (Ω) + κkΠW u − ukL2 (Ω) kΠW Ψ − ψh kL2 (Ω) + κkεuh kL2 (Ω) kΠW Ψ − ΨkL2 (Ω) .
(16)
The first term on the right hand side of (16) can be estimated by means of (6a) with lu = 1 and lq = 0, ∗ κkq − q h kL2 (Ω) kΠV Φ − ΦkL2 (Ω) ≤ Cκhkq − q h kL2 (Ω) |Φ|H 1 (Ω) + max hK τK |Ψ|H 2 (Ω) , K∈Th
Similarly, using the projection error bound |Ψ − Ψh |H 1 (Ω) ≤ Ch|Ψ|H 2 (Ω) we find for the second term that kΠV q − qkL2 (Ω) k∇Ψ − ∇Ψh kL2 (Ω) ≤ ChkΠV q − qkL2 (Ω) |Ψ|H 2 (Ω) . To estimate the third term in (16), we apply the triangle inequality, (6b) with lu = 0 and lq = 0, and the projection error bound kΨ − ψh kL2 (Ω) ≤ Ch|Ψ|H 1 (Ω) , which yields κkΠW u − ukL2 (Ω) kΠW Ψ − ψh kL2 (Ω) ≤ κkΠW u − ukL2 (Ω) kΠW Ψ − ΨkL2 (Ω) + kΨ − ψh kL2 (Ω) 1 | div Φ|L2 (Ω) . ≤ CκhkΠW u − ukL2 (Ω) |Ψ|H 1 (Ω) + max minK∈Th τK Similarly, we find for the last term on the right hand side of (16) that 1 2 (Ω) . | div Φ| κkεuh kL2 (Ω) kΠW Ψ − ΨkL2 (Ω) ≤ Cκhkεuh kL2 (Ω) |Ψ|H 1 (Ω) + L max minK∈Th τK Combining these estimates using (12) with Θ = εuh we obtain that (16) reduces to kεuh k2L2 (Ω) ≤ CCκ κhMτmax kq − q h kL2 (Ω) kεuh kL2 (Ω) + CCκ hkΠV q − qkL2 (Ω) kεuh kL2 (Ω) + CCκ κhMτmin kΠW u − ukL2 (Ω) kεuh kL2 (Ω) + CCκ κhMτmin kεuh k2L2 (Ω) . Therefore, assuming that Cκ κhMτmin is small enough, 1 kεuh kL2 (Ω) ≤ CCκ κh Mτmax kεqh kL2 (Ω) + + Mτmax kΠV q − qkL2 (Ω) + Mτmin kΠW u − ukL2 (Ω) . κ
(17)
Now recalling Lemma 3.2, using the arithmetic-geometric mean inequality with some α > 0, we find from (17) that i q q q u u u u u b u (εh ,εh )Th − hτ (εuh − εb h ), εh − εh i∂Th = (ΠV q − q), εh )Th − ((ΠW u − u), εh )Th + (εh , εh )Th κ α 1 1 3 ≤ kεqh k2L2 (Ω) + kΠV q − qk2L2 (Ω) + kΠW u − uk2L2 (Ω) + kεuh k2L2 (Ω) 2 2α 2 2 α 1 1 q ≤ + CCκ2 κ2 h2 (Mτmax )2 kεh k2L2 (Ω) + + CCκ2 κ2 h2 2 + (Mτmax )2 kΠV q − qk2L2 (Ω) 2 2α κ 1 + + CCκ2 κ2 h2 (Mτmin )2 kΠW u − uk2L2 (Ω) . 2 8
Therefore, choosing α small enough, assuming that Cκ κhMτmax is small, and recalling that Cκ κhMτmin is small, we obtain that c u 2 2 2 2 u b u (18) kεqh k2L2 (Ω) + hτ (εuh − εb h ), εh − εh i∂Th ≤ C (1 + Cκ h )kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) , κ for some constants c > 0 and C > 0, which implies (14). To show (15), we recall (9e) and use an inverse inequality (element by element) to obtain u bh · nkh ≤ kεqh · nkh + kτ (εuh − εb kPM (q · n) − q h )kh max 1/2 u u b u 1/2 ≤ Ckεqh kL2 (Ω) + max (hK τK ) hτ (εuh − εb h ), εh − εh i∂Th K∈Th κ c u max 1/2 u − εb u i1/2 . ≤ C max 1, max (hK τK kεqh kL2 (Ω) + hτ (εuh − εb ), ε ) h h h ∂Th c K∈Th κ
Therewith, (15) follows from (18). Theorem 4.2 (Convergence of the scalar variable). Suppose that the assumptions of Theorem 4.1 hold. If Cκ κhMτmin
and
Cκ κhMτmax
(as defined in (13)) are sufficiently small, then kΠW u − uh kL2 (Ω) ≤ CCκ κh
1 κ
+ Mτmax kΠV q − qkL2 (Ω) + (Mτmin + Mτmax )kΠW u − ukL2 (Ω)
(19)
and kPM u − u bh kh ≤ Cκh
1 + Cκ h + Cκ
1 κ
+ Mτmax
kΠV q − qkL2 (Ω)
+ 1+
Cκ (Mτmin
+
Mτmax )
kΠW u − uk
L2 (Ω)
. (20)
Remark 4.3. Collecting together the results of Theorems 4.1 and 4.2 together with the projection estimates (6a) and (6b) we can obtain explicit bounds for the error in terms of the mesh size. A practical case is when τ = O(1) and u ∈ H k+2 (Ω) (convergence for less regular solutions follows in the obvious way). Then we obtain ku − uh kL2 (Ω) + kq − q h kL2 (Ω) = O(hk+1 ) . Similarly, Theorem 4.2 guarantees super-convergence of the projected error of the scalar variable and of the numerical trace, which means for τ = O(1) and u ∈ H k+2 (Ω) that kΠW u − uh kL2 (Ω) = O(hk+2 )
and
kPM u − u bh kh = O(hk+2 ) ,
respectively. Proof. Combining (17) and (14), and recalling that Cκ κhMτmax is small, we find that kΠW u − uh kL2 (Ω) ≤ CCκ κh Mτmax (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) 1 + + Mτmax kΠV q − qkL2 (Ω) + Mτmin kΠW u − ukL2 (Ω) κ 1 ≤ CCκ κh (1 + κMτmax )kΠV q − qkL2 (Ω) + (Mτmin + Mτmax )kΠW u − ukL2 (Ω) , κ which is (19). 9
u For (20) we follow the proof of [6, Thm. 4.1] and select r ∈ P k (K) such that r · n = εb h on ∂K and 1/2 b u krkL2 (K) ≤ ChK kεh kL2 (∂K) . Now substituting hK r into (9a) and applying an inverse inequality (see [1, Lmm. 4.5.3]) we obtain q u u 2 hK kεb h kL2 (∂K) = hK (i κ(ΠV q − q), r)K + hK (εh , div r)K − hK (i κεh , r)K ≤ CκhK krkL2 (K) kΠV q − qkL2 (K) + kεqh kL2 (K) + CkrkL2 (K) kεuh kL2 (K) 1/2 u q u ≤ ChK kεb h kL2 (∂K) κhK kΠV q − qkL2 (K) + κhK kεh kL2 (K) + kεh kL2 (K) .
Therefore, using (19) and (14) we find that kPM u − u bh k2h ≤ C κ2 h2 kΠV q − qk2L2 (Ω) + κ2 h2 kεqh k2L2 (Ω) + kεuh k2L2 (Ω) ≤ C κ2 h2 kΠV q − qk2L2 (Ω) + κ2 h2 (1 + Cκ h)2 kΠV q − qk2L2 (Ω) + kΠW u − uk2L2 (Ω) +
Cκ2 κ2 h2
≤ Cκ2 h2
1 κ
+
Mτmax
1 + Cκ2 h2 + Cκ2
2
1
kΠV q −
+ Mτmax
qk2L2 (Ω)
+
(Mτmin
+
Mτmax )2 kΠW u
−
uk2L2 (Ω)
2
kΠV q − qk2L2 (Ω) 2 min max 2 2 + 1 + Cκ (Mτ + Mτ ) kΠW u − ukL2 (Ω) , κ
which yields (20). Remark 4.4. The super-convergence of the approximation uh of the scalar variable u to the projection ΠW u established in Theorem 4.2 can be exploited to obtain an approximation u∗h of improved accuracy. As outlined in section 5 of [6], the postprocessing schemes for mixed finite elements proposed by Stenberg [11] can be suitably adapted and applied to the hybridizable DG method (see also [3] and Cockburn, Guzmán, and Wang [5]). Following the analysis in [11] it can be shown that (for k ≥ 1) the postprocessed solution u∗h converges of one order higher than uh .
5
Numerical results
In this section we probe the error estimates established in Theorem 4.1 and Theorem 4.2 by numerical examples for convex and non-convex two-dimensional domains, and for solutions of differing regularity. The HDG-method has been implemented for piecewise constant, piecewise linear, and piecewise quadratic finite element spaces in R2 . We use order 6 quadrature rules on the triangles as well as on the edges to evaluate the system matrix and the right hand side corresponding to (4). The mesh generation and all computations are done in MATLAB. Example 5.1. In the first example we consider a square domain Ω =]−1, 1[×]−1, 1[⊂ R2 and structured meshes Th with N = 8, 32, 128, 512, and 2048 triangles as shown in Figure 1 (left) for N = 128. We choose the boundary conditions g and the source term f in (4) such that the exact solution of (1) is given by p u(x) = Jξ κ (x + 1)2 + y 2 cos ξ arctan(y/(x + 1)) for ξ = 1, ξ = 2/3, and ξ = 3/2, respectively, where Jξ denotes the Bessel function of the first kind and order ξ. It is well known (see, e.g., Hiptmair, Moiola, and Perugia [8] and Grisvard [7] for details) that u(x) is smooth for ξ ∈ N, while its derivative has a singularity at (−1, 0) for ξ ∈ / N. More precisely, u ∈ H 2 (Ω) 3 1 2 but u ∈ / H (Ω) for ξ = 3/2 and u ∈ H (Ω) but u ∈ / H (Ω) for ξ = 2/3. 10
1
1
0.5
0.5
0
0
ï0.5
ï0.5
ï1 ï1
ï0.5
0
0.5
ï1 ï1
1
ï0.5
0
0.5
1
Figure 1: Visualization of structured meshes used in Example 5.1 (left) and Example 5.2 (right). Figure 2 shows the absolute error in the scalar variable ku − uh kL2 (Ω) , the absolute error in the flux variable kq − q h kL2 (Ω) , and the absolute error in the numerical trace ku − u bh kh as a function of κh/2π for ξ = 1, ξ = 2/3, and ξ = 3/2 (top down) and wave number κ = 4, i.e., the wave length is λ = 2π/κ ≈ 1.57. The plots in the left column correspond to k = 1, i.e., piecewise linear approximations, while the right column has been obtained for k = 2, i.e., piecewise quadratic approximations. As predicted by (14), (19), and (20) together with (6) we observe optimal convergence in qh , uh , and u bh for the smooth solution, i.e., ξ = 1, using linear as well as quadratic elements. The singular solutions with ξ = 2/3 and ξ = 3/2 show a slightly different behavior: For linear elements and ξ = 2/3, i.e., u ∈ H 1 (Ω) but u ∈ / H 2 (Ω), we loose one order of convergence in the flux variable, while the scalar variables still converge at optimal order. On the other hand, for quadratic elements and ξ = 2/3 the convergence in all three variables kq − q h kL2 (Ω) , ku − uh kL2 (Ω) , and ku − u bh kh slows down. If ξ = 3/2, i.e., u ∈ H 2 (Ω) but u ∈ / H 3 (Ω), the convergence is optimal in all three variables for the piecewise linear approximation while the convergence in the flux variable slows down for the piecewise quadratic approximation. This behavior is actually better than our theoretical results predict: By (6) our estimates guarantee optimal convergence of ku − uh kL2 (Ω) only if the exact solution u is in H k+1 (Ω) and optimal convergence order in kq − q h kL2 (Ω) only if u ∈ H k+2 (Ω). For the piecewise constant approximation, i.e., k = 0, we observed (results not shown) optimal convergence for all three choices of ξ and sufficiently small h after a comparatively large pre-asymptotic area of slower convergence. Next we study the dependence of the convergence rate on the wave number κ. To this end we repeat the previous experiment for ξ = 1, quadratic elements, and four different wave numbers κ = 1, κ = 4, κ = 10, and κ = 16. Figure 3 shows ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh as a function of κh/2π. In all four examples we observe optimal convergence order for sufficiently small values of h. With increasing wave number the absolute error increases and a pre-asymptotic region of slower convergence for large values of h appears, as is to be expected because of pollution error in the h version. If we order the degrees of freedom counting first those for uh and q h and then those for u ˆh we obtain a block linear system A B x b = y c BT C where x denotes the vector of degrees of freedom for uh and q h while y contains the degrees of freedom for the trace u ˆh . Because of the use of a discontinuous Galerkin method, the matrix A is block diagonal (assuming degrees of freedom are ordered element by element). Thus assuming it is invertible (numerically the case in our study) we can cheaply eliminate x and obtain the Schur complement system Sy = d where S = (C − B T A−1 B) and d = c − B T A−1 b. Of course in practical calculations the matrix S would be computed element by element so that it would never be necessary to assemble A or B. The Schur 11
1
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
error
error
10 ï1
10
ï2
10
ï3
10 ï2
10
ï4
10 ï3
10
ï5
0.05
0.1
10
0.5
0.05
0.1
gh/2/ 1
1
10
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3) O(h)
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
0
10
ï1
error
10 error
0.5 gh/2/
ï1
10
ï2
10
ï3
10 ï2
10
ï4
10 ï3
10
ï5
0.05
0.1
10
0.5
0.05
0.1
gh/2/
0.5 gh/2/
1
1
10
10
2
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
10 ï1
error
error
10
ï2
10
ï2
10
ï3
10 ï3
10
ï4
10
ï4
10
ï5
0.05
0.1
10
0.5
0.05
gh/2/
0.1
0.5 gh/2/
Figure 2: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, ξ = 2/3, and ξ = 3/2 (top down) and wave number κ = 4 for linear elements (left) and quadratic elements (right).
12
1
2
10
10
0
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
10 ï2
error
error
10
ï2
10
ï4
10
ï3
10 ï6
10
ï4
10
ï8
ï5
10
0.05 gh/2/
10
0.1
2
3
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
2
10
0
error
error
10
ï1
10
ï2
0
10
ï1
10
10
ï3
ï2
10
10
ï4
0.1
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
1
10
10
0.5
10 2
10
0.1 gh/2/
10
1
0.05
ï3
10
0.5 gh/2/
0.5 gh/2/
Figure 3: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1 and κ = 1, κ = 4, κ = 10, and κ = 16 (from left to right and top down) for quadratic elements.
13
κ=1 κ=4 κ = 10 κ = 16
N =8 2.98e+02 9.76e+01 1.11e+02 6.14e+01
N = 32 9.84e+02 3.98e+02 2.39e+02 2.95e+02
N = 128 2.91e+03 1.44e+03 3.64e+03 6.70e+02
N = 512 9.55e+03 5.60e+03 1.75e+04 1.67e+04
N = 2048 3.40e+04 2.31e+04 7.82e+04 1.37e+05
Table 1: Condition numbers of the Schur complement matrix S with respect to the wave number κ and the number of triangles N for quadratic elements. 1
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
10
ï1
ï1
10 error
error
10
ï2
10
ï3
ï2
10
ï3
10
10
ï4
ï4
10
10
ï5
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
ï5
0.05
0.1
10
0.5 gh/2/
0.05
0.1
0.5 gh/2/
Figure 4: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, wave number κ = 4, and τ = h, τ = 1, and τ = 1/h (from left to right) for quadratic elements.
complement matrix S (equivalently, the matrix for the condensed hybridized system) is sparse but the system would typically be solved iteratively in practice so it is of interest to examine the condition number of S as a function of the number of triangles N and wave number κ as shown in Table 1 for quadradic elements. In the examples considered so far the coupling parameter τ from (4e) has been chosen to be τ = 1. To analyze the dependence of the numerical solution on this parameter, we show in Figure 4 the convergence behavior for the smooth solution with ξ = 1 and wave number κ = 4 for τ = h and τ = 1/h. These plots should be compared with the corresponding plot for τ = 1 in Figure 2 (top right). While we observed optimal convergence for τ = 1, the convergence of ku − uh kL2 (Ω) slows down for τ = h, while for τ = 1/h the convergence speed of the flux variable kq − q h kL2 (Ω) is reduced. This reflects exactly the approximation properties of the projection operator Πh stated in (6) that have been used in the convergence analysis. The behavior for τ = 1/h clearly follows from (14), (19), and (20), but this is not so obvious for τ = h. To exemplify the super-convergence properties of the HDG-scheme established in Theorem 4.2, we show in Figure 5 plots of kPM u − u bh kh corresponding to the smooth solution with ξ = 1, wave number κ = 4, and τ = 1 for linear elements (left) and quadratic elements (right). Super-convergence is clearly observed: As predicted in (20) the convergence order of kPM u − u bh kh is O(h3 ) for linear elements and O(h4 ) for quadratic elements. Example 5.2. In the second example we consider a non-convex domain as shown in Figure 1 (right). The boundary conditions g and the source term f in (4) are chosen similarly to the first example such that the exact solution of (1) is given by p u(x) = Jξ κ x2 + y 2 cos ξ arctan(y/x))
14
1
0
10
10 L2ïerror of u L2ïerror of q projected hnormïerror of uhat O(h2) O(h3)
0
10
L2ïerror of u L2ïerror of q projected hnormïerror of uhat O(h3) O(h4)
ï1
10
ï2
10
ï1
error
error
10
ï3
10
ï2
10
ï4
10 ï3
10
ï5
10
ï4
10
ï6
0.05
0.1
10
0.5 gh/2/
0.05
0.1
0.5 gh/2/
Figure 5: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and kPM u − u bh kh with ξ = 1, wave number κ = 4, and τ = 1 for linear elements (left) and quadratic elements (right).
for ξ = 1, ξ = 2/3, and ξ = 3/2. The smoothness properties of these solutions are exactly the same as in the first example, but for ξ = 2/3, and ξ = 3/2 the singularity in the derivative of the exact solution is now at (0, 0). Due to the non-convexity of this domain, Theorem 4.1 and Theorem 4.2, which have both been obtained using duality arguments, do not even apply for the smooth solution. Nevertheless, we compute piecewise quadratic numerical approximations of the solution to (1) with wave number κ = 4 and coupling parameter τ = 1 for regular structured meshes with N = 24, 96, 384, and 1536 elements (see Figure 1 (right) for a visualization with N = 96). Plots of the corresponding numerical errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh for ξ = 1, ξ = 2/3, and ξ = 3/2 can be found in Figure 6. We observe the same qualitative behavior as for the convex domain: For ξ = 1 we have optimal convergence in all three variables, while if ξ = 3/2 the convergence of the flux variable slows, and if ξ = 2/3 the flux variable, the scalar variable and the numerical trace no longer converge at optimal order. Since we assume that the convergence order for ξ = 3/2 and ξ = 2/3 degrades because of to the singularity of the derivative of the exact solution at (0, 0), we consider in our final experiment a sequence of unstructured locally refined meshes with N = 197, 254, 417, and 629 elements as shown in Figure 7. Corresponding graphs of the absolute errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh as a function of the number of triangles N can be found for κ = 4, quadratic elements, and ξ = 2/3 (left) as well as for ξ = 3/2 (right) in Figure 8. Indeed the errors decrease when we refine the grid locally around the singularity. However convergence is slower for the smoother (ξ = 3/2) solution in the right panel compared to that in the left panel (ξ = 2/3). This is because the locally refined grids are chosen for the more singular problem and are under refined away from the singularity when ξ = 3/2.
Concluding Remarks We have show that the HDG scheme of [4] can be applied to the Helmholtz equation. In particular techniques from [6] can be used to prove error estimates, and our numerical results demonstrate that the method can be applied even when the restrictive assumptions of our results do not hold. Indeed the numerical results suggest that it should be possible to prove error estimates for non-convex domains and improve on the regularity requirements of the solution in order to obtain a particular rate. Clearly the main open problem is to prove h − p convergence results, since in practice we need to adjust p and h in order to control dispersion element by element (see Melenk and Sauter [9]). Our numerical results suggest that HDG methods have great promise to provide a flexible and efficient technique for solving linear time harmonic wave propagation problems.
15
0
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
ï1
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3) O(h)
0
10
ï2
10
ï1
error
error
10 ï3
10
ï2
10 ï4
10
ï3
10
ï5
10
ï6
ï4
10
0.05
10
0.1 gh/2/
0.05
0.1 gh/2/
0
10
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
ï1
10
ï2
error
10
ï3
10
ï4
10
ï5
10
0.05
0.1 gh/2/
Figure 6: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, ξ = 2/3, and ξ = 3/2 (from left to right and top down) for a non-convex domain and quadratic elements.
1
1
1
1
0.5
0.5
0.5
0.5
0
0
0
0
ï0.5
ï0.5
ï0.5
ï0.5
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
Figure 7: Visualization of the locally refined meshes used in Example 5.2.
16
0
0.5
1
ï1
ï1
10
10
ï2
ï2
10
error
error
10
ï3
10
ï4
10
ï5
10
ï4
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(Nï3/2)
200
300 400 number of triangles N
ï3
10
ï5
500
10
600
200
L2ïerror of u L2ïerror of q hnormïerror of uhat O(Nï3/2) 300 400 number of triangles H
500
600
Figure 8: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh for the non-convex domain and locally refined grids with quadratic elements and ξ = 2/3 (left) as well as ξ = 3/2 (right).
Acknowledgements We would like to thank anonymous referees for comments resulting in major improvements to this paper. The research of PM was supported in part by grant number FA9550-08-1-0138 from the US AFOSR.
References [1]
S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, 3rd ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008.
[2]
P. G. Ciarlet, The finite element method for elliptic problems, Studies in Mathematics and its Applications, vol. 4, NorthHolland Publishing Co., Amsterdam-New York-Oxford, 1978.
[3]
B Cockburn, B Dong, and J Guzman, A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems, Math. Comp. 77 (2008), 1887-1916.
[4]
B. Cockburn, J. Gopalakrishnan, and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365.
[5]
B. Cockburn, J. Guzmán, and H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comp. 78 (2009), no. 265, 1–24.
[6]
B. Cockburn, J. Gopalakrishnan, and F.-J. Sayas, A projection-based error analysis of HDG methods, Math. Comp. 79 (2010), no. 271, 1351–1367.
[7]
P. Grisvard, Elliptic problems in nonsmooth domains, Monographs and Studies in Mathematics, vol. 24, Pitman, Boston, MA, 1985.
[8]
R. Hiptmair, A. Moiola, and I. Perugia, Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the p-version, SAM-ETH Zürich Report 2009-20 (2009).
[9]
J. M. Melenk and S. Sauter, Convergence analysis for finite element discretizations of the Helmholtz equation with Dirichletto-Neumann boundary conditions, Math. Comp. 79 (2010), 1871–914.
[10] A. H. Schatz, An observation concerning Ritz-Galerkin methods with indefinite bilinear forms, Math. Comp. 28 (1974), 959–962. [11] R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO Modél. Math. Anal. Numér. 25 (1991), no. 1, 151–167.
17
Peter Monk†
Department of Mathematical Sciences, University of Delaware, Newark, Delaware 19716, U.S.A.
Abstract Finite element methods for acoustic wave propagation problems at higher frequency result in very large matrices due to the need to resolve the wave. This problem is made worse by discontinuous Galerkin methods that typically have more degrees of freedom than similar conforming methods. However hybridizable discontinuous Galerkin methods offer an attractive alternative because degrees of freedom in each triangle can be cheaply removed from the global computation and the method reduces to solving only for degrees of freedom on the skeleton of the mesh. In this paper we derive new error estimates for a hybridizable discontinuous Galerkin scheme applied to the Helmholtz equation. We also provide extensive numerical results that probe the optimality of these results. An interesting observation is that, after eliminating the internal element degrees of freedom, the condition number of the condensed hybridized system is seen to be almost independent of the wave number.
Keywords
discontinuous Galerkin methods, hybridization, error analysis, Helmholtz equation.
AMS subject classification: Short title:
1
65N15, 65N30, 35J05.
HDG for the Helmholtz equation
Introduction
Discontinuous Galerkin (DG) methods have several attractive features when applied to the Helmholtz equation. In particular the polynomial degree can be changed from element to element to compensate for different element sizes in order to control dispersion, and some types of non-standard meshes can be handled. However, they possess an important disadvantage: the dimension of the standard DG piecewise polynomial space is much larger than the dimension of the corresponding conforming space. This increase in dimension is particularly troublesome for wave propagation problems because as the wave-number κ increases (i.e., as the wave length decreases), the dimension of the space increases faster than O(κn ) in n = 2, 3 dimensions and hence becomes very large for high-frequency problems. In addition there is currently, to our knowledge, no optimal solver for the discrete linear system resulting from a finite element discretization of the Helmholtz equation. Since many practical problems involve high-frequency sound (for example simulations of ultra-sound propagation), DG methods appear less practical than conforming elements. Recently, however, hybridizable DG (HDG) methods have been developed that are the focus of this paper (see Cockburn, Gopalakrishnan, and Lazarov [4], Cockburn, Dong and Guzman [3]). HDG methods are DG ∗ e-mail: † e-mail:
[email protected] [email protected] (corresponding author).
1
methods in which it is possible to introduce new variables on the boundary of elements such that the solution inside each element can be computed in terms of these new variables. In particular, element by element, volume degrees of freedom can be cheaply parameterized by the surface degrees. The linear system can then be reduced, at least logically, to solving only for the unknowns on the skeleton of the mesh (i.e. the edges or faces of the mesh depending on the dimension of the space for the Helmholtz equation). Thus the size of the linear system is considerably reduced compared to a standard DG scheme, and objections based on the bloated dimension of DG methods are no longer relevant. This paper is motivated by the analysis of HDG methods via projections due to Cockburn, Gopalakrishnan, and Sayas in [6]. That analysis, applied to the HDG solution of the Poisson problem, shows that the scalar variable in the discrete problem possesses super-convergence to an appropriate projection of the exact solution in the L2 norm. This is attractive from the point of view of the Helmholtz equation. The standard approach to error analysis of the H 1 -conforming finite element approximation of the Helmholtz equation of Schatz [10] uses duality and L2 estimates for the scalar variable. This is then combined with a Gårding inequality to show that, as the mesh size h approaches zero, the conforming finite element approximation approaches the elliptic projection of the solution. The analysis in [6] suggests that a modification of Schatz’ analysis can be used, and we shall indeed show this is the case. Besides proving error estimates, we also provide some numerical results when the method is applied to the interior Dirichlet problem for the Helmholtz equation (away from a Dirichlet eigenvalue, also called a resonance of the domain). Our numerical results show that the predicted asymptotic convergence rates are observed at practical mesh sizes, but suggest that the regularity required of the solution by our theory is too high. We also test the method for non-smooth solutions and on a non-convex domain. The numerical results show that it can be used even when our analysis is not applicable, with the expected rates of convergence. This suggests that optimal convergence can be proved even for non-convex domain as we might expect. An important observation is that, after eliminating the internal element degrees of freedom, the condensed hybridized system for the boundary degrees of freedom has a condition number that is roughly constant as the wave number increases provided we stay away from resonances (until the mesh becomes too coarse to resolve the wave). It may be that hybridized methods offer an approach that will allow the development of wave number robust linear system solvers. The major contributions of our paper are the first h-version analysis of an HDG scheme for the Helmholtz equation which involves also complex valued functions appropriate to this problem. In addition our numerical results point to specific issues related to the optimality of this analysis. Finally our observation concerning the wave number dependence of the condition number of the condensed hybridized system implies that the condensed hybridized system displays better properties with respect to wave number than the global system. Our analysis, like that in [6], is devoted to the h-version of the method. The most desirable extension to our analysis would include the dependence of the estimates on the degree p of the polynomials on each element. This will be the subject of future work. In addition, because of the use of a specific estimate for the dual problem, we need to assume that the domain is convex. Obviously this restriction needs to be lifted in the future. The plan of the paper is as follows. In the next section (Section 2) we state precisely the problem to be solved and define the HDG methods. We also collect some results from [6] used in our paper. In Section 3 we derive stability (a Gårding-type equality), consistency and duality estimates for the Helmholtz problem. The error estimates for the HDG scheme are obtained in Section 4. In Section 5 we present our numerical results, and finally we draw some conclusions.
2
The hybridizable discontinuous Galerkin method
We start by describing in detail the model problem we shall analyze: the interior Dirichlet problem for the Helmholtz equation. We denote by Ω ⊂ Rn , n = 2, 3, a Lipschitz polyhedral domain with boundary ∂Ω and
2
seek the complex valued solution u of the following problem ∆u + κ2 u = f˜ u=g
in Ω,
(1a)
on ∂Ω,
(1b)
with wave number κ > 0, source term f˜ ∈ L2 (Ω), and boundary values g ∈ H 1/2 (∂Ω). The HDG scheme is based on a first order formulation of the above problems so we introduce q := (i /κ)∇u and f := (i f˜)/κ. Then the Helmholtz equation (1) can be rewritten in mixed form as finding (u, q) such that i κq + ∇u = 0
in Ω,
(2a)
i κu + div q = f
in Ω,
(2b)
on ∂Ω.
(2c)
u=g
Assuming that κ2 is not a Dirichlet eigenvalue for Ω, existence and uniqueness of solutions u ∈ H 1 (Ω) to (1) (or equivalently of solutions (q, u) ∈ H (div, Ω) × H 1 (Ω) to (2)) is well known and these solutions satisfy standard elliptic regularity results, for example kqkH (div,Ω) + kukH 1 (Ω) ≤ C(kf kL2 (Ω) + kgkH 1/2 (∂Ω) )
(3)
where C depends on κ and Ω. Here and throughout C denotes a generic constant, the value of which might be different at different occurrences. The development of the HDG scheme starts by meshing the domain. We consider a subdivision of Ω into a finite element mesh of regular tetrahedra K with faces F in R3 (or triangles K and edges F in R2 ) and denote the resulting collection of triangles/tetrahedra by Th , the collection of edges/faces by Eh , and the collection of element boundaries by ∂Th := {∂K | K ∈ Th }. Regularity here means that the elements are non-degenerate [2]. On each element K and each edge/face F we consider the local spaces of polynomials of degree smaller or equal to k ≥ 0, V (K) := P k (K) := (Pk (K))n ,
W (K) := Pk (K),
and M (F ) := Pk (F ),
where Pk (S) denotes the space of polynomials of total degree at most k in m variables on S where S ⊂ Rm . The corresponding global finite element spaces are given by V h := v ∈ L2 (Ω) v|K ∈ V (K) for all K ∈ Th , Wh := w ∈ L2 (Ω) w|K ∈ W (K) for all K ∈ Th , Mh := µ ∈ L2 (Eh ) µ|F ∈ M (F ) for all F ∈ Eh }, Q where L2 (Ω) := (L2 (Ω))n and L2 (Eh ) := F ∈Eh L2 (F ). Similarly, we will write H s (Ω) := (H s (Ω))n for any s > 0. On these spaces we use the bilinear forms X X X (v, w)Th := (v, w)K , (v, w)Th := (v, w)K , and hv, wi∂Th := hv, wi∂K , K∈Th
R
K∈Th
K∈Th
R
R
with (v, w)K := K v · w dx, (v, w)K := K vw dx, and hv, wi∂K := ∂K vw ds. The hybridizable discontinuous Galerkin method considered in this work belongs to the class of hybridizable local discontinuous Galerkin (LDG-H) methods (cf., e.g., [4]) and yields approximations uh ∈ Wh of u, q h ∈ V h of q, and a numerical trace u bh ∈ Mh approximating u on Eh , which satisfy uh , r · ni∂Th = 0, (i κq h , r)Th − (uh , div r)Th + hb
(4a)
q h · n, wi∂Th = (f, w)Th , (i κuh , w)Th − (q h , ∇w)Th + hb
(4b)
hb uh , µi∂Ω = hg, µi∂Ω , hb q h · n, µi∂Th \∂Ω = 0,
(4c) (4d)
3
for all r ∈ V h , w ∈ Wh , and µ ∈ Mh where the overbar denotes complex conjugation. Here the numerical bh is given by flux q bh = qh + τ (uh − u q bh )n on ∂Th (4e) for some non-negative stabilization function τ : ∂Th → R that is assumed to be constant on each edge/face F ∈ Eh , and n is the outward normal to K ∈ Th . From (2a) we observe that u and q are dimensionally equivalent, and we will show in the next section that this definition yields a consistent scheme. The error analysis for this method presented in the following two sections consists of a variation and extension of arguments used in [6] to analyze a corresponding LDG-H method for a second-order uniformly elliptic problem. It is based on a projection operator Πh inspired by the particular form of the numerical bh from (4e). This projection is defined by flux q Πh : L2 (Ω) × L2 (Ω) ⊃ dom(Πh ) → V h × Wh ,
Πh (q, u) := (ΠV q, ΠW u),
where for any K ∈ Th and all edges/faces F of K the functions ΠV q and ΠW u satisfy (ΠV q, v)K = (q, v)K
for all v ∈ P k−1 (K),
(5a)
(ΠW u, w)K = (u, w)K
for all w ∈ Pk−1 (K),
(5b)
for all µ ∈ Pk (F ).
(5c)
hΠV q · n + τ ΠW u, µiF = hq · n + τ u, µiF
The domain dom(Πh ) of Πh is such that the right hand sides of (5) are well defined, and we note that by (3) solutions to (2) belong to dom(Πh ). max := max τ |∂K > 0, it has been shown in Using the assumed regularity of the mesh, and denoting by τK [6, Thm. 2.1] that the system (5) is uniquely solvable and that its solution (ΠV q, ΠW u) has the following approximation properties: There exists a constant C independent of K, κ and τ such that l +1
kΠV q − qkL2 (K) ≤ ChKq
∗ |q|H lq +1 (K) + ChlKu +1 τK |u|H lu +1 (K) ,
(6a)
l +1
kΠW u − ukL2 (K) ≤ ChlKu +1 |u|H lu +1 (K) + C
hKq max | div q|H lq (K) , τK
(6b)
∗ := max τ |∂K\F ∗ with F ∗ being an edge/face of K at which τ |∂K is maximum. for lu , lq ∈ [0, k], where τK We will also use the following weak commutativity property which is proved in Proposition 2.1 of [6]. For any w ∈ Wh and any (Φ, Ψ) ∈ dom(Πh ),
(w, div Φ)K = (w, div ΠV Φ)K + hw, τ (ΠW Ψ − Ψ)i∂K .
3
(7)
Consistency, Gårding-type identity, and duality
This section is devoted to analyzing the key properties of the scheme that allow us to prove the error estimates in the next section. The arguments in this section rest heavily on those in [6] suitably modified for the non-coercive Helmholtz equation and complex valued solution functions. We denote by PM the L2 -orthogonal projection onto Mh and note that since τ is piecewise constant,
τ (PM u − u), µ ∂T = 0 for all µ ∈ Mh . (8) h
Lemma 3.1 (Consistency). Introducing the projected errors εqh := ΠV q − q h , εuh := ΠW u − uh , and u εb bh , we have h := PM u − u u (i κεqh , r)Th − (εuh , div r)Th + hεb h , r · ni∂Th = (i κ(ΠV q − q), r)Th ,
(i κεuh , w)Th
−
(εqh , ∇w)Th
+ hb εh · n, wi∂Th = (i κ(ΠW u − u), w)Th . u hεb h , µi∂Ω = 0, hb εh · n, µi∂Th \∂Ω = 0, 4
(9a) (9b) (9c) (9d)
for all r ∈ V h , w ∈ Wh , and µ ∈ Mh . Here, u b bh · n εh · n := εqh · n + τ (εuh − εb h ) = PM (q · n) − q
on ∂Th \ ∂Ω.
(9e)
Proof. Recalling that the exact solution (q, u) of (2) satisfies (i κq, r)Th − (u, div r)Th + hu, r · ni∂Th = 0, (i κu, w)Th − (q, ∇w)Th + hq · n, wi∂Th = (f, w)Th for all r ∈ V h and w ∈ Wh , the definitions of Πh and PM yield that (i κΠV q, r)Th − (ΠW u, div r)Th + hPM u, r · ni∂Th = (i κ(ΠV q − q), r)Th (i κΠW u, w)Th − (ΠV q, ∇w)Th + hΠV q · n + τ (ΠW u − u), wi∂Th = (f, w)Th + (i κ(ΠW u − u), w)Th for all r ∈ V h and w ∈ Wh . Now subtracting (4a) from the first equation we obtain (i κ(ΠV q − q h ), r)Th − (ΠW u − uh , div r)Th + hPM u − u ˆh , r · ni∂Th = (i κ(ΠV q − q), r)Th , which yields (9a). Similarly, subtracting (4b) from the second equation above shows that bh ) · n + τ (ΠW u − u), wi∂Th = (i κ(ΠW u − u), w)Th . (i κ(ΠW u − uh ), w)Th − (ΠV q − q h , ∇w)Th + h(ΠV q − q bh in (4e), we find that Recalling the definition of b εh in (9e) and of the numerical flux q b εh · n = (ΠV q − q h ) · n + τ (ΠW u − uh − PM u + u bh ) bh ) · n + τ (uh − u = (ΠV q − q bh ) + τ (ΠW u − u) − τ (uh + PM u − u bh − u), which together with (8) yields (9b). To show (9c) we observe that (4c), (8), and (2c) imply u hεb h , µi∂Ω = hPM u − g, µi∂Ω = hu − g, µi∂Ω = 0.
Finally, using the definition of b εh · n in (9e), (5c), (8), and (4e) we obtain hb εh · n, µi∂Th \∂Ω = h(ΠV q − q h ) · n + τ (ΠW u − uh − PM u + u bh ), µi∂Th \∂Ω = h(q − q h ) · n + τ (u − uh − u + u bh ), µi∂Th \∂Ω q h · n, µi∂Th \∂Ω . = hq · n, µi∂Th \∂Ω − hb bh · n ∈ Mh . Furthermore, since q ∈ H (div, Ω) and q bh satisfies This proves (9e) because b εh · n ∈ Mh and q (4d), both terms on the right hand side vanish, which implies (9d). The following Gårding-type identity is a little more complex than the coercivity bound for the Poisson problem. In addition compared to conforming methods we see the presence of coupling terms on the boundary. From this estimate we see that choosing τ to be real valued gives well controlled terms on element boundaries. Lemma 3.2 (Gårding-type identity). The projected errors satisfy
q u u u b u i κ (εqh , εqh )Th − (εuh , εuh )Th + τ (εuh − εb h ), εh − εh ∂T = i κ(ΠV q − q), εh T + i κ(ΠW u − u), εh T . h
εqh
h
h
εuh
Proof. If we choose r = in (9a), w = in the complex conjugate of (9b), µ = −b εh · n in (9c), and b u µ = −εh in the complex conjugate of (9d) and add these equations we obtain that q u q u u u (i κεqh , εqh )Th − (εuh , div εqh )Th + hεb εh · n, εuh i∂Th h , εh · ni∂Th + (i κεh , εh )Th − (εh , ∇εh )Th + hb q u u u − hεb εh · ni∂Ω − hb εh · n, εb h, b h i∂Th \∂Ω = (i κ(ΠV q − q), εh )Th + (i κ(ΠW u − u), εh )Th .
5
Therefore, (i κεqh , εqh )Th + (i κεuh , εuh )Th + Rh = (i κ(ΠV q − q), εqh )Th + (i κ(ΠW u − u), εuh )Th , with q u q u u εh · n, εuh i∂Th − hεb εh · ni∂Th . Rh = −(εuh , div εqh )Th + hεb h , εh · ni∂Th − (εh , ∇εh )Th + hb h, b
Integrating by parts and applying (9e) yields q u q u u Rh = (∇εuh , εqh )Th − hεuh , εqh · ni∂Th + hεb εh · n, εuh − εb h , εh · ni∂Th − (εh , ∇εh )Th + hb h i∂Th q u u i + hb ε · n, εu − εb i = −hε · n, εu − εb h
h
h ∂Th
h
h
h ∂Th
u u ), εu − εb = hτ (εuh − εb h h i∂Th . h
To obtain error estimates for q h and uh in Section 4 we use a duality argument. To this end, given Θ ∈ L2 (Ω), we introduce the dual problem − i κΦ + ∇Ψ = 0
in Ω,
(10a)
− i κΨ + div Φ = Θ
in Ω,
(10b)
on ∂Ω
(10c)
Ψ=0
and note that (10) is uniquely solvable since we assumed that κ is not a Dirichlet eigenvalue for Ω. Lemma 3.3 (Duality). Given Θ ∈ L2 (Ω) we denote by (Φ, Ψ) the solution to (10). Then, for any Ψh ∈ Wh and any ψh ∈ Wh such that ψh |K ∈ Pk−1 (K) for all K ∈ Th , (εuh , Θ)Th = i κ(q − q h ), ΠV Φ − Φ T + ΠV q − q, ∇Ψ − ∇Ψh T h (11) h + i κ(ΠW u − u), ΠW Ψ − ψh T − i κεuh , ΠW Ψ − Ψ T . h
h
Proof. Using (10b), and the weak commutativity property (7), we obtain (εuh , Θ)Th = (εuh , − i κΨ + div Φ)Th = (εuh , − i κΨ)Th + (εuh , div ΠV Φ)Th + hεuh , τ (ΠW Ψ − Ψ)i∂Th Hence using and (9a) this equality can be rewritten as u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κεqh , ΠV Φ)Th + hεb h , ΠV Φ · ni∂Th − (i κ(ΠV q − q), ΠV Φ)Th
+ hεuh , τ (ΠW Ψ − Ψ)i∂Th u u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hεb h , (ΠV Φ − Φ) · ni∂Th + hεh , τ (ΠW Ψ − Ψ)i∂Th , u where we used in the last step that Φ · n is continuous across all edges/faces and that εb h = 0 on ∂Ω by (9d). Then, using (5c) and (8) this becomes u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hεuh − εb h , τ (ΠW Ψ − Ψ)i∂Th u u b u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th − hτ (εh − εh ), PM Ψi∂Th .
Using the discrete error equations (9e), and (9d), the boundary condition (10c), and the properties of PM in (8), we obtain q u (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th + hεh · n, PM Ψi∂Th q u = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hτ (εuh − εb h ), ΠW Ψi∂Th + hεh · n, Ψi∂Th .
6
Applying (9e), (9b), and integration by parts, we find that εh · n − εqh · n, ΠW Ψi∂Th + hεqh · n, Ψi∂Th (εuh , Θ)Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + hb = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (i κεuh , ΠW Ψ)Th + (εqh , ∇ΠW Ψ)Th − hεqh · n, ΠW Ψi∂Th + hεqh · n, Ψi∂Th = (εuh , − i κΨ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (i κεuh , ΠW Ψ)Th − (div εqh , ΠW Ψ)Th + hεqh · n, Ψi∂Th . Next using (5b), and integration by parts once more, we obtain (εuh , Θ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th − (div εqh , Ψ)Th + hεqh · n, Ψi∂Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (εqh , ∇Ψ)Th . Finally, we use (10a), (5a), and (5b) to obtain (εuh , Θ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th + (i κ(q − q h ), Φ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (ΠV q − q h , ∇Ψ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th − (q − q h , ∇Ψ)Th + (i κ(ΠW u − u), ΠW Ψ)Th + (ΠV q − q h , ∇Ψ)Th = −(i κεuh , ΠW Ψ − Ψ)Th + (i κ(q − q h ), ΠV Φ − Φ)Th + (ΠV q − q, ∇Ψ − ∇Ψh )Th + (i κ(ΠW u − u), ΠW Ψ − ψh )Th which is the desired identity.
4
The error estimates
Before we state the first error estimates in Theorem 4.1 below,Qwe introduce the mesh parameter h := maxK∈Th hK and the mesh dependent norm k · kh on L2 (∂Th ) := K∈Th L2 (∂K) by X kµk2h := hK kµk2∂K for any µ ∈ L2 (∂Th ). K∈Th
Theorem 4.1 (Flux error estimates). Assume k ≥ 1. Given f ∈ L2 (Ω) and g ∈ H 1/2 (∂Ω) we denote by q and u the exact solution of (2), and by q h , uh , and u bh the solution of (4). Furthermore, for any Θ ∈ L2 (Ω) let (Φ, Ψ) be the solution of the dual problem (10) and suppose that kΦkH 1 (Ω) + kΨkH 2 (Ω) ≤ Cκ kΘkL2 (Ω) , where the constant Cκ is independent of Θ but depends on the wave number κ. If 1 ∗ Cκ κhMτmin := Cκ κh max 1, and Cκ κhMτmax := Cκ κh max(1, max hK τK ) max K∈Th minK∈Th τK
(12)
(13)
are sufficiently small, then kΠV q − q h kL2 (Ω) ≤ C (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω)
(14)
and bh · nkh kPM (q · n) − q
κ max 1/2 ≤ C max 1, max (hK τK ) (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) , (15) c K∈Th
where c and C are independent of κ and h for hκ small enough. 7
Note that the regularity assumption (12) on the solution of the dual problem used in Theorem 4.1 is, e.g., well known to hold if Ω is convex (see Brenner and Scott [1, p 139]). Proof. Substituting Θ = εuh into (10) we find from (11) that for any Ψh ∈ Wh and any ψh ∈ Wh such that ψh |K ∈ Pk−1 (K) for all K ∈ Th , (εuh , εuh )Th = (i κ(q − q h ), ΠV Φ − Φ)Th + (ΠV q − q, ∇Ψ − ∇Ψh )Th + (i κ(ΠW u − u), ΠW Ψ − ψh )Th − (i κεuh , ΠW Ψ − Ψ)Th . So, using the Cauchy-Schwartz inequality we obtain that kεuh k2L2 (Ω) ≤ C κkq − q h kL2 (Ω) kΠV Φ − ΦkL2 (Ω) + kΠV q − qkL2 (Ω) k∇Ψ − ∇Ψh kL2 (Ω) + κkΠW u − ukL2 (Ω) kΠW Ψ − ψh kL2 (Ω) + κkεuh kL2 (Ω) kΠW Ψ − ΨkL2 (Ω) .
(16)
The first term on the right hand side of (16) can be estimated by means of (6a) with lu = 1 and lq = 0, ∗ κkq − q h kL2 (Ω) kΠV Φ − ΦkL2 (Ω) ≤ Cκhkq − q h kL2 (Ω) |Φ|H 1 (Ω) + max hK τK |Ψ|H 2 (Ω) , K∈Th
Similarly, using the projection error bound |Ψ − Ψh |H 1 (Ω) ≤ Ch|Ψ|H 2 (Ω) we find for the second term that kΠV q − qkL2 (Ω) k∇Ψ − ∇Ψh kL2 (Ω) ≤ ChkΠV q − qkL2 (Ω) |Ψ|H 2 (Ω) . To estimate the third term in (16), we apply the triangle inequality, (6b) with lu = 0 and lq = 0, and the projection error bound kΨ − ψh kL2 (Ω) ≤ Ch|Ψ|H 1 (Ω) , which yields κkΠW u − ukL2 (Ω) kΠW Ψ − ψh kL2 (Ω) ≤ κkΠW u − ukL2 (Ω) kΠW Ψ − ΨkL2 (Ω) + kΨ − ψh kL2 (Ω) 1 | div Φ|L2 (Ω) . ≤ CκhkΠW u − ukL2 (Ω) |Ψ|H 1 (Ω) + max minK∈Th τK Similarly, we find for the last term on the right hand side of (16) that 1 2 (Ω) . | div Φ| κkεuh kL2 (Ω) kΠW Ψ − ΨkL2 (Ω) ≤ Cκhkεuh kL2 (Ω) |Ψ|H 1 (Ω) + L max minK∈Th τK Combining these estimates using (12) with Θ = εuh we obtain that (16) reduces to kεuh k2L2 (Ω) ≤ CCκ κhMτmax kq − q h kL2 (Ω) kεuh kL2 (Ω) + CCκ hkΠV q − qkL2 (Ω) kεuh kL2 (Ω) + CCκ κhMτmin kΠW u − ukL2 (Ω) kεuh kL2 (Ω) + CCκ κhMτmin kεuh k2L2 (Ω) . Therefore, assuming that Cκ κhMτmin is small enough, 1 kεuh kL2 (Ω) ≤ CCκ κh Mτmax kεqh kL2 (Ω) + + Mτmax kΠV q − qkL2 (Ω) + Mτmin kΠW u − ukL2 (Ω) . κ
(17)
Now recalling Lemma 3.2, using the arithmetic-geometric mean inequality with some α > 0, we find from (17) that i q q q u u u u u b u (εh ,εh )Th − hτ (εuh − εb h ), εh − εh i∂Th = (ΠV q − q), εh )Th − ((ΠW u − u), εh )Th + (εh , εh )Th κ α 1 1 3 ≤ kεqh k2L2 (Ω) + kΠV q − qk2L2 (Ω) + kΠW u − uk2L2 (Ω) + kεuh k2L2 (Ω) 2 2α 2 2 α 1 1 q ≤ + CCκ2 κ2 h2 (Mτmax )2 kεh k2L2 (Ω) + + CCκ2 κ2 h2 2 + (Mτmax )2 kΠV q − qk2L2 (Ω) 2 2α κ 1 + + CCκ2 κ2 h2 (Mτmin )2 kΠW u − uk2L2 (Ω) . 2 8
Therefore, choosing α small enough, assuming that Cκ κhMτmax is small, and recalling that Cκ κhMτmin is small, we obtain that c u 2 2 2 2 u b u (18) kεqh k2L2 (Ω) + hτ (εuh − εb h ), εh − εh i∂Th ≤ C (1 + Cκ h )kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) , κ for some constants c > 0 and C > 0, which implies (14). To show (15), we recall (9e) and use an inverse inequality (element by element) to obtain u bh · nkh ≤ kεqh · nkh + kτ (εuh − εb kPM (q · n) − q h )kh max 1/2 u u b u 1/2 ≤ Ckεqh kL2 (Ω) + max (hK τK ) hτ (εuh − εb h ), εh − εh i∂Th K∈Th κ c u max 1/2 u − εb u i1/2 . ≤ C max 1, max (hK τK kεqh kL2 (Ω) + hτ (εuh − εb ), ε ) h h h ∂Th c K∈Th κ
Therewith, (15) follows from (18). Theorem 4.2 (Convergence of the scalar variable). Suppose that the assumptions of Theorem 4.1 hold. If Cκ κhMτmin
and
Cκ κhMτmax
(as defined in (13)) are sufficiently small, then kΠW u − uh kL2 (Ω) ≤ CCκ κh
1 κ
+ Mτmax kΠV q − qkL2 (Ω) + (Mτmin + Mτmax )kΠW u − ukL2 (Ω)
(19)
and kPM u − u bh kh ≤ Cκh
1 + Cκ h + Cκ
1 κ
+ Mτmax
kΠV q − qkL2 (Ω)
+ 1+
Cκ (Mτmin
+
Mτmax )
kΠW u − uk
L2 (Ω)
. (20)
Remark 4.3. Collecting together the results of Theorems 4.1 and 4.2 together with the projection estimates (6a) and (6b) we can obtain explicit bounds for the error in terms of the mesh size. A practical case is when τ = O(1) and u ∈ H k+2 (Ω) (convergence for less regular solutions follows in the obvious way). Then we obtain ku − uh kL2 (Ω) + kq − q h kL2 (Ω) = O(hk+1 ) . Similarly, Theorem 4.2 guarantees super-convergence of the projected error of the scalar variable and of the numerical trace, which means for τ = O(1) and u ∈ H k+2 (Ω) that kΠW u − uh kL2 (Ω) = O(hk+2 )
and
kPM u − u bh kh = O(hk+2 ) ,
respectively. Proof. Combining (17) and (14), and recalling that Cκ κhMτmax is small, we find that kΠW u − uh kL2 (Ω) ≤ CCκ κh Mτmax (1 + Cκ h)kΠV q − qkL2 (Ω) + kΠW u − ukL2 (Ω) 1 + + Mτmax kΠV q − qkL2 (Ω) + Mτmin kΠW u − ukL2 (Ω) κ 1 ≤ CCκ κh (1 + κMτmax )kΠV q − qkL2 (Ω) + (Mτmin + Mτmax )kΠW u − ukL2 (Ω) , κ which is (19). 9
u For (20) we follow the proof of [6, Thm. 4.1] and select r ∈ P k (K) such that r · n = εb h on ∂K and 1/2 b u krkL2 (K) ≤ ChK kεh kL2 (∂K) . Now substituting hK r into (9a) and applying an inverse inequality (see [1, Lmm. 4.5.3]) we obtain q u u 2 hK kεb h kL2 (∂K) = hK (i κ(ΠV q − q), r)K + hK (εh , div r)K − hK (i κεh , r)K ≤ CκhK krkL2 (K) kΠV q − qkL2 (K) + kεqh kL2 (K) + CkrkL2 (K) kεuh kL2 (K) 1/2 u q u ≤ ChK kεb h kL2 (∂K) κhK kΠV q − qkL2 (K) + κhK kεh kL2 (K) + kεh kL2 (K) .
Therefore, using (19) and (14) we find that kPM u − u bh k2h ≤ C κ2 h2 kΠV q − qk2L2 (Ω) + κ2 h2 kεqh k2L2 (Ω) + kεuh k2L2 (Ω) ≤ C κ2 h2 kΠV q − qk2L2 (Ω) + κ2 h2 (1 + Cκ h)2 kΠV q − qk2L2 (Ω) + kΠW u − uk2L2 (Ω) +
Cκ2 κ2 h2
≤ Cκ2 h2
1 κ
+
Mτmax
1 + Cκ2 h2 + Cκ2
2
1
kΠV q −
+ Mτmax
qk2L2 (Ω)
+
(Mτmin
+
Mτmax )2 kΠW u
−
uk2L2 (Ω)
2
kΠV q − qk2L2 (Ω) 2 min max 2 2 + 1 + Cκ (Mτ + Mτ ) kΠW u − ukL2 (Ω) , κ
which yields (20). Remark 4.4. The super-convergence of the approximation uh of the scalar variable u to the projection ΠW u established in Theorem 4.2 can be exploited to obtain an approximation u∗h of improved accuracy. As outlined in section 5 of [6], the postprocessing schemes for mixed finite elements proposed by Stenberg [11] can be suitably adapted and applied to the hybridizable DG method (see also [3] and Cockburn, Guzmán, and Wang [5]). Following the analysis in [11] it can be shown that (for k ≥ 1) the postprocessed solution u∗h converges of one order higher than uh .
5
Numerical results
In this section we probe the error estimates established in Theorem 4.1 and Theorem 4.2 by numerical examples for convex and non-convex two-dimensional domains, and for solutions of differing regularity. The HDG-method has been implemented for piecewise constant, piecewise linear, and piecewise quadratic finite element spaces in R2 . We use order 6 quadrature rules on the triangles as well as on the edges to evaluate the system matrix and the right hand side corresponding to (4). The mesh generation and all computations are done in MATLAB. Example 5.1. In the first example we consider a square domain Ω =]−1, 1[×]−1, 1[⊂ R2 and structured meshes Th with N = 8, 32, 128, 512, and 2048 triangles as shown in Figure 1 (left) for N = 128. We choose the boundary conditions g and the source term f in (4) such that the exact solution of (1) is given by p u(x) = Jξ κ (x + 1)2 + y 2 cos ξ arctan(y/(x + 1)) for ξ = 1, ξ = 2/3, and ξ = 3/2, respectively, where Jξ denotes the Bessel function of the first kind and order ξ. It is well known (see, e.g., Hiptmair, Moiola, and Perugia [8] and Grisvard [7] for details) that u(x) is smooth for ξ ∈ N, while its derivative has a singularity at (−1, 0) for ξ ∈ / N. More precisely, u ∈ H 2 (Ω) 3 1 2 but u ∈ / H (Ω) for ξ = 3/2 and u ∈ H (Ω) but u ∈ / H (Ω) for ξ = 2/3. 10
1
1
0.5
0.5
0
0
ï0.5
ï0.5
ï1 ï1
ï0.5
0
0.5
ï1 ï1
1
ï0.5
0
0.5
1
Figure 1: Visualization of structured meshes used in Example 5.1 (left) and Example 5.2 (right). Figure 2 shows the absolute error in the scalar variable ku − uh kL2 (Ω) , the absolute error in the flux variable kq − q h kL2 (Ω) , and the absolute error in the numerical trace ku − u bh kh as a function of κh/2π for ξ = 1, ξ = 2/3, and ξ = 3/2 (top down) and wave number κ = 4, i.e., the wave length is λ = 2π/κ ≈ 1.57. The plots in the left column correspond to k = 1, i.e., piecewise linear approximations, while the right column has been obtained for k = 2, i.e., piecewise quadratic approximations. As predicted by (14), (19), and (20) together with (6) we observe optimal convergence in qh , uh , and u bh for the smooth solution, i.e., ξ = 1, using linear as well as quadratic elements. The singular solutions with ξ = 2/3 and ξ = 3/2 show a slightly different behavior: For linear elements and ξ = 2/3, i.e., u ∈ H 1 (Ω) but u ∈ / H 2 (Ω), we loose one order of convergence in the flux variable, while the scalar variables still converge at optimal order. On the other hand, for quadratic elements and ξ = 2/3 the convergence in all three variables kq − q h kL2 (Ω) , ku − uh kL2 (Ω) , and ku − u bh kh slows down. If ξ = 3/2, i.e., u ∈ H 2 (Ω) but u ∈ / H 3 (Ω), the convergence is optimal in all three variables for the piecewise linear approximation while the convergence in the flux variable slows down for the piecewise quadratic approximation. This behavior is actually better than our theoretical results predict: By (6) our estimates guarantee optimal convergence of ku − uh kL2 (Ω) only if the exact solution u is in H k+1 (Ω) and optimal convergence order in kq − q h kL2 (Ω) only if u ∈ H k+2 (Ω). For the piecewise constant approximation, i.e., k = 0, we observed (results not shown) optimal convergence for all three choices of ξ and sufficiently small h after a comparatively large pre-asymptotic area of slower convergence. Next we study the dependence of the convergence rate on the wave number κ. To this end we repeat the previous experiment for ξ = 1, quadratic elements, and four different wave numbers κ = 1, κ = 4, κ = 10, and κ = 16. Figure 3 shows ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh as a function of κh/2π. In all four examples we observe optimal convergence order for sufficiently small values of h. With increasing wave number the absolute error increases and a pre-asymptotic region of slower convergence for large values of h appears, as is to be expected because of pollution error in the h version. If we order the degrees of freedom counting first those for uh and q h and then those for u ˆh we obtain a block linear system A B x b = y c BT C where x denotes the vector of degrees of freedom for uh and q h while y contains the degrees of freedom for the trace u ˆh . Because of the use of a discontinuous Galerkin method, the matrix A is block diagonal (assuming degrees of freedom are ordered element by element). Thus assuming it is invertible (numerically the case in our study) we can cheaply eliminate x and obtain the Schur complement system Sy = d where S = (C − B T A−1 B) and d = c − B T A−1 b. Of course in practical calculations the matrix S would be computed element by element so that it would never be necessary to assemble A or B. The Schur 11
1
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
error
error
10 ï1
10
ï2
10
ï3
10 ï2
10
ï4
10 ï3
10
ï5
0.05
0.1
10
0.5
0.05
0.1
gh/2/ 1
1
10
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3) O(h)
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
0
10
ï1
error
10 error
0.5 gh/2/
ï1
10
ï2
10
ï3
10 ï2
10
ï4
10 ï3
10
ï5
0.05
0.1
10
0.5
0.05
0.1
gh/2/
0.5 gh/2/
1
1
10
10
2
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h2)
0
10
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
10 ï1
error
error
10
ï2
10
ï2
10
ï3
10 ï3
10
ï4
10
ï4
10
ï5
0.05
0.1
10
0.5
0.05
gh/2/
0.1
0.5 gh/2/
Figure 2: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, ξ = 2/3, and ξ = 3/2 (top down) and wave number κ = 4 for linear elements (left) and quadratic elements (right).
12
1
2
10
10
0
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
ï1
10 ï2
error
error
10
ï2
10
ï4
10
ï3
10 ï6
10
ï4
10
ï8
ï5
10
0.05 gh/2/
10
0.1
2
3
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
2
10
0
error
error
10
ï1
10
ï2
0
10
ï1
10
10
ï3
ï2
10
10
ï4
0.1
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
1
10
10
0.5
10 2
10
0.1 gh/2/
10
1
0.05
ï3
10
0.5 gh/2/
0.5 gh/2/
Figure 3: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1 and κ = 1, κ = 4, κ = 10, and κ = 16 (from left to right and top down) for quadratic elements.
13
κ=1 κ=4 κ = 10 κ = 16
N =8 2.98e+02 9.76e+01 1.11e+02 6.14e+01
N = 32 9.84e+02 3.98e+02 2.39e+02 2.95e+02
N = 128 2.91e+03 1.44e+03 3.64e+03 6.70e+02
N = 512 9.55e+03 5.60e+03 1.75e+04 1.67e+04
N = 2048 3.40e+04 2.31e+04 7.82e+04 1.37e+05
Table 1: Condition numbers of the Schur complement matrix S with respect to the wave number κ and the number of triangles N for quadratic elements. 1
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
10
10
ï1
ï1
10 error
error
10
ï2
10
ï3
ï2
10
ï3
10
10
ï4
ï4
10
10
ï5
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
0
ï5
0.05
0.1
10
0.5 gh/2/
0.05
0.1
0.5 gh/2/
Figure 4: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, wave number κ = 4, and τ = h, τ = 1, and τ = 1/h (from left to right) for quadratic elements.
complement matrix S (equivalently, the matrix for the condensed hybridized system) is sparse but the system would typically be solved iteratively in practice so it is of interest to examine the condition number of S as a function of the number of triangles N and wave number κ as shown in Table 1 for quadradic elements. In the examples considered so far the coupling parameter τ from (4e) has been chosen to be τ = 1. To analyze the dependence of the numerical solution on this parameter, we show in Figure 4 the convergence behavior for the smooth solution with ξ = 1 and wave number κ = 4 for τ = h and τ = 1/h. These plots should be compared with the corresponding plot for τ = 1 in Figure 2 (top right). While we observed optimal convergence for τ = 1, the convergence of ku − uh kL2 (Ω) slows down for τ = h, while for τ = 1/h the convergence speed of the flux variable kq − q h kL2 (Ω) is reduced. This reflects exactly the approximation properties of the projection operator Πh stated in (6) that have been used in the convergence analysis. The behavior for τ = 1/h clearly follows from (14), (19), and (20), but this is not so obvious for τ = h. To exemplify the super-convergence properties of the HDG-scheme established in Theorem 4.2, we show in Figure 5 plots of kPM u − u bh kh corresponding to the smooth solution with ξ = 1, wave number κ = 4, and τ = 1 for linear elements (left) and quadratic elements (right). Super-convergence is clearly observed: As predicted in (20) the convergence order of kPM u − u bh kh is O(h3 ) for linear elements and O(h4 ) for quadratic elements. Example 5.2. In the second example we consider a non-convex domain as shown in Figure 1 (right). The boundary conditions g and the source term f in (4) are chosen similarly to the first example such that the exact solution of (1) is given by p u(x) = Jξ κ x2 + y 2 cos ξ arctan(y/x))
14
1
0
10
10 L2ïerror of u L2ïerror of q projected hnormïerror of uhat O(h2) O(h3)
0
10
L2ïerror of u L2ïerror of q projected hnormïerror of uhat O(h3) O(h4)
ï1
10
ï2
10
ï1
error
error
10
ï3
10
ï2
10
ï4
10 ï3
10
ï5
10
ï4
10
ï6
0.05
0.1
10
0.5 gh/2/
0.05
0.1
0.5 gh/2/
Figure 5: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and kPM u − u bh kh with ξ = 1, wave number κ = 4, and τ = 1 for linear elements (left) and quadratic elements (right).
for ξ = 1, ξ = 2/3, and ξ = 3/2. The smoothness properties of these solutions are exactly the same as in the first example, but for ξ = 2/3, and ξ = 3/2 the singularity in the derivative of the exact solution is now at (0, 0). Due to the non-convexity of this domain, Theorem 4.1 and Theorem 4.2, which have both been obtained using duality arguments, do not even apply for the smooth solution. Nevertheless, we compute piecewise quadratic numerical approximations of the solution to (1) with wave number κ = 4 and coupling parameter τ = 1 for regular structured meshes with N = 24, 96, 384, and 1536 elements (see Figure 1 (right) for a visualization with N = 96). Plots of the corresponding numerical errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh for ξ = 1, ξ = 2/3, and ξ = 3/2 can be found in Figure 6. We observe the same qualitative behavior as for the convex domain: For ξ = 1 we have optimal convergence in all three variables, while if ξ = 3/2 the convergence of the flux variable slows, and if ξ = 2/3 the flux variable, the scalar variable and the numerical trace no longer converge at optimal order. Since we assume that the convergence order for ξ = 3/2 and ξ = 2/3 degrades because of to the singularity of the derivative of the exact solution at (0, 0), we consider in our final experiment a sequence of unstructured locally refined meshes with N = 197, 254, 417, and 629 elements as shown in Figure 7. Corresponding graphs of the absolute errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh as a function of the number of triangles N can be found for κ = 4, quadratic elements, and ξ = 2/3 (left) as well as for ξ = 3/2 (right) in Figure 8. Indeed the errors decrease when we refine the grid locally around the singularity. However convergence is slower for the smoother (ξ = 3/2) solution in the right panel compared to that in the left panel (ξ = 2/3). This is because the locally refined grids are chosen for the more singular problem and are under refined away from the singularity when ξ = 3/2.
Concluding Remarks We have show that the HDG scheme of [4] can be applied to the Helmholtz equation. In particular techniques from [6] can be used to prove error estimates, and our numerical results demonstrate that the method can be applied even when the restrictive assumptions of our results do not hold. Indeed the numerical results suggest that it should be possible to prove error estimates for non-convex domains and improve on the regularity requirements of the solution in order to obtain a particular rate. Clearly the main open problem is to prove h − p convergence results, since in practice we need to adjust p and h in order to control dispersion element by element (see Melenk and Sauter [9]). Our numerical results suggest that HDG methods have great promise to provide a flexible and efficient technique for solving linear time harmonic wave propagation problems.
15
0
1
10
10 L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
ï1
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(h3) O(h)
0
10
ï2
10
ï1
error
error
10 ï3
10
ï2
10 ï4
10
ï3
10
ï5
10
ï6
ï4
10
0.05
10
0.1 gh/2/
0.05
0.1 gh/2/
0
10
2
L ïerror of u L2ïerror of q hnormïerror of uhat O(h3)
ï1
10
ï2
error
10
ï3
10
ï4
10
ï5
10
0.05
0.1 gh/2/
Figure 6: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh with ξ = 1, ξ = 2/3, and ξ = 3/2 (from left to right and top down) for a non-convex domain and quadratic elements.
1
1
1
1
0.5
0.5
0.5
0.5
0
0
0
0
ï0.5
ï0.5
ï0.5
ï0.5
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
0
0.5
1
ï1 ï1
ï0.5
Figure 7: Visualization of the locally refined meshes used in Example 5.2.
16
0
0.5
1
ï1
ï1
10
10
ï2
ï2
10
error
error
10
ï3
10
ï4
10
ï5
10
ï4
10
L2ïerror of u L2ïerror of q hnormïerror of uhat O(Nï3/2)
200
300 400 number of triangles N
ï3
10
ï5
500
10
600
200
L2ïerror of u L2ïerror of q hnormïerror of uhat O(Nï3/2) 300 400 number of triangles H
500
600
Figure 8: Errors ku − uh kL2 (Ω) , kq − q h kL2 (Ω) , and ku − u bh kh for the non-convex domain and locally refined grids with quadratic elements and ξ = 2/3 (left) as well as ξ = 3/2 (right).
Acknowledgements We would like to thank anonymous referees for comments resulting in major improvements to this paper. The research of PM was supported in part by grant number FA9550-08-1-0138 from the US AFOSR.
References [1]
S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, 3rd ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008.
[2]
P. G. Ciarlet, The finite element method for elliptic problems, Studies in Mathematics and its Applications, vol. 4, NorthHolland Publishing Co., Amsterdam-New York-Oxford, 1978.
[3]
B Cockburn, B Dong, and J Guzman, A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems, Math. Comp. 77 (2008), 1887-1916.
[4]
B. Cockburn, J. Gopalakrishnan, and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365.
[5]
B. Cockburn, J. Guzmán, and H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comp. 78 (2009), no. 265, 1–24.
[6]
B. Cockburn, J. Gopalakrishnan, and F.-J. Sayas, A projection-based error analysis of HDG methods, Math. Comp. 79 (2010), no. 271, 1351–1367.
[7]
P. Grisvard, Elliptic problems in nonsmooth domains, Monographs and Studies in Mathematics, vol. 24, Pitman, Boston, MA, 1985.
[8]
R. Hiptmair, A. Moiola, and I. Perugia, Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the p-version, SAM-ETH Zürich Report 2009-20 (2009).
[9]
J. M. Melenk and S. Sauter, Convergence analysis for finite element discretizations of the Helmholtz equation with Dirichletto-Neumann boundary conditions, Math. Comp. 79 (2010), 1871–914.
[10] A. H. Schatz, An observation concerning Ritz-Galerkin methods with indefinite bilinear forms, Math. Comp. 28 (1974), 959–962. [11] R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO Modél. Math. Anal. Numér. 25 (1991), no. 1, 151–167.
17
