DPG method with optimal test functionsfor a transmission problem Supported by CONICYT through FONDECYT projects 1110324, 3140614 and Anillo ACT1118 (ANANUM).

# DPG method with optimal test functions for a transmission problem ††thanks: Supported by CONICYT through FONDECYT projects 1110324, 3140614 and Anillo ACT1118 (ANANUM).

Norbert Heuer    Michael Karkulik Facultad de Matemáticas, Pontificia Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, Macul, Santiago, Chile, email: {nheuer,mkarkulik}@mat.puc.cl
###### Abstract

We propose and analyze a numerical method to solve an elliptic transmission problem in full space. The method consists of a variational formulation involving standard boundary integral operators on the coupling interface and an ultra-weak formulation in the interior. To guarantee the discrete inf-sup condition, the system is discretized by the DPG method with optimal test functions. We prove that principal unknowns are approximated quasi-optimally. Numerical experiments for problems with smooth and singular solutions confirm optimal convergence orders.

Key words: transmission problem, coupling, DPG method with optimal test functions, ultra-weak formulation, boundary elements, finite elements

AMS Subject Classification: 65N30, 35J20, 65N38

## 1 Introduction

In recent years, Demkowicz and Gopalakrishnan have established the discontinuous Petrov-Galerkin method with optimal test functions as a method that is designed to be stable [17, 16]. In particular, it aims at robust discretizations of singularly perturbed problems, cf. [10, 11, 12, 14, 18]. In general terms, this method consists in applying Petrov-Galerkin approximations with optimal test functions to (in most cases) ultra-weak variational formulations (see [13] for an early use of ultra-weak formulations). Such a variational formulation is obtained by element-wise integration by parts on some partition and the replacement of appearing terms on the elements’ boundaries with new unknowns (see [5] for the idea of introducing independent boundary unknowns).

Until recently the DPG method with optimal test functions has been studied only for problems on bounded domains. In [24], we considered boundary value and screen problems of Neumann type which can be reduced to a hypersingular boundary integral equation. This includes the case of a PDE on an unbounded domain. In this paper we study for the first time a DPG strategy with optimal test functions for a transmission problem. This problem is of academic nature (Poisson equation) and set in the full space. We expect that our fairly general approach is applicable to more practical cases of transmission problems with singular perturbations on bounded subdomains. This is feasible whenever there is a DPG finite element technique available for the corresponding singularly perturbed problem on the subdomain.

Transmission problems often appear in the modeling of multiple physical phenomena, and therefore require the combination of two, possibly different, numerical methods. A popular approach is the coupling of finite elements and boundary elements. Whereas finite elements can be used in a relatively straightforward way to solve nonlinear PDEs with space-dependent coefficients and sources, it is a challenging task to apply them on unbounded domains. On the other hand, boundary elements can deal naturally with unbounded domains but, by construction, work best for linear, homogeneous PDEs with constant coefficients.

There are different approaches for the coupling of finite and boundary elements. In this paper we consider the so-called nonsymmetric or one-equation coupling, also referred to as Johnson-Nédélec coupling. It has the practical advantage of involving only two of the four classical boundary integral operators. This coupling was mathematically analyzed first in [8, 9, 26]. The mathematical proofs in the mentioned works require one of the involved boundary integral operators (the double-layer operator) to be compact. In case of the Laplace equation this property does not hold on polyhedral domains and, in the case of linear elasticity, not even on domains with smooth boundary. In [30], Sayas proved the well-posedness of the Johnson-Nédélec coupling without using compactness arguments and thus gave a mathematical justification of the nonsymmetric coupling even on polyhedral domains. Since then, different authors have re-considered nonsymmetric couplings, see [1, 19, 21, 29, 31]. For an extensive discussion of this topic, we refer to [20].

In this paper, we extend the nonsymmetric coupling of Johnson-Nédélec to a DPG method with optimal test functions. More specifically, given a coupled system of PDEs (one in a bounded domain, one in the exterior), we use an ultra-weak formulation for the interior part coupled to classical boundary integral equations for the exterior problem. The whole system is discretized by the Petrov-Galerkin method with optimal test functions.

An alternative approach would be to couple standard boundary elements with a DPG scheme restricted to the interior problem, i.e., optimal test functions are only used for the interior (finite element) part. In this case, however, several difficulties arise. For instance, it is unclear how to choose the remaining test functions to generate a square system. The design and analysis of a DPG-BEM coupling is left for future research. In contrast, the approach of computing optimal test functions for the whole system (as in this paper) appears to be the most generic one and deserves a thorough analysis.

An outline of this paper is as follows. In Section 2, we introduce the model problem and present the mathematical framework. We also summarize some results related to the DPG method, to be used subsequently. In Subsection 2.4 we formulate the method and state the main results. Theorem 4 establishes stability of the continuous variational formulation and Theorem 5 shows the quasi-optimality of conforming discretizations (so-called Céa-lemma). Proofs of the two theorems are given in Section 4. Principal part of these proofs involves the analysis of adjoint problems. This analysis is made in Section 3. In Section 5 we present some numerical results that underline our theory. Some conclusions are drawn in Section 6.

## 2 Mathematical setting and main results

Let , , be a bounded Lipschitz domain with boundary . Our model transmission problem is as follows: given volume data and jumps , , find and such that

 −div∇u =f in Ω (1a) −div∇uc =0 in Ωc:=Rd∖¯¯¯¯Ω (1b) u−uc =u0 on Γ (1c) ∂∂nΩ(u−uc) =ϕ0 on Γ (1d) uc(x) =O(|x|−1) as |x|→∞. (1e)

The normal vector on points in direction of . Here, , , , and the spaces are of standard Sobolev type (some more details are given in the next section). For , we assume that and . The scaling condition on is to ensure the ellipticity of the single layer operator, and the compatibility condition on the data and is needed in order to use the radiation condition in the form (1e).

### 2.1 Abstract DPG method

We briefly recall the abstract framework of the DPG method with optimal test functions, cf. [16, 17, 34]. We state this in a form that will be convenient for the forthcoming analysis. The continuous framework is provided by the following result which is a consequence of the open mapping theorem and the properties of conjugate operators, cf. [33, Chapters II.5, VII.1]. In this manuscript, all suprema are taken over the indicated sets except .

###### Lemma 1.

Denote by and two reflexive Banach spaces. Let be a bijective and bounded linear operator with conjugate operator . Define

 ∥u∥U,opt:=supv∈V⟨Bu,v⟩∥v∥V and ∥v∥V,opt:=supu∈U⟨Bu,v⟩∥u∥U.

Then both operators and are isomorphisms, and and define norms in and which are equivalent to and , respectively. Furthermore, there holds

 ∥u∥U=supv∈V⟨Bu,v⟩∥v∥V,opt and ∥v∥V=supu∈U⟨Bu,v⟩∥u∥U,opt. (2)

For a proof of (2) we refer to [34, Prop. 2.1]. Now suppose additionally that is a Hilbert space. Given a bilinear form , and a linear functional , we aim to

 find u∈U such that b(u,v)=ℓ(v) for all v∈V. (3)

For an approximation space , the Petrov-Galerkin method with optimal test functions is to

 find uhp∈Uhp such that b(uhp,vhp)=ℓ(vhp)∀vhp∈Θ(Uhp). (4)

Here, is the trial-to-test operator defined by

 ⟨Θu,v⟩V=b(u,v)∀v∈V

and is the inner product in which induces the norm . The following result is a consequence of the Babuška-Brezzi theory [2, 7, 32]. For a proof of the best approximation property see [17, Thm. 2.2].

###### Lemma 2.

Define the operator by and suppose that the assumptions from Lemma 1 hold. Then, (3) and (4) have unique solutions and , respectively. Furthermore, they satisfy

 ∥u∥U,opt≤∥ℓ∥V′ (5)

and

 ∥u−uhp∥U,opt=infu′hp∈Uhp∥u−u′hp∥U,opt. (6)

### 2.2 Sobolev spaces

We use the standard Sobolev spaces , , , , for Lipschitz domains . Vector-valued spaces and their elements will be denoted by bold symbols. In addition, we use spaces on the boundaries of Lipschitz domains . Denoting by the trace operator, we define

 H1/2(∂ω):={γωu:u∈H1(ω)} and its dual H−1/2(∂ω):=[H1/2(∂ω)]′

and equip them with the canonical trace norm and dual norm, respectively. Here, duality is understood with respect to the extended inner product . The inner product will be denoted by . Let denote a disjoint partition of into open Lipschitz sets , i.e., . The set of all boundaries of all elements is the skeleton . By we mean the outer normal vector on for a Lipschitz set . On a partition, we use product spaces and , equipped with respective product norms. The symbols and denote the -piecewise gradient and divergence operators. We use spaces on the skeleton of , namely

 H1/2(S) :={ˆu∈ΠT∈TH1/2(∂T):∃w∈H1(Ω) such that ˆu|∂T=w|∂T∀T∈T}, H−1/2(S) :={ˆσ∈ΠT∈TH−1/2(∂T):∃q∈H(div,Ω) % such that ˆσ|∂T=q⋅nT∀T∈T}.

These spaces are equipped with the norms

 ∥ˆu∥H1/2(S) :=inf{∥w∥H1(Ω):w∈H1(Ω) such % that ˆu|∂T=w|∂T∀T∈T}, ∥ˆσ∥H−1/2(S) :=inf{∥q∥H(div,Ω):q∈H(div,Ω) such that ˆσ|∂T=q⋅n|T∀T∈T}.

Note that we think of the skeleton not as one geometric object, but rather as the set of boundaries of all elements. Consequently, we have defined and not as canonical trace spaces but as product spaces of trace spaces. This subtle difference simplifies the subsequent analysis. For two functions and we use the notation

 ⟨ˆu,ˆσ⟩S:=∑T∈T⟨ˆu|T,ˆσ|T⟩∂T.

Note that for and integration by parts shows that

 ⟨v,ˆσ⟩S=⟨v,ˆσ⟩Γ,

and so the above left-hand side makes sense also for functions and . We use an analogous duality pairing for functions and . For and we define norms of their jumps across by duality,

Here, on . We will need the following estimates.

###### Lemma 3.

There is a constant which only depends on such that

 ∥[vn]∥1/2,S ≤∥v∥H1(T) for all v∈H1(T), (7) ∥[τ⋅n]∥1/2,S ≤∥τ∥H(div,T) for all τ∈H(div,T). (8)

Furthermore, there holds

 ∥ˆu∥H1/2(Γ) ≤∥ˆu∥H1/2(S) for all ˆu∈H1/2(S), (9) ∥ˆσ∥H−1/2(Γ) ≤∥ˆσ∥H1/2(S) for % all ˆσ∈H−1/2(S). (10)
###### Proof.

The estimates (7)–(8) follow straightforwardly by integration by parts, e.g., cf. [16, Section 4.4] for (8). The estimate (9) follows by definition of the norms in and ,

 ∥ˆu∥H1/2(Γ)=inf{∥w∥H1(Γ):w∈H1(Ω) such that ˆu=γ∂Ωw}≤∥ˆu∥H1/2(S).

The estimate (10) follows the same way, using that is equivalently described as the space of normal components on of functions in , and that , cf. [22, Cor. 2.8]. ∎

### 2.3 Boundary integral operators

In order to incorporate the PDE given in the exterior domain , the classical boundary integral operators will be used. The fundamental solution

 G(z):=⎧⎨⎩−12πlog|z| for d=2,14π1|z| for d>2,

of the Laplacian gives rise to the two potential operators and defined by

 ˜Vϕ(x):=∫ΓG(x−y)ϕ(y)dsy, and ˜Kv(x):=∫Γ∂nΩ(y)G(x−y)v(y)dsy for x∈Rd∖Γ.

Then, boundary integral operators are defined as (single layer operator) and (double layer operator) with adjoint . The operators , , , and are bounded, and there holds the representation formula

 V(∂nΩuc)+(1/2−K)(uc|Γ)=0 (11)

for solutions of the exterior PDE (1b). We refer to [15, 25, 27, 28] for proofs and more details regarding boundary integral equations and the above operators.

### 2.4 Nonsymmetric coupling with ultra-weak formulation and main results

The trial space of our variational formulations will be . This space is a Hilbert space with norm

 ∥(σ,u,ˆu,ˆσ)∥2U(T):=∥σ∥2L2(Ω)+∥u∥2L2(Ω)+∥ˆu∥2H1/2(S)+∥ˆσ∥2H−1/2(S).

In addition, we will need the space , which is a Hilbert space with norm

 ∥(σ,u,ˆu,ˆσ)∥2U(Ω) :=∥σ∥2L2(Ω)+∥u∥2L2(Ω)+∥ˆu∥2H1/2(Γ)+∥ˆσ∥2H−1/2(Γ).

Note that the canonical restrictions of and show that can be viewed as an element of . Using this restriction, we can regard as a subspace of . However, is only a seminorm on . The test space of our formulation will be , being a Hilbert space with norm

 ∥v∥2V(T):=∥v∥2H1(T)+∥τ∥2H(div,T)+∥ψ∥2H−1/2(Γ).

In addition, we will need the space , which is a Hilbert space with norm

 ∥v∥2V(Ω):=∥v∥2H1(Ω)+∥τ∥2H(div,Ω)+∥ψ∥2H−1/2(Γ).

Note that . The variational formulation that we will analyze is the following Johnson-Nédélec type coupling: find such that

 (σ,∇Tv)Ω−⟨ˆσ,[v]⟩S =(f,v)Ω (12a) (σ,τ)Ω+(u,divTτ)Ω−⟨ˆu,[τ⋅n]⟩S =0 (12b) ⟨Vˆσ,ψ⟩Γ+⟨(1/2−K)ˆu,ψ⟩Γ (12c)

for appropriate test functions . The equations (12a)– (12b) are obtained by treating the interior PDE (1a) as in the DPG-finite element method, cf. [16], i.e., writing it as a first order system, testing with appropriate functions, integrating by parts piecewise, and replacing the appearing boundary terms by new unknowns and . These new unknowns already involve the interior trace and normal derivative of on , which are coupled to the exterior problem by using the interface conditions (1c)–(1d) in the representation formula (11). In contrast to that, in the classical nonsymmetric coupling, cf. (17) below, the unknowns are and , where is the normal derivative of on .

The bilinear form on the left-hand side of (12) will be called , and the linear form on the right-hand side will be called . We will use two different formulations which differ in the underlying spaces; the one we actually analyze and solve numerically is

 find u∈U(T) such that b(u,v)=ℓ(v) for all v∈V(T). (13)

The second one is only of theoretical interest and will be needed in the proofs of the main theorems, it is

 find u∈U(Ω) such that b(u,v)=ℓ(v) for all v∈V(Ω). (14)

The first result of this work states unique solvability and stability of the variational formulation (13). The proof will be given in Section 4 below.

###### Theorem 4.

For given , , and and given partition of , the variational formulation (13) has a unique solution . Furthermore,

 ∥(σ,u,ˆu,ˆσ)∥U(Ω)≲∥f∥L2(Ω)+∥u0∥H1/2(Γ)+∥ϕ0∥H−1/2(Γ).

The hidden constant in only depends on .

The second main result of this work is the following quasi optimality result of the Petrov-Galerkin method with optimal test functions associated with the norm .

###### Theorem 5.

Suppose that is a discrete subspace. Then, the discrete formulation (4), where , has a unique solution . Furthermore,

 ∥(σ−σhp,u−uhp, ˆu−ˆuhp,ˆσ−ˆσhp)∥U(Ω)≲ inf(σ′hp,u′hp,ˆu′hp,ˆσ′hp)∈Uhp(T)∥(σ−σ′hp,u−u′hp,ˆu−ˆu′hp,ˆσ−ˆσ′hp)∥U(T),

where is the exact solution of (13). The hidden constant in only depends on .

###### Remark 6.

Note that the norm in the stability estimate and on the left-hand side of the quasi-optimality is a weaker norm in . This norm does not control the parts of and on the inner parts of the skeleton , i.e., the parts which are not on the boundary . However, we have full control of the Cauchy data and on the boundary , which are the only ingredients to solve the exterior problem.

For the proofs of Theorems 4 and 5 we will apply the results of Lemmas 1 and 2, hence we have to check the assumptions of Lemma 1. The boundedness of the bilinear form is shown in Lemma 7. The bijectivity of the operator that corresponds to the bilinear form will be proved in Section 4 below. In a first step, this will yield the results of Lemma 2, i.e., stability and quasi-optimality in the norm defined in Lemma 1. To obtain the results in the main theorems, it remains to relate the norm to the norms and . This will be done by characterizing the optimal test norms (optimal to ) and (optimal to ) and by relating them to the norm (optimal to ). These norm equivalences are the topic of Section 3.2.

## 3 Technical results

We start by showing the boundedness of the bilinear form .

###### Lemma 7.

It holds that

 b(u,v) ≲∥u∥U(T)∥v∥V(T) for all u∈U(T),v∈V(T), (15) b(u,v) ≲∥u∥U(Ω)∥v∥V(Ω) for all u∈U(Ω),v∈V(Ω). (16)

The hidden constant only depends on .

###### Proof.

The proof of (15) follows immediately by the Cauchy-Schwarz inequality and standard boundedness properties. Note that by definition of the norms and and (7) and (8), it holds that

 |⟨ˆσ,v⟩S| ≤∥ˆσ∥H−1/2(Γ)∥[vn]∥1/2,S≤∥ˆσ∥H−1/2(Γ)∥v∥H1(T), |⟨ˆu,τ⋅n⟩| ≤∥[τ⋅n]∥1/2,S∥ˆu∥H1/2(S)≤∥τ∥H(div,T)∥ˆu∥H1/2(S).

Furthermore, the boundedness of the operators , and (9), (10) yield

 |⟨Vˆσ,ψ⟩Γ| ≲∥ˆσ∥H−1/2(Γ)∥ψ∥H−1/2(Γ)≤∥ˆσ∥H−1/2(S)∥ψ∥H−1/2(Γ) |⟨(1/2−K)ˆu,ψ⟩Γ| ≲∥ˆu∥H1/2(Γ)∥ψ∥H−1/2(Γ)≤∥ˆu∥H1/2(S)∥ψ∥H−1/2(Γ).

For the proof of (16), note in addition that

 |⟨ˆσ,v⟩S|≤|⟨ˆσ,v⟩Γ|≤∥ˆσ∥H−1/2(Γ)∥v∥H1/2(Γ)≤∥ˆσ∥H−1/2(Γ)∥v∥H1(Ω)

due to the definition of the norm . The part is treated the same way. ∎

### 3.1 Johnson-Nédélec coupling

The aim of this subsection is to show that our new formulation is equivalent to the classical nonsymmetric coupling. As we will see in Section 4, this implies, in particular, injectivity of the operator that corresponds to the bilinear form .

The transmission problem (1) can be written equivalently as: Given , find such that

 (17a) ⟨(1/2−K)u,ψ⟩Γ+⟨Vϕ,ψ⟩Γ =⟨ψ,λ⟩Γ (17b)

for all . Proof of unique solvability is not straightforward as Problem (17) is not elliptic, and was addressed recently in [30] and also in [1, 29, 31]. Following the approach of [1], the following stability result can be shown.

###### Lemma 8.

For given , the variational formulation (17) has a unique solution , and

 ∥u∥H1(Ω)+∥ϕ∥H−1/2(Γ)≲∥G∥L2(Ω)+∥F∥L2(Ω)+∥ρ∥H−1/2(Γ)+∥λ∥H1/2(Γ).
###### Proof.

Denote the bilinear form on the left-hand side of (17) as and the linear functional on the right-hand side as . Consider the problem of finding such that

 (18)

for all . According to [1, Thm. 14], a solution of (17) also solves (18) and vice versa. Furthermore, [1, Thm. 15] states that the bilinear form on the left-hand side of (18) is continuous and elliptic on . The norm of the linear functional on the right-hand side of (18) is bounded by

 ∥G∥L2(Ω)+∥F∥L2(Ω)+∥ρ∥H−1/2(Γ)+∥λ∥H1/2(Γ).

We finish the proof by application of the Lax-Milgram lemma. ∎

The next lemma shows that our new formulation (13) is equivalent to the classical formulation (17).

###### Lemma 9.

Let , , and . Define , , , and . Then there hold the following statements:

• Suppose that is a solution of (17). Then , and is a solution of (13) and of (14).

• Suppose that respectively is a solution of (13), respectively (14). Then is a solution of (17) and as well as , on , respectively on .

###### Proof.

We first show (i). Using in (17a) shows that

 div∇u=−f∈L2(Ω). (19)

Hence, and eventually . This yields . Identity (19) and integration by parts implies (12a) for all . Integration by parts also shows (12b) for all . From (19) and (17a) we conclude that

Using the symmetry of , this leads us to

 ⟨Vϕ,ψ⟩Γ=⟨V∇u⋅n,ψ⟩Γ−⟨Vϕ0,ψ⟩Γ for all ψ∈H−1/2(Γ)

If we plug the last identity into (17b), we obtain exactly (12c) for all . In total, is a solution of (13). Furthermore, it is also a solution of (14). This follows immediately as by the canonical restriction, and as .

Now we show (ii). To that end, denote by a solution of (13). Equation (12b) first shows that and hence , as well as . As , we have , and from (12a) follows (17a). Likewise, (17b) follows from (12c) using . The case that is a solution of (14) follows analogously. ∎

We additionally need the following stronger result, which shows the surjectivity of the operator associated to our bilinear form.

###### Lemma 10.

For every there exists with for all .

###### Proof.

Using the Riesz representation theorem, we write

 ℓ(v)=(S,v)Ω+(S,∇Tv)Ω+(T,τ)Ω+(T,divTτ)Ω+⟨μ,ψ⟩Γ

with and and . According to Lemma 8, there is a unique solution of (17) with right-hand side data , , arbitrary and . From (17a) it follows that with . Now define , , , and . Integration by parts shows

 (σ,∇Tv)Ω−⟨ˆσ,[τ⋅n]⟩S =(S,v)Ω+(S,∇Tv)Ω, (σ,τ)Ω+(u,divTτ)Ω−⟨ˆu,[v]⟩S =(T,τ)Ω+(T,divTτ)Ω.

Furthermore, (17a) shows that . Therefore, by definition of and (17b), we conclude that . We have thus shown that for all . ∎

### 3.2 Bielak-MacCamy coupling and norm equivalences

In order to relate the norm to a norm of our choice, we will investigate norm equivalences in the test spaces. To that end define seminorms in and by

 ∥v∥V(Ω) :=∥τ+∇Tv∥L2(Ω)+∥divTτ∥L2(Ω) +∥(1/2−K′)ψ−τ⋅n∥H−1/2(Γ)+∥Vψ−v∥H1/2(Γ), ∥v∥V(T),opt :=∥τ+∇Tv∥L2(Ω)+∥divTτ∥L2(Ω) +∥[(τ−E(1/2−K′)ψ)⋅n]∥1/2,S+∥[(v−˜Vψ)n]∥1/2,S.

Here, is a bounded and linear extension operator, i.e., . See [22, Cor. 2.8] for an explicit construction of . Equivalence of norms in the test space amounts to an analysis of the adjoint problem. In case of the nonsymmetric coupling, the adjoint problem is the so-called Bielak-MacCamy coupling, which first appeared in [3]. Given , it consists in finding such that

 (∇v,∇u)Ω+⟨(1/2−K′)ψ,u⟩Γ