1 Introduction

# Optimal adaptivity for non-symmetric FEM/BEM coupling

## Abstract.

We develop a framework which allows us to prove the essential general quasi-orthogonality for the non-symmetric Johnson-Nédélec finite element/boundary element coupling. General quasi-orthogonality was first proposed in [8] as a necessary ingredient of optimality proofs and is the major difficulty on the way to prove rate optimal convergence of adaptive algorithms for many strongly non-symmetric problems. The proof exploits a new connection between the general quasi-orthogonality and -factorization of infinite matrices. We then derive that a standard adaptive algorithm for the Johnson-Nédélec coupling converges with optimal rates. The developed techniques are fairly general and can most likely be applied to other problems like Stokes equation.

Supported by the Australian Research Council (ARC) under grant number DE170100222 and by the Austrian Research Fund (FWF) under grant number P27005.

## 1. Introduction

The theory of rate optimal adaptive algorithms for finite element methods originated in the seminal paper [32] by Stevenson and was further improved in [11] by Cascon, Kreuzer, Nochetto, and Siebert. These papers prove essentially, that a standard adaptive algorithm of the form

 \frameboxSolve⟶\frameboxEstimate⟶\frameboxMark⟶\frameboxRefine

generates asymptotically optimal meshes for the approximation of the solution of a Poisson problem. The new ideas sparked a multitude of papers applying and extending the techniques to different problems, see e.g., [26, 12] for conforming methods, [28, 6, 7, 9, 27] for nonconforming methods, [13, 10, 24] for mixed formulations, and [21, 33, 4, 18, 19] for boundary element methods (the list is not exhausted, see also [8] and the references therein). All the mentioned results, however, focus on symmetric problems in the sense that the underlying equation induces a symmetric operator. The first proof of rate optimality for a non-symmetric problem which does not rely on additional assumptions is given in [20] for a general second order elliptic operator with non-vanishing diffusion coefficient of the form

 −div(A∇u)+b⋅∇u+cu=f.

This approach, however, relies heavily on the fact that the non-symmetric part of the operator is only a compact perturbation (one differentiation instead of two for the diffusion part). The present work aims to shed some light on the completely unexplored world of rate optimality for strongly non-symmetric problems (meaning that the non-symmetric part of the operator is not substantially smaller in any sense). Although the work is focused on the particular model problem of Johnson-Nédélec FEM/BEM coupling, we believe that the developed techniques will be very useful for many other non-symmetric problems of the form for a non-symmetric operator and a right-hand side .

For the particular case of FEM/BEM coupling, only convergence of the adaptive algorithm is known. This was first proven rigorously in [2] for the standard residual based error estimator which was first derived in [5].

Using the abstract framework for rate optimality developed in [8], we observe that the major obstacle is the general quasi-orthogonality property introduced in [8]. The property is a generalization of the usual orthogonality property

 ∥u−uℓ+1∥2+∥uℓ+1−uℓ∥2=∥u−uℓ∥2, (1.1)

for increasingly accurate nested Galerkin approximations of the exact solution . The orthogonality (1.1) follows immediately from the well-known Galerkin orthogonality, if the underlying problem induces a symmetric bilinear form and thus a Hilbert (energy-) norm . If the problem is non-symmetric, however, (1.1) fails to hold (even in approximate forms usually called quasi-orthogonality) and this breaks all existing optimality proofs. In [8], we prove that general quasi-orthogonality is the weakest possible orthogonality condition in the sense that it is necessary to prove optimality. While rigorously stated in Section 2.6 below, general quasi-orthogonality roughly implies that the approximation error has a decomposition of the form

 ∥u−uℓ∥2≃∥uℓ+1−uℓ∥2+∥uℓ+1−uℓ+2∥2+…for all ℓ∈N

for nested Galerkin approximations of . While this property seems hard to prove by itself, we discover an interesting connection to the -factorization of infinite matrices in this work.

This connection can be formulated as follows: Assume that there exists a Riesz basis of the underlying Hilbert space such that the problem can be equivalently stated as a matrix equation

 Mx=GwithM∈RN×N,G∈RN,

where for all , , and . If the matrix has an -factorization for lower/upper-triangular infinite matrices such that are bounded operators, then general quasi-orthogonality holds.

We exploit the stated connection by constructing a suitable Riesz basis for the particular FEM/BEM coupling and then proving that a bounded factorization exists. The construction of the Riesz basis is quite challenging, since for the FEM/BEM coupling, we need compatible basis functions in and . This requires an extension of the well-known Scott-Zhang projection to make operators on different level commute with each other. To prove that the resulting matrix has the desired -factorization, we rely on techniques from wavelet methods, which prove that is exponentially decaying (in the sense of Jaffard class matrices). While being probably an artifact of the proof, we are forced to introduce a grading condition on the adaptively generated meshes, as was also done in [16].

This opens the door to prove rate optimality of the adaptive algorithm. It turns out that all the other requirements for rate optimality (formulated in [8]) can be shown by combining arguments from FEM and BEM (for both methods, rate optimality has been proved already).

The strategy of the proof is very general, but the details a tailored to the present FEM/BEM coupling problem. We are confident that similar techniques can be used to prove optimality of adaptive algorithms for the Stokes equation and many other related non-symmetric problems. The author would like to note that, to the best of his knowledge, the theory of -factorization of infinite matrices currently cannot answer very interesting (and for this work very useful) questions like: Which positive definite matrices have a bounded factorization? This is the reason why this paper is quite technical despite the simple underlying idea. It is possible that advances in this direction could improve (e.g., by removing the grading condition) and simplify the present result.

As an interesting side result, the Riesz basis constructed in Section 6 below can, in principle, be computed and used for actual implementations. This brings the benefit of uniformly bounded condition numbers of the involved matrices without preconditioning.

### 1.1. Outline of the paper

The main result is given in Theorem 7.1 in Section 7. In Section 3.1, we introduce the class of Jaffard matrices, and show that they admit a bounded -factorization under certain conditions. In Section 3.2, we show that general quasi-orthogonality is equivalent to the fact that a certain (infinite) system matrix of the problem at hand has a bounded -factorization. This observation is the key element of the paper. The remainder of the work is devoted to building a system matrix for the Johnson-Nédélec coupling, which fits into this framework of Jaffard class matrices. Therefore, we construct a local wavelet basis in Section 6. To that end, we use a new quasi-interpolation operator from Section 5 which is based on the classical Scott-Zhang projection. In Section 4 we construct certain metrics which characterize the exponential decay of the system matrix. Finally, Section 7.5 constructs the system matrix.

### 1.2. Notation

We use to denote the cardinality of a set . Moreover, denotes the identity matrix , for all and for all . The dimension is only specified when not clear from the context. The standard space of squared summable sequences is denoted by . We denote the -norm by , whereas the operator norm for operators on is denoted by . Operators on are often identified with infinite matrices and we use the norms and .

## 2. General assumptions

### 2.1. Preliminaries

In the following, is a polygonal domain with boundary . Given a Lipschitz domain , we denote by the usual Sobolev spaces for . For non-integer values of , we use real interpolation to define . Their dual spaces are defined by extending the -scalar product. Given , we define as the trace space of for all . Again, the dual space is defined via the extended -scalar product.

Remark.  There is no reason for the author to believe that the methods developed in this work are restricted to the 2D case. However, the technical difficulties are already substantial for and thus we decided to restrict to this case for clarity of presentation. ∎

### 2.2. Variational form

The main goal of this paper is to prove optimality of FEM/BEM coupling. However, most of the methods work in a much broader context. Therefore, we start with an abstract variational problem and go back to the concrete application in Section 7. To that end, suppose is a separable Hilbert space. Moreover, suppose that is a nested sequence of subspaces, i.e.,

 Xℓ⊆Xk⊆Xfor all ℓ≤k∈N.

Assume that is a bounded bilinear form, which is additionally elliptic, i.e.,

 infx∈Xa(x,x)∥x∥2X=c0>0. (2.1)

For , define and for all as the unique solutions of

 a(u,v)=f(v)for all v∈Xanda(uℓ,v)=f(v)for all v∈Xℓ. (2.2)

We further assume that is a space of functions on the domain and that the subspaces correspond to some triangulations discussed in detail below.

### 2.3. Mesh refinement

Let be a triangulation of into compact triangles which resolves the corners of . Given two triangulations , we write for some if is generated from by refinement of all via newest vertex bisection. We write if is generated from by a finite number of iterated newest-vertex-bisection refinements and we denote the set of all possible refinements by . Given , we call a local refinement of , if there exists such that . Given for some , denotes the number of bisections necessary to generate from a parent element in .

We define as the set of nodes of and as the set of edges of .

We define as the mesh-size function by for all .

Given , we define the patch

 ω(T,T):={T′∈T:T∩T′≠∅}.

Given a subset , we define the patch

 ω(Ω′,T):={T∈T:(⋃ω(T,T))∘∩Ω′≠∅},

where denotes the interior of a set. Note that in case of , the two definitions coincide. The extended patches are defined iteratively by

 ω1(Ω′,T):=ω(Ω′,T),andωk(Ω′,T):=ω(⋃ωk−1(Ω′,T),T).
###### Definition 2.1.

We consider an auxiliary sequence of uniform refinements such that and

 Extra open brace or missing close brace

which means that each element of is bisected -times to obtain . There exist constants which depend on and on such that

 C−1baseC−ℓmesh≤diam(T)≤CbaseC−ℓmesh

for all and all . We choose sufficiently large such that , where is defined in Lemma 5.5 below.

Given a triangulation , we assume that we can compute an error estimator . In the application to the FEM-BEM coupling below, we have to restrict to adaptive triangulations with mild grading in the sense that there exists such that

This condition is necessary for the present proof and also appears in [16] to prove optimal convergence in the -norm. By , we denote all triangulations which satisfy (2.3) for a given . Lemma 2.3 below shows that the restriction does not alter the optimal convergence rate. Numerical experiments (see, e.g., [2]) suggest that the restriction is not even necessary for optimal convergence rate, and thus might just be an artifact of the proof. In the following, we assume that is sufficiently large to satisfy all the conditions in the proofs below.

We assume that the sequence is generated by an adaptive algorithm of the form

###### Algorithm 2.2.

Input: , , , , .
For do:

1. Compute .

2. Compute error estimator for all .

3. Mark set of minimal cardinality such that

 ∑T∈MℓηT(Tℓ)≥θ∑T∈TℓηT(Tℓ). (2.4)
4. Refine at least the elements of to obtain .

5. Refine additional elements to ensure that satisfies (2.3).

Output: sequence of meshes and corresponding solutions .

### 2.5. Rate optimality

We aim to analyze the best possible algebraic convergence rate which can be obtained by the adaptive algorithm. This is mathematically characterized as follows: For the exact solution , we define an approximation class by

 u∈Asdef.⟺∥u∥As:=supN∈NminT∈T#T−#T0≤NNsη(T)<∞. (2.5)

By definition, a convergence rate is theoretically possible if the optimal meshes are chosen. In view of mildly graded triangulations, we define

In Lemma (2.3) below, we show that in many situations . In the spirit of [8], rate optimality of the adaptive algorithm means that there exists a constant such that

 C−1opt∥u∥As≤supℓ∈N0η(Tℓ)(#Tℓ−#T0+1)−s≤Copt∥u∥As,

for all with .

### 2.6. The Axioms

To formulate the axioms below, we define for a given triangulation the corresponding space as well as the discrete solution of

 a(uT,v)=⟨f,v⟩for all v∈XT.

As proved in [8], we need to check the axioms (A1)–(A4) to ensure rate optimality for a given adaptive algorithm: There exist constant , , , , , and such that

1. Stability on non-refined elements: For all refinements of a triangulation , for all subsets of non-refined elements, it holds that

 ∣∣(∑T∈SηT(ˆT)2)1/2−(∑T∈SηT(T)2)1/2∣∣≤Cstab∥uT−uˆT∥X.
2. Reduction property on refined elements: Any refinement of a triangulation satisfies

 ∑T∈ˆT∖TηT(ˆT)2≤qred∑T∈T∖ˆTηT(T)2+Cred∥uT−uˆT∥2X.
3. General quasi-orthogonality: For one sufficiently small the output of Algorithm 2.2 satisfies for all

 ℓ+N∑k=ℓ(∥uk+1−uk∥2X−Cqoε∥u−uk∥2X)≤Cqo∥u−uℓ∥2X. (2.6)
4. Discrete reliability: For all refinements of a triangulation , there exists a subset with and such that

 ∥uˆT−uT∥2X≤C2dlr∑T∈R(T,ˆT)ηT(T)2.

In some situations, it might even be possible to prove the stronger form of (A3), namely

 C−1qo∥u−uℓ∥2X≤∞∑k=ℓ∥uk+1−uk∥2X≤Cqo∥u−uℓ∥2X\rm for all ℓ∈N. (2.7)

The main obstacle is the general quasi-orthogonality (A3) and its proof for the particular FEM/BEM coupling below takes up the majority of this work. The other axioms follow from the combination of techniques for FEM and BEM.

###### Lemma 2.3.

Under (A1)–(A2) and (A4), there holds .

###### Proof.

The assumption in [16, Section A.3] and (2.3) are equivalent up to shape regularity. The definition of and implies . Let being generated from by iterated refinements for . Since every newest-vertex bisection refinement generates at least two sons, we have . The result [16, Theorem 4] shows that by refining all elements in the and additionally making sure that (2.3) holds, we find with . The constant depend only on and . From [8, Lemma 3.4], we see that (A1)–(A2)&(A4) imply quasi-monotonicity of in the sense

This shows

and concludes the proof. ∎

## 3. General quasi-orthogonality and Lu-factorization

In this section, we establish the link between general quasi-orthogonality (A3) and -factorization of infinite matrices. To that end, we first introduce exponentially decaying matrices.

### 3.1. Jaffard class matrices

Jaffard class matrices generalize the notion of matrices which decay exponentially away from the diagonal. The generalization allows to replace the distance between indices by a general metric . This class was introduced and analyzed in [23].

###### Definition 3.1 (Jaffard class).

We say that an infinite matrix is of Jaffard class, for some metric and some if for all there exists such that

 |Mij|≤C(γ′)exp(−γ′d(i,j))for all i,j∈N. (3.1)

Moreover, the metric must satisfy for all

 supi∈N∑j∈Nexp(−εd(i,j))<∞. (3.2)

We also write to state the existence of parameters such that .

###### Definition 3.2 (banded matrix).

We say that an infinite matrix is banded with respect to some metric if there exists a bandwidth such that

 d(i,j)>b⟹Mij=0for all i,j∈N. (3.3)

In this case, we write . Note that we do not require to satisfy (3.2). We also write or to state that the missing parameters exist.

The following technical lemmas prove some straightforward facts about infinite matrices.

###### Lemma 3.3.

Let , , for some with respect to some metric and respective bandwidths . Then, there holds

 n∑i=1m∏j=1Mi,j∈B(d,m∑j=1bj).
###### Proof.

Obviously, is closed under summation. The definition of the matrix product shows

 (Mi,1Mi,2)ij=∞∑k=1(Mi,1)ik(Mi,2)kj=∑k∈Nd(i,k)≤b1,d(k,j)≤b2(Mi,1)ik(Mi,2)kj,

and hence if . Induction on proves . This concludes the proof. ∎

###### Lemma 3.4.

Let , then is a bounded operator (the modulus is understood entry wise).

###### Proof.

Given , there holds for all with (3.1) and (3.2)

 supi∈N∑j∈N|Mij|≲supi∈N∑j∈Nexp(−γ′d(i,j))<∞.

Analogously, we obtain . The standard interpolation estimate concludes the proof. ∎

###### Lemma 3.5.

Let such that is additionally elliptic, i.e.,

 supx∈ℓ2Mx⋅x∥x∥2ℓ2=Cell>0. (3.4)

Then,

 ¯¯¯¯¯¯M∈RN×Nwith¯¯¯¯¯¯Mij:=supn∈Nn≥max{i,j}|(M|{1,…,n}×{1,…,n})−1|ij

is of exponential class and thus a bounded operator . The constant depends only on , , and , whereas for all , depends only on an upper bound for and on .

###### Proof.

The result [23, Proposition 2] shows that . Inspection of the proof reveals that depends only on , , and depends only on an upper bound for from Definition 3.1 and on . Therefore, we have for all

 |¯¯¯¯¯¯Mij|≤˜C(γ′)exp(−γ′d(i,j))for all i,j∈N

and hence . Lemma 3.4 concludes the proof. ∎

### 3.2. LU-factorization

We say that a matrix has an -factorization if for matrices such that

 Lij=0 and Uji=0for all i,j∈N,i≤j.

Given a block structure in the sense that there exist numbers with and for all , we denote matrix blocks by

 M(i,j):=M|{ni,…,ni+1−1}×{nj,…,nj+1−1}∈R(ni+1−ni)×(nj+1−nj).

By , we denote the restriction of to the first blocks.

###### Lemma 3.6.

Let such that is bounded and elliptic (3.4). Moreover, let . Assume a block structure . Then, given , there exists a bandwidth such that for all , there exists such that

 ∥M−1−R∥2+supk∈N∥M[k]−1−Rk∥2≤ε.

If is additionally block-banded in the sense for all , then, and will additionally be block-banded with bandwidth . If is block-diagonal, also and will be block-diagonal. The bandwidth depends only on , , , and .

###### Proof.

Let or . Since is elliptic with constant , we obtain for and

 ∥x−αAx∥2ℓ2 =∥x∥2ℓ2−2α⟨x,Ax⟩ℓ2+α2∥Ax∥2ℓ2 ≤(1−2αCell+α2∥A∥22)∥x∥2ℓ2≤(1−C2ell/∥M∥22)∥x∥2ℓ2.

This shows . We obtain

 A−1=α(αA)−1=α(I−(I−αA))−1=α∞∑k=0(I−αA)k.

We define () by for some such that . Since , Lemma 3.3 shows that is banded as well. The bandwidth depends only on , , and . If is additionally block-banded, also and will be block-banded with bandwidth . Hence will be block-banded with bandwidth . The same argumentation proves the statement for block-diagonal . This concludes the proof. ∎

The following results prove that block-banded matrices hand down some structure to their -factors. This is used in Section 7 to construct suitable hierarchical bases for FEM/BEM coupling.

###### Lemma 3.7.

Let such that is bounded and is elliptic (3.4). Assume a block structure in the sense that there exist numbers with and for all . Moreover, let be block-banded in the sense for for some . Then, the block--factorization for block-upper/block-lower triangular matrices such that for all exists, is block-banded with bandwidth , and satisfies

 ∥L∥2+∥U∥2+∥L−1∥2+∥U−1∥2<∞.

Moreover, the block-diagonal matrix , is bounded and elliptic (3.4) with bounded and elliptic inverse.

###### Proof.

Since has invertible principal sub matrices, it is well-known that the (block)-LU-factorization exists. The triangular structure of the factors implies that also has a block--factorization. Since and are invertible by definition, is invertible for all . The block-triangular structure guarantees that , , and hence existence of as matrices in . Moreover, it is well-known that and are block-banded with bandwidth . By definition, we have . Since is lower-triangular with normalized block-diagonal (only identities in the diagonal blocks), the same is true for . Therefore, we obtain

 M[k]−1(i,k)=k∑r=1U[k]−1(i,r)L[k]−1(r,k)=U[k]−1(i,k)=U−1(i,k), (3.5)

where the last identity follows from the fact that is upper-block triangular.

Since is bounded and elliptic (3.4), also is bounded and elliptic. Therefore, we see that (since is block-triangular) is bounded and elliptic and thus conclude that and also are bounded and elliptic. Moreover, we see that . Hence, there holds

 supi,j∈N∥L(i,j)∥2≤supi,j∈N∑|k−i|≤b0∥M(i,k)U−1(k,j)∥2≤2b0∥M∥2supi,j∥U−1(i,j)∥2<∞.

Since is block-banded with bandwidth , Lemma 8.4 shows . This implies .

Let be the analogous block-LU-factorization for the transposed matrix (note that is still elliptic, bounded, and banded). Since normalized -factorizations are unique, we see that

 ˜L=UTD−Tand˜U=DTLT (3.6)

Repeating the above arguments shows . With boundedness and ellipticity of , (3.6) shows . Altogether, we prove the statement and conclude the proof. ∎

###### Lemma 3.8.

Under the assumptions of Lemma 3.7, assume that additionally . Given , there exists and block-upper triangular which is additionally block-banded in the sense for such that

 ∥U−1−U−1ε∥2≤ε.

The approximation is invertible with bounded inverse such that . Moreover, there exists block-diagonal which is bounded and elliptic such that .

###### Proof.

Lemma 3.6 shows that there exist which are block-banded with bandwidth such that for all . Since is bounded and elliptic, also is bounded and elliptic (3.4). Choosing sufficiently small, we ensure that also and are bounded an elliptic with uniform constant.

Inspired by (3.5), we define a first approximation to by

 T(i,j):={0for j

This ensures that and that is block-banded with bandwidth . Additionally, we obtain

 supi,j∈N∥T(i,j)−U−1(i,j)∥2≲ε. (3.8)

We define an approximation to (which is block-banded with bandwidth ) by

 S(i,j):=⎧⎨⎩0for ji,Ifor j=i,(MT)(i,j)for i−b0≤j

The definition and (3.8) imply

 ∥L(i,j)−S(i,j)∥2 ≤∥i+b0∑k=i−b0M(i,k)(U−1(k,j)−T(k,j))∥2 (3.9) ≤i+b0∑k=i−b0∥M∥2∥U−1(k,j)−T(k,j))∥2≲∥M∥2b0ε. (3.10)

Since both and are block-banded with bandwidth , Lemma 8.4 shows even

 ∥L−S∥2≲ε, (3.11)

where the hidden constant is independent of . Moreover, Lemma 3.3 shows that , for some which depends only on and .

Recall from above with , and is block-banded with bandwidth . This allows to define by

 U−1ε:=RS.

We obtain from the definition and with (3.11)

 ∥U−1ε−U