Error estimates and a two grid scheme for approximating transmission eigenvaluesThis work was supported by the National Science Foundation of China (Grant Nos.11561014, 11201093).

Error estimates and a two grid scheme for approximating transmission eigenvalues††thanks: This work was supported by the National Science Foundation of China (Grant Nos.11561014, 11201093).

YIDU YANG, JIAYU HAN, AND HAI BI School of Mathematics and Computer Science,Guizhou Normal University,GuiYang,550001, China (ydyang@gznu.edu.cn,hanjiayu126@126.com,bihaimath@gznu.edu.cn).
Abstract

In this paper, using the linearization technique we write the Helmholtz transmission eigenvalue problem as an equivalent nonselfadjoint linear eigenvalue problem whose left-hand side term is a selfadjoint, continuous and coercive sesquilinear form. To solve the resulting nonselfadjoint eigenvalue problem, we give an conforming finite element discretization and establish a two grid discretization scheme. We present a complete error analysis for both discretization schemes, and theoretical analysis and numerical experiments show that the methods presented in this paper can efficiently compute real and complex transmission eigenvalues.

Key words. transmission eigenvalues, finite element method, two grid discretization scheme, error estimates

AMS subject classifications. 65N25, 65N30, 65N15

1 Introduction

Transmission eigenvalues not only have wide physical applications, for example, they can be used to obtain estimates for the material properties of the scattering object [6, 7], but also have theoretical importance in the uniqueness and reconstruction in inverse scattering theory [13]. Many papers such as [7, 13, 23, 28, 29] study the existence of transmission eigenvalues, and [7, 8, 15] explore upper and lower bounds for the index of refraction from knowledge of the transmission eigenvalues. The attention from the computational mathematics community is also increasing (see, e.g., [1, 9, 14, 19, 21, 22, 26, 31, 32]).
In this paper, we consider the Helmholtz transmission eigenvalue problem: Find , , such that

 Δw+k2n(x)w=0,   in D, \hb@xt@.01(1.1) Δσ+k2σ=0,   in D, \hb@xt@.01(1.2) w−σ=0,   on ∂D, \hb@xt@.01(1.3) ∂w∂ν−∂σ∂ν=0,   on ∂D, \hb@xt@.01(1.4)

where is a bounded simply connected set containing an inhomogeneous medium, and is the unit outward normal to .
From [8, 29] we know that for , the weak formulation for the transmission eigenvalue problem (LABEL:s1.1)-(LABEL:s1.4) can be stated as follows: Find , , such that

 (1n(x)−1(Δu+k2u),Δv+¯¯¯¯¯k2n(x)v)0=0,   ∀v∈ H20(D), \hb@xt@.01(1.5)

where denotes the inner product. We denote as usual, then (LABEL:s1.5) is a quadratic eigenvalue problem.
In recent years, based on this weak formulation, various numerical methods have been presented to solve the transmission eigenvalue problem. The first numerical treatment appeared in [14] where three finite element methods were proposed for the Helmholtz transmission eigenvalues which have been further developed in [1, 9, 19, 21, 22, 31]. Among them [14, 22, 31] studied the conforming finite element method, [9, 14, 21] mixed finite element methods, [1] spectral-element method and [19] the Galerkin-type numerical method. Inspired by these works, this paper further studies the conforming finite element method and has three features as follows:
(1) A complete error analysis is presented. Due to the fact that the problem is neither elliptic nor self-adjoint, once its error analysis was viewed as a difficult task. Sun [31] first proved an error estimate for the conforming finite element approximation of (LABEL:s1.5). Following him Ji et al. [22] presented an accurate error estimate, but Theorems 1-2 therein are only valid for real eigenvalues and rely on the conditions of Lemma 3.2 in [31] which are not easy to verify in general, especially for multiple eigenvalues. In this paper, we use the linearization technique in [33] to write the weak formulation (LABEL:s1.5) as an equivalent nonselfadjoint linear eigenvalue problem (LABEL:s2.6). Then the conforming finite element approximation eigenpair of (LABEL:s2.6) is exactly the one of (LABEL:s1.5). Fortunately, the left-hand side term of (LABEL:s2.6) is a selfadjoint, continuous and coercive sesquilinear form, thus we can use Babuska-Osborn’s spectral approximation theory [2] to give a complete error analysis for the conforming finite element approximation of (LABEL:s2.6), namely the error analysis of (LABEL:s1.5). The difficulty of the error estimates in low norms lies in the nonsymmetry of the right-hand side term of (LABEL:s2.6) that involves derivatives. We introduce two auxiliary problems and overcome this difficulty by the Nitche technique in a subtle way. Our theoretical results are proved under general conditions and are valid for arbitrary real and complex eigenvalues.
(2) A finite element discretization with good algebraic properties is presented. We give an conforming finite element discretization (see (LABEL:s3.1) or (LABEL:s5.1)) that will lead to a positive definite Hermitian and block diagonal stiff matrix. Then we use this discretization to solve the transmission eigenvalue problem numerically and obtain both real and complex transmission eigenvalues of high accuracy as expected. Similar discretizations already exist in the literatures. Colton et al. [14] established the conforming finite element discretization for the formulation (LABEL:s1.5) (see (4.3) in [14]), by starting with discretizing (LABEL:s1.5) into a quadratic algebra eigenvalue problem and then linearizing it, which is in the opposite order of our treatment. Though our discrete form (LABEL:s5.1) is formally different from the one in [14] for the resulting mass matrix in [14] is positive definite Hermitian and block diagonal, the two discrete forms are equivalent from the linear algebra point of view. Gintides et al. [19] has used square root of matrices, which is a standard tool, and a Hilbert basis in as test functions to discretizite (LABEL:s1.5) into a linear eigenvalue problem (see (2.14) in [19]) and proved the convergence for the approximate eigenvalues.
(3) A two grid discretization scheme is proposed. As we know, it is difficult to solve nonselfadjoint eigenvalue problems in general. So, Sun [31] adopted iterative methods to work on a series of generalized Hermitian problems and finally to solve for the real transmission eigenvalues. Ji et al. [22] combined the iterative methods in [31] with the extended finite element method to establish a multigrid algorithm for solving the real transmission eigenvalues. This paper establishes a two grid discretization scheme directly for the nonselfadjoint eigenvalue problem (LABEL:s2.6) which is suitable for arbitrary real and complex eigenvalues. The two grid discretization scheme is an efficient approach which was first introduced by Xu [34] for nonsymmetric or indefinite problems, and successfully applied to eigenvalue problems later (see, e.g., [17, 24, 35, 36, 37] and the references therein). With the two grid discretization scheme, the solution of the transmission eigenvalue problem on a fine grid is reduced to the solution of the primal and dual eigenvalue problem on a much coarser grid and the solutions of two linear algebraic systems with the same positive definite Hermitian and block diagonal coefficient matrix on the fine grid , and the resulting solution still maintains an asymptotically optimal accuracy. The key to our theoretical analysis is to find a and prove that (see Step 1 of Scheme 4.1) has a positive lower bound uniformly with respect to the mesh size , which is also critical for further establishing multigrid method and adaptive algorithm. We successfully apply spectral approximation theory to solve this problem (see Lemma 4.1 and Remark 4.1).
The conforming finite element discretization is easy to implement under the package of iFEM [11] with MATLAB, and the numerical results indicate that this discretization is efficient for computing transmission eigenvalues. Moreover, we can further improve the computational efficiency by using the two grid discretization scheme.
Regarding the basic theory of finite element methods, we refer to [2, 4, 12, 27].
Throughout this paper, denotes a positive constant independent of the mesh size , which may not be the same constant in different places. For simplicity, we use the symbol to mean that .

2 A linear formulation of (LABEL:s1.5)

In this paper, we suppose that the index of refraction satisfies the following assumption

 1+δ≤infDn(x)≤n(x)≤supDn(x)<∞,

for some constant , although, with obvious changes, the theoretical analysis in this paper also holds for strictly less than .
Inspired by the works in [9, 14, 21], we now apply the linearization technique in [33] to write the weak formulation (LABEL:s1.5) as an equivalent linear formulation.
From (LABEL:s1.5) we derive that

 (1n−1Δu,Δv)0+k2(1n−1u,Δv)0 +k2(Δu,nn−1v)0+k4(nn−1u,v)0=0,   ∀v∈ H20(D). \hb@xt@.01(2.1)

Introduce an auxiliary variable

 ω=k2u, \hb@xt@.01(2.2)

then

 (ω,z)0=k2(u,z)0,   ∀z∈L2(D). \hb@xt@.01(2.3)

Thus, combining (LABEL:s2.1) and (LABEL:s2.3), we arrive at a linear formulation: Find and nontrivial such that

 (1n−1Δu,Δv)0=−k2(1n−1u,Δv)0 −k2(Δu,nn−1v)0−k2(nn−1ω,v)0,   ∀v∈ H20(D), \hb@xt@.01(2.4) (ω,z)0=k2(u,z)0,   ∀z∈L2(D). \hb@xt@.01(2.5)

This is a nonselfadjoint linear eigenvalue problem.
Suppose that is an eigenpair of (LABEL:s1.5), then it is obvious that is an eigenpair of (LABEL:s2.4)-(LABEL:s2.5). On the other hand, suppose that satisfies (LABEL:s2.4)-(LABEL:s2.5), from (LABEL:s2.5) we get , and substituting it into (LABEL:s2.4) we get (LABEL:s2.1)(or (LABEL:s1.5)). The above argument indicates that (LABEL:s2.4)-(LABEL:s2.5) and (LABEL:s1.5) are equivalent.
Let be the “negative space”, with norm given by

 ∥f∥−l=sup0≠v∈Hl0(D)|(f,v)0|∥v∥l,   l=1,2.

Define the Hilbert space with norm , and define with norm , , . It’s obvious that compactly (see[5, 35] ).
Denote . Let

 A((u,ω),(v,z))=(1n−1Δu,Δv)0+(ω,z)0, B((u,ω),(v,z))=−(1n−1u,Δv)0−(Δu,nn−1v)0−(nn−1ω,v)0+(u,z)0,

then (LABEL:s2.4)-(LABEL:s2.5) can be rewritten as: Find , such that

 A((u,ω),(v,z))=λB((u,ω),(v,z)),   ∀(v,z)∈H. \hb@xt@.01(2.6)

Thus we get the following theorem.

Theorem 2.1

If is an eigenpair of (LABEL:s1.5) and , then satisfies (LABEL:s2.6), and if satisfies (LABEL:s2.6), then is an eigenpair of (LABEL:s1.5) and .

By calculation we derive for any ,

 A((u,ω),(v,z))=∫D1n−1Δu¯¯¯¯¯¯¯Δv+ω¯¯¯zdx=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∫D1n−1Δv¯¯¯¯¯¯¯¯Δu+z¯¯¯ωdx=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯A((v,z),(u,ω)), |A((u,ω),(v,z))|≲∥Δu∥0∥Δv∥0+∥ω∥0∥z∥0≲∥(u,ω)∥H∥(v,z)∥H, A((u,ω),(u,ω))=∫D1n−1Δu¯¯¯¯¯¯¯¯Δu+ω¯¯¯ωdx≳∥(u,ω)∥2H,

i.e., is a selfadjoint, continuous and coercive sesquilinear form on .
We use and as an inner product and norm on , respectively. Then is equivalent to the norm .
When , , , we deduce

 |B((f,g),(v,z))|=|−(1n−1f,Δv)0−(Δf,nn−1v)0−(nn−1g,v)0+(f,z)0| ≲∥1n−1f∥0∥v∥2+∥nn−1Δf∥−1∥v∥1+∥nn−1g∥−1∥v∥1+∥f∥0∥z∥0 ≲(∥f∥0+∥nn−1Δf∥−1+∥nn−1g∥−1)∥(v,z)∥H. \hb@xt@.01(2.7)

When , , we have

 |B((f,g),(v,z))|=|(∇(1n−1f),∇v)0+(∇f,∇(nn−1v))0−(nn−1g,v)0+(f,z)0| ≲∥f∥1∥v∥1+∥f∥1∥v∥1+∥g∥−1∥v∥1+∥f∥1∥z∥−1 ≲∥(f,g)∥H1∥(v,z)∥H. \hb@xt@.01(2.8)

And when , , we have

 |B((f,g),(v,z))|=|−(1n−1f,Δv)0−(f,Δ(nn−1v))0−(nn−1g,v)0+(f,z)0| ≲∥f∥0∥v∥2+∥f∥0∥v∥2+∥g∥−2∥v∥2+∥f∥0∥z∥0 ≲∥(f,g)∥H0∥(v,z)∥H. \hb@xt@.01(2.9)

We can see from (LABEL:s2.7)-(LABEL:s2.9) that for any given (), is a continuous linear form on .
The source problem associated with (LABEL:s2.6) is given by: Find such that

 A((ψ,φ),(v,z))=B((f,g),(v,z)),   ∀(v,z)∈H. \hb@xt@.01(2.10)

From the Lax-Milgram theorem we know that (LABEL:s2.10) admits a unique solution. Therefore, we define the corresponding solution operators (s=0,1,2) by

 A(T(f,g),(v,z))=B((f,g),(v,z)),   ∀(v,z)∈H. \hb@xt@.01(2.11)

Then (LABEL:s2.6) has the equivalent operator form:

 T(u,ω)=λ−1(u,ω). \hb@xt@.01(2.12)

Namely, if is an eigenpair of (LABEL:s2.6), then is an eigenpair of (LABEL:s2.12). Conversely, if is an eigenpair of (LABEL:s2.12), then is an eigenpair of (LABEL:s2.6).

Theorem 2.2

Suppose , then is compact; suppose , then is compact.

Proof. When , let in (LABEL:s2.11), then from (LABEL:s2.7) we have

 ∥T(f,g)∥2H≲A(T(f,g),T(f,g))=B((f,g),T(f,g)) ≲(∥f∥0+∥nn−1Δf∥−1+∥nn−1g∥−1)∥T(f,g)∥H,

thus

 ∥T(f,g)∥H≲∥f∥0+∥nn−1Δf∥−1+∥nn−1g∥−1,   ∀(f,g)∈H. \hb@xt@.01(2.13)

Let be any bounded set in , because of the compact embedding and , then is relatively compact in , is relatively compact in and is relatively compact in . And from (LABEL:s2.13) we conclude that is compact.
When , in (LABEL:s2.11), let , then from (LABEL:s2.8) and (LABEL:s2.9) we have

 ∥T(f,g)∥2H≲A(T(f,g),T(f,g))=B((f,g),T(f,g))≲∥(f,g)∥Hs∥T(f,g)∥H,

thus

 ∥T(f,g)∥H≲∥(f,g)∥Hs,   ∀(f,g)∈Hs (s=0,1), \hb@xt@.01(2.14)

which implies that is continuous. Due to the compact embedding , is compact.
Consider the dual problem of (LABEL:s2.6): Find , such that

 A((v,z),(u∗,ω∗))=¯¯¯¯¯λ∗B((v,z),(u∗,ω∗)),   ∀(v,z)∈H. \hb@xt@.01(2.15)

Define the corresponding solution operators by

 A((v,z),T∗(f,g))=B((v,z),(f,g)),   ∀(v,z)∈H. \hb@xt@.01(2.16)

Then (LABEL:s2.15) has the equivalent operator form:

 T∗(u∗,ω∗)=λ∗−1(u∗,ω∗). \hb@xt@.01(2.17)

It can be proved that is the adjoint operator of in the sense of inner product . In fact, from (LABEL:s2.11) and (LABEL:s2.16) we have

 A(T(f,g),(v,z))=B((f,g),(v,z)) =A((f,g),T∗(v,z)),   ∀(f,g),(v,z)∈H. \hb@xt@.01(2.18)

Hence the primal and dual eigenvalues are connected via .

3 The H2 conforming finite element approximation and its error estimates

Let be a shape-regular grid of with mesh size and be a piecewise polynomial space defined on ; for example, is the finite element space associated with one of the Argyris element, the Bell element or the Bogner-Fox-Schmit element (BFS element) (see [12]). Thanks to (LABEL:s2.2), we choose . Then be a conforming finite element space. For the three finite element spaces mentioned above, since are dense in both and , it is obvious that the following condition (C1) holds:
(C1)  If , then as ,

 inf(v,z)∈Hh∥(ψ,φ)−(v,z)∥H→0.

From the operator interpolation theory (see e.g. [4]), the following (C2) is also valid:
(C2)  If , then

 infv∈Sh∥ψ−v∥s≲h2+r−s∥ψ∥2+r,   s=0,1,2.

Throughout this paper, we suppose that (C1) is valid.
The conforming finite element approximation of (LABEL:s2.6) is given by the following: Find , such that

 A((uh,ωh),(v,z))=λhB((uh,ωh),(v,z)),   ∀(v,z)∈Hh. \hb@xt@.01(3.1)

Consider the approximate source problem: Find such that

 A((ψh,φh),(v,z))=B((f,g),(v,z)),   ∀(v,z)∈Hh. \hb@xt@.01(3.2)

We introduce the corresponding solution operator: (s=0,1):

 A(Th(f,g),(v,z))=B((f,g),(v,z)),   ∀(v,z)∈Hh. \hb@xt@.01(3.3)

Then (LABEL:s3.1) has the operator form:

 Th(uh,ωh)=λ−1h(uh,ωh). \hb@xt@.01(3.4)

Define the projection operators and by

 (1n−1Δ(u−P1hu),Δv)0=0,   ∀v∈Sh, \hb@xt@.01(3.5) (ω−P2hω,z)0=0,   ∀z∈Sh. \hb@xt@.01(3.6)

Let

 Ph(u,ω)=(P1hu,P2hω),   ∀(u,ω)∈H.

Then , and

 A((u,ω)−Ph(u,ω),(v,z))=A((u,ω)−(P1hu,P2hω),(v,z)) =(1n−1Δ(u−P1hu),Δv)0+(ω−P2hω,z)0=0,   ∀(v,z)∈Hh, \hb@xt@.01(3.7)

i.e., is the Ritz projection.
Combining (LABEL:s3.7), (LABEL:s2.11) with (LABEL:s3.3), we deduce for any that

 A(PhT(u,ω)−Th(u,ω),(v,z))=A(PhT(u,ω)−T(u,ω),(v,z)) +A(T(u,ω)−Th(u,ω),(v,z))=0,   ∀(v,z)∈Hh,

thus we get

 Th=PhT. \hb@xt@.01(3.8)
Theorem 3.1

Let , then

 ∥T−Th∥H→0, \hb@xt@.01(3.9)

and let , then

 ∥T−Th∥Hs→0,   s=0,1. \hb@xt@.01(3.10)

Proof. When , for any , from (C1) we have

Since is compact, is relatively compact. Thus by the definition of operator norm we have

 ∥T−Th∥H=sup(f,g)∈H,∥(f,g)∥H=1∥(T−Th)(f,g)∥H =sup(f,g)∈H,∥(f,g)∥H=1∥(I−Ph)T(f,g)∥H→0.

When , for , we have

 ∥(I−Ph)T(f,g)∥Hs≲∥(I−Ph)T(f,g)∥H→0.

Since is compact, is relatively compact, thus we have

 ∥T−Th∥Hs=sup(f,g)∈Hs,∥(f,g)∥Hs=1∥(T−Th)(f,g)∥Hs =sup(f,g)∈Hs,∥(f,g)∥Hs=1∥(I−Ph)T(f,g)∥Hs→0.

The proof is completed.

The conforming finite element approximation of (LABEL:s2.15) is given by: Find , such that

 A((v,z),(u∗h,ω∗h))=¯¯¯¯¯¯λ∗hB((v,z),(u∗h,ω∗h)),   ∀(v,z)∈Hh. \hb@xt@.01(3.11)

Define the solution operator satisfying

 A((v,z),T∗h(f,g))=B((v,z),(f,g)),   ∀ (v,z)∈Hh. \hb@xt@.01(3.12)

(LABEL:s3.11) has the following equivalent operator form

 T∗h(u∗h,ω∗h)=λ∗−1h(u∗h,ω∗h). \hb@xt@.01(3.13)

It can be proved that is the adjoint operator of in the sense of inner product . In fact, from (LABEL:s3.3) and (LABEL:s3.12) we have

 A(Th(f,g),(v,z))=B((f,g),(v,z)) =A((f,g),T∗h(v,z)),   ∀(f,g),(v,z)∈Hh. \hb@xt@.01(3.14)

Hence, the primal and dual eigenvalues are connected via .
We need the following lemma (see Lemma 5 on page 1091 of [18]).

Lemma 3.2

Let . Let be an enumeration of the eigenvalues of , each repeated according to its multiplicity. Then there exist enumerations of the eigenvalues of , with repetitions according to multiplicity, such that  ().

In this paper we suppose that and satisfy the above lemma, and let be the kth eigenvalue with the algebraic multiplicity and the ascent , . Since , eigenvalues of (LABEL:s3.1) will converge to .
Let be the spectral projection associated with and , then is the space of generalized eigenvectors associated with and , where denotes the range and denotes the null space. Let be the spectral projection associated with and the eigenvalues ; let be the space of generalized eigenvectors associated with and , and let , then if is small enough. In view of the dual problem (LABEL:s2.15) and (LABEL:s3.11), the definitions of , , , , and are analogous to , , , , and (see [2]).
Given two closed subspaces and , denote

 δ(R,U)=sup(u,ω)∈R∥(u,ω)∥H=1inf(v,z)∈U∥(u,ω)−(v,z)∥H, θ(R,U)s=sup(u,ω)∈R∥(u,ω)∥Hs=1inf(v,z)∈U∥(u,ω)−(v,z)∥Hs,   s=0,1.

We define the gaps between and in as

 ˆδ(R(E),R(Eh))=max{δ(R(E),R(Eh)),δ(R(Eh),R(E))},

and in as

 ˆθ(R(E),R(Eh))s=max{θ(R(E),R(Eh))s,θ(R(Eh),R(E))s}.

Define

 εh(λ)=sup(u,ω)∈R(E)∥(u,ω)∥H=1inf(v,z)∈Hh∥(u,ω)−(v,z)∥H, ε∗h(λ∗)=sup(u∗,ω∗)∈R(E∗)∥(u∗,ω∗)∥H=1inf(v,z)∈Hh∥(u∗,ω∗)−(v,z)∥H.

It follows directly from (C1) that

 εh(λ)→0 (h→0),   ε∗h(λ∗)→0 (h→0).

Suppose that , then from (C2) we get

 εh(λ)≲hr,   ε∗h(λ∗)≲hr. \hb@xt@.01(3.15)

Further suppose that , is the Argyris finite element space, then from the interpolation theory we have

 εh(λ)≲h4,   ε∗h(λ∗)≲h4. \hb@xt@.01(3.16)

Note that when the functions in and are piecewise smooth (LABEL:s3.15) and (LABEL:s3.16) are also valid.
Thanks to [2], we get the following Theorem 3.3.

Theorem 3.3

Suppose , then

 ˆδ(R(E),R(Eh))≲εh(λ), \hb@xt@.01(3.17) |λ−(1qk+q−1∑j=kλ−1j,h)−1|≲εh(λ)ε∗h(λ∗), \hb@xt@.01(3.18) |λ−λj,h|≲[εh(λ)ε∗h(λ∗)]1α,   j=k,k