# A hybridized discontinuous Galerkin method for Poisson-type problems with sign-changing coefficients

In this paper, we present a hybridized discontinuous Galerkin (HDG) method for Poisson-type problems with sign-changing coefficients. We introduce a sign-changing stabilization parameter that results in a stable HDG method independent of domain geometry and the ratio of the negative and positive coefficients. Since the Poisson-type problem with sign-changing coefficients is not elliptic, standard techniques with a duality argument to analyze the HDG method cannot be applied. Hence, we present a novel error analysis exploiting the stabilized saddle-point problem structure of the HDG method. Numerical experiments in two dimensions and for varying polynomial degree verify our theoretical results.

## Authors

• 6 publications
• 18 publications
09/30/2020

### Superconvergence of Discontinuous Galerkin methods for Elliptic Boundary Value Problems

In this paper, we present a unified analysis of the superconvergence pro...
02/25/2020

### A generalized finite element method for problems with sign-changing coefficients

Problems with sign-changing coefficients occur, for instance, in the stu...
05/20/2022

### An optimal control-based numerical method for scalar transmission problems with sign-changing coefficients

In this work, we present a new numerical method for solving the scalar t...
01/29/2021

### Homogeneous multigrid for embedded discontinuous Galerkin methods

We introduce a homogeneous multigrid method in the sense that it uses th...
01/12/2020

### Computer-assisted analysis of the sign-change structure for elliptic problems

In this paper, a method is proposed for rigorously analyzing the sign-ch...
06/21/2021

### A posteriori goal-oriented bounds for the Poisson problem using potential and equilibrated flux reconstructions: application to the hybridizable discontinuous Galerkin method

We present a general framework to compute upper and lower bounds for lin...
12/14/2020

### Simultaneous approximation terms and functional accuracy for diffusion problems discretized with multidimensional summation-by-parts operators

Several types of simultaneous approximation term (SAT) for diffusion pro...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Partial differential equations (PDE) with sign-changing coefficients are used to model physical phenomena of meta-materials, for example, wave transmission problems between classical materials and meta-materials. Therefore, there is an emerging interest in the development of numerical methods for such PDEs. However, the bilinear forms for problems with sign-changing coefficients are not coercive, so standard techniques relying on coercivity cannot be applied. The lack of coercivity indeed poses fundamental challenges to study the well-posedness of such PDEs, as well as the development of numerical methods.

A popular approach in the study of numerical methods for PDEs with sign-changing coefficients is the “-coercivity” approach. This approach utilizes a linear operator that recovers the coercivity of the bilinear form. The existence of such an operator is a natural assumption because the existence of is equivalent to the well-posedness of the PDE in the sense of Banach–Nečas–Babuška in the functional analysis framework [7]. The application of the operator to develop numerical methods was first proposed in [2]. Since then, the -coercivity approach has been used in the development of numerical methods for Poisson-type problems [7], Helmholtz-type problems [8], and Maxwell-type problems [2, 5, 4].

In numerical methods that use -coercivity, the operator cannot be used directly, because the domain and the range of are not the discrete spaces in general. For this reason, is introduced, a discrete approximation of , which has the coercivity recovery property in the discrete setting. The existence of such is non-trivial and difficult to guarantee because the form of is unknown in general. For general geometry and meshes sufficient conditions for the existence of are proposed in [7]

, but these include non-quantitative constants such as a norm of the discrete interpolation of

and the ratio of the positive and negative coefficients (the contrast). It is furthermore shown in [7] that locally symmetric meshes across the transmission interface, where the coefficient changes sign, improve the quality of numerical solutions. However, the generation of such meshes is a non-trivial restriction in mesh generation for complex interfaces.

Hybridization (or static condensation) was originally introduced for mixed finite element methods to reduce computational cost [15], and it was further analyzed in [1]. More recently, hybridization was combined with the discontinuous Galerkin (DG) method. This resulted in the hybridizable discontinuous Galerkin (HDG) method which was shown to be more efficient than standard DG methods. The HDG method was systematically presented in [9] for elliptic partial differential equations. Since its introduction, it has been extended to many different problems. These include HDG methods for the Helmholtz problem, e.g., [12, 17] and the Maxwell problem, e.g., [6, 14].

In this paper, we develop a numerical method for Poisson-type problems with a sign-changing coefficient that avoids the -coercivity approach. The numerical method we present in this paper is always well-posed without any conditions on the domain geometry or the ratio between the negative and positive coefficients. This is achieved by employing two key ideas for the problem: a mixed formulation of the problem and using a hybridized discontinuous Galerkin (HDG) method with sign-changing stabilization parameters. We will see that the sign-changing stabilization parameter and the discontinuous test function space in discontinuous Galerkin methods allow us to obtain stable numerical methods without any non-quantitative assumptions on mesh size or coefficients.

We also carry out an error analysis of the HDG method for this problem. The error analysis of HDG methods typically uses the Aubin–Nitsche duality argument with the full elliptic regularity assumption. However, this assumption is not feasible in our setting, because the PDE is not elliptic. By revealing a stabilized saddle point problem structure of HDG methods, we circumvent this obstacle and we can show an error estimate without using the duality argument. We note that such analysis has been applied also to HDG methods for standard Poisson-type problems to avoid the full elliptic regularity assumption in the a priori error analysis (see

[11]).

The outline of this paper is as follows. In Section 2 we introduce the Poisson problem with a sign-changing coefficient. We introduce the HDG method and discuss its well-posedness in Section 3. An a priori error analysis is presented in Section 4. In Section 5 we show that a superconvergence result can be obtained when a suitable regularity assumption is available. Our analysis is verified by numerical examples in Section 6, while conclusions are drawn in Section 7.

## 2 The Poisson problem with a sign-changing coefficient

Let , , be a polygonal domain that is divided into two subdomains and such that and . The boundaries of , and are denoted by, respectively, , and . The interface separating the domains is denoted by and we define furthermore and . We assume that and are disjoint subsets of such that

. The outward unit normal vector field on

is denoted by .

Let be a scalar function defined as

 σ:={σ+on Ω+,σ−on Ω−, (1)

where and are constants. The contrast is defined as . Throughout this paper we assume that

 0<σmin≤σ+,−σ−≤σmax<+∞. (2)

We consider the following partial differential equation (PDE) for the scalar : equationparentequation

 ∇⋅\delσ∇u =f in Ω, (3a) u =uD on ΓD, (3b) −σ∇u⋅n =uN on ΓN, (3c)

where and are given boundary data and is a given source term. This is not a standard second order PDE, because the sign of is indefinite. Let where is the standard Sobolev space of functions such that the -norm of their gradient is bounded. The variational formulation for the problem (when ) is given by: Find such that

 ∫Ωσ∇u⋅∇v\difx=∫Ωfv\difx∀v∈H10(Ω). (4)

It is known that this problem is not always well-posed, for example, when , the problem is ill-posed. The conditions of well-posedness of the problem depends on the values of and the geometry of and . For more information on the well-posedness of this problem we refer to [7, 3].

## 3 The HDG method

Introducing the auxiliary variable , Eq. 3 can be written as a system of first-order equations: equationparentequation

 σ−1q+∇u =0 in Ω, (5a) ∇⋅q =−f in Ω, (5b) u =uD on ΓD, (5c) −σ∇u⋅n =uN on ΓN. (5d)

In this section we introduce the HDG method for Eq. 5.

### 3.1 Preliminaries

Let be a triangulation of and be the set of faces in the triangulation with co-dimension 1. By we denote the subset of such that its elements are in . is defined similarly. We assume that is a conforming triangulation with respect to , i.e., there is a subset of such that . For later reference we define and as the subsets of interior facets of such that the facets are in and , respectively. We also define as the facets on , i.e., . In summary, is a disjoint union of , , , . We also define as the subset of such that .

Let , denote the diameter of by , and let . For scalar functions and vector functions we denote, respectively, and . Furthermore, for two functions and with well-defined traces on , we define , where is the unit outward normal vector field on . Additionally, we define , , and , with similar definitions on and .

We use the standard notation for Sobolev spaces based on the -norm, i.e., , , denotes the Sobolev space based on the -norm with -differentiability on the domain . We refer to [13] for a rigorous definition of this space. The norm on is denoted by . When we drop the subscript .

To define the HDG method let

 Vh(K)=Pk(K;Rd),Wh(K)=Pk(K),Mh(F)=Pk(F),

where and are the spaces of scalar and -dimensional vector valued polynomials of degree at most on a domain . We will use the following finite element spaces:

 Vh =\cbr[0]r∈L2(Th;Rd):r|K∈Vh(K),∀K∈Th, (6) Wh =\cbr[0]v∈L2(Th):v|K∈Wh(K),∀K∈Th, (7) Mh =\cbr[0]μ∈L2(Fh):μ|F∈Mh(F),∀F∈Fh,μ|FDh=0. (8)

### 3.2 The discrete formulation

For a sufficiently regular solution of the mixed problem Eq. 5, and for the restriction of on , we can derive a system of variational equations from Eq. 5 with test functions in as equationparentequation

 \del[0]σ−1q,rΩ+\del[0]∇u,rΩ−⟨u−¯u,r⋅n⟩∂Th =0 r∈Vh, (9a) \del[0]q,∇vΩ−⟨q⋅n+τ(u−¯u),v⟩∂Th =\del[0]f,vΩ v∈Wh, (9b) ⟨q⋅n+τ(u−¯u),¯v⟩∂Th∖FDh =⟨uN,¯v⟩FNh ¯v∈Mh. (9c)

Here we use to denote a piecewise constant function adapted to ; in the integration on is the trace of . Our HDG method for Eq. 5 is the discrete counterpart of this system of variational equations, i.e., we find such that equationparentequation

 \del[0]σ−1qh,rΩ+\del[0]∇uh,rΩ−⟨uh−¯uh,r⋅n⟩∂Th r∈Vh, (10a) \del[0]qh,∇vΩ−⟨qh⋅n+τ(uh−¯uh),v⟩∂Th =−⟨τuD,v⟩FDh+\del[0]f,vΩ v∈Wh, (10b) ⟨qh⋅n+τ(uh−¯uh),¯v⟩∂Th =⟨uN,¯v⟩FNh ¯v∈Mh. (10c)

Here equations Eq. 9a, Eq. 9b and the discrete equations Eq. 10a, Eq. 10b look inconsistent. This is because we impose the restriction on so that we can take the same trial and test function spaces in our numerical method. This choice will be useful to reveal a stabilized saddle point structure of our numerical method, and also allow us to obtain optimal error estimates without the Aubin–Nitsche duality argument. To clarify the consistency between Eq. 9a, Eq. 9b and Eq. 10a, Eq. 10b, we point out that Eq. 10a, Eq. 10b with are equivalent to the discretized forms of Eq. 9a, Eq. 9b with defined by

 ~Mh=\cbr[0]μ∈L2(Fh):μ|F∈Mh(F) ∀F∈Fh, ⟨μ,λ⟩FDh=⟨uD,λ⟩FDh ∀λ∈Mh(FDh).

For later use we define

 ah(q,r):=\del[0]σ−1q,rΩ,bh(v,¯v;r):=\del[0]∇v,rΩ−⟨v−¯v,r⋅n⟩∂Th,ch(u,¯u;v,¯v):=⟨τ(u−¯u),v−¯v⟩∂Th. (11)

Then the sum of the left-hand sides of Eq. 10 can be rewritten as

 Bh(qh,uh,¯uh;r,v,¯v):=ah(qh,r)+bh(uh,¯uh;r)+bh(v,¯v;qh)−ch(uh,¯uh;v,¯v). (12)

For brevity, but without loss of generality, we assume and in the remainder of this paper.

In the following lemma we prove well-posedness of the HDG method Eq. 10.

###### Lemma 1 (Well-posedness)

Let . If on and on , there exists a unique solution to Eq. 10.

###### Proof

To show well-posedness of Eq. 10 assume that . We will show that , , vanish.

Let

be the characteristic function that has the value 1 in domain

and 0 elsewhere. Then let , , in Eq. 10, and add all the equations. The evaluation of Eq. 10 with the first components of these test functions gives

 \del[0]σ−1qh,qhΩ++⟨τ(uh−¯uh),uh−¯uh⟩∂T+h=0.

Similarly, the evaluation of Eq. 10 with the second components of the above test functions gives

 −\del[0]σ−1qh,qhΩ−−⟨τ(uh−¯uh),uh−¯uh⟩∂T−h=0.

Since and are positive and negative on and , respectively, we can conclude that on and on . This also implies that is continuous on .

Since and it follows from Eq. 10a that , i.e., is constant on . Then on because on .

We consider the following norms on and :

 \norm[0]r2Vh =\del[0]σ−1r,rΩ+−\del[0]σ−1r,rΩ−, \norm[0](v,¯v)2Wh×Mh =\del[0]v,vΩ+⟨|τ|\del[0]v−¯v,v−¯v⟩∂Th,

where is the absolute value of . Additionally, we define with norm

 \norm[0](r,v,¯v)Xh:=\del[1]\norm[0]r2Vh+\norm[0](v,¯v)2Wh×Mh12.

Furthermore, in the remainder of this paper we let be a constant independent of the mesh size .

###### Lemma 2

There exists which depends only on and such that for there exists satisfying and with on .

###### Proof

We will use a result in [16, p.176] which claims that for any given and satisfying , there exists such that , on , and .

We first note that if has mean value zero, then there exists such that and by taking in the above result.

Consider now a decomposition with where is the Lebesgue measure of , and where is the component of with mean-value zero. It is clear that . As just noted, there exists such that and . Since is a polygonal/polyhedral domain, one can find a function whose support resides on , , and . We can find , a renormalization of , satisfying , and it holds that

 \normϕ0H12(∂Ω)≤\normvmL1(Ω)c0≤C\normvmL2(Ω),

with a constant that depends on and . Applying the result in [16, p.176] there exists such that , on , and . Then is a function in satisfying the desired conditions.

###### Theorem 1 (inf-sup condition)

Suppose that . If on and on , and for every , with positive constants and independent of , then there exists a positive constant independent of such that

 inf(p,z,¯z)∈Xhsup(r,v,¯v)∈XhBh(p,z,¯z;r,v,¯v)\norm[0](p,z,¯z)Xh\norm[0](r,v,¯v)Xh≥β>0. (13)

###### Proof

We first prove the following weak inf-sup condition: There exist independent of such that

 inf(z,¯z)∈Wh×Mhsupr∈Vhbh(z,¯z;r)\normrVh≥C0\normzΩ−C1⟨|τ|(z−¯z),z−¯z⟩12∂Th. (14)

We show this by proving an equivalent condition, i.e., there exist independent of such that for any there exists such that and

 bh(z,¯z;r)≥C′0\normz2Ω−C′1⟨|τ|(z−¯z),z−¯z⟩12∂Th\normzΩ. (15)

By Lemma 2 there exists with such that and with depending only on and . We use to denote the element-wise orthogonal projection of into . It then holds that

 \del[0]w,rK=\del[0]Πhw,rK∀r∈Pk−1(K;Rd),

for all . Then

 −\normz2Ω=\del[0]z,∇⋅wΩ=−\del[0]∇z,wΩ+⟨z,w⋅n⟩∂Th=−\del[0]∇z,ΠhwΩ+⟨z−¯z,w⋅n⟩∂Th=−bh\del[0]z,¯z;Πhw+⟨z−¯z,(w−Πhw)⋅n⟩∂Th, (16)

where the third equality holds because and are single-valued on facets, on by the definition of , and on .

The element-wise trace inequality and the Cauchy-Schwarz inequality give

 (17)

Combining Eq. 16 and Eq. 17 we obtain

 bh\del[0]z,¯z;Πhw≥\normz2Ω−C′⟨|τ|(z−¯z),z−¯z⟩12∂Th\normzΩ. (18)

Moreover,

 \normΠhwΩ≤\normwΩ≤\normw1,Ω≤CΩ\normzΩ (19)

holds, so the weak inf-sup condition Eq. 15 follows.

To prove Eq. 13 suppose that is an element satisfying Eq. 15 for given . For given we take , , in Eq. 12, with to be determined later. Then by Eq. 15 and Young’s inequality,

 Bh(p,z,¯z;r,v,¯v) =\normp2Vh+ε\del[0]σ−1p,r0Ω+εbh(z,¯z;r0) +⟨|τ|(z−¯z),z−¯z⟩∂Th ≥\normp2Vh−12\normp2Vh−12ε2(C′2)2\normz2Ω +εC′0\normz2Ω−εC′1⟨|τ|(z−¯z),z−¯z⟩12∂Th\normzΩ +⟨|τ|(z−¯z),z−¯z⟩∂Th ≥\normp2Vh−12\normp2Vh−12ε2(C′2)2\normz2Ω +εC′0\normz2Ω−12⟨|τ|(z−¯z),z−¯z⟩∂Th −12ε2(C′1)2\normz2Ω+⟨|τ|(z−¯z),z−¯z⟩∂Th.

Choosing sufficiently small, we obtain

 Bh(p,z,¯z;r,v,¯v)≥C\del[1]\normp2Vh+\normz2Ω+⟨|τ|(z−¯z),z−¯z⟩∂Th.

Finally, it is not difficult to show, by the choice of , that

 \norm(r,v,¯v)Xh≤C\norm(p,z,¯z)Xh.

This proves Eq. 13.

## 4 Error analysis

In this section we present the a priori error estimates of our HDG method. For this we recall an interpolation result (cf. [10, Theorem 2.1]): For suppose that is a piecewise constant function which is nonnegative if and nonpositive if such that . Then there exists a bounded interpolation operator from to such that equationparentequation

 \del[0]ΠVr,~rK =\del[0]r,~rK ∀~r∈\sbr[0]Pk−1(K)d, (20a) \del[0]ΠWv,~vK =\del[0]v,~vK ∀~v∈Pk−1(K), (20b) ⟨ΠVr⋅n+τΠWv,λ⟩F =⟨r⋅n+τv,λ⟩F ∀λ∈Pk(F), (20c)

for and , and

 \normr−ΠVrK ≤Chkr+1K|r|kr+1,K+hkv+1Kτ∗K|v|kv+1,K, (21) \normv−ΠWvK ≤Chkv+1K|v|kv+1,K+hkr+1KτmaxK|∇⋅r|kr,K, (22)

hold with . Here,

 τmaxK={maxF⊂∂Kτ|F if τ≥0,maxF⊂∂K(−τ|F) if τ≤0, τ∗K={maxF⊂∂K∖F∗τ|F if τ≥0,maxF⊂∂K∖F∗(−τ|F) if τ≤0,

where is the face where is attained, and the implicit constants are independent of and . We remark that the proof in [10] is only for nonnegative . The proof, however, can easily be adapted to nonpositive , so we omit its proof here. We will also require , the facet-wise projection from to .

For brevity of notation, we will denote the difference between an unknown and its approximation by . It will be convenient to also split the error into interpolation and approximation errors:

 eq =eIq+ehq:=(q−ΠVq)+(ΠVq−qh), eu =eIu+ehu:=(u−ΠWu)+(ΠWu−uh), (23) e¯u =eI