## Computation of the interior transmission eigenvalues for elastic scattering in an inhomogeneous medium containing an obstacle

### Wei-Chen Chang

^{∗}

### Tiexiang Li

^{†}

### Wen-Wei Lin

^{‡}

### Jenn-Nan Wang

^{§}

Abstract

In this work, we study the interior transmission eigenvalues for elas- tic scattering in an inhomogeneous medium containing an obstacle. This problem is related to the reconstruction of the support of the inhomo- geneity without the knowledge of the embedded obstacle by the far-field data or the invisibility cloaking of an obstacle. Our goal is to provide an efficient numerical algorithm to compute as many positive interior trans- mission eigenvalues as possible. We consider two cases of medium jumps:

Case 1, where C0 = C1, ρ0 6= ρ1, and Case 2, where C0 6= C1, ρ0 = ρ1

with either Dirichlet or Neumann boundary conditions on the boundary of the embedded obstacle. The partial differential equation (PDE) problem is reduced to a generalized eigenvalue problem (GEP) for matrices by the finite element method (FEM). We will apply the Jacobi-Davidson (JD) algorithm to solve the GEP. Case 1 requires special attention because of the large number of zero eigenvalues, which depends on the discretization size. To compute the positive eigenvalues effectively, it is necessary to deflate the zeros to infinity at the beginning of the algorithm.

Keywords — Interior transmission eigenvalues, elastic waves, Jacobi-Davidson method, nonequivalence deflation.

### 1 Introduction

In this paper, we study the interior transmission eigenvalue problem (ITEP) for elastic waves propagating outside of an obstacle. Our aim is to design an

∗Department of Mathematics, National Taiwan University, Taipei 106, Taiwan. Email:

d04221003@ntu.edu.tw

†School of Mathematics, Southeast University, Nanjing 211189 & Nanjing Center for Ap- plied Mathematics, Nanjing 211135, China. Email:txli@seu.edu.cn

‡Department of Applied Mathematics, Yang Ming Chiao Tung University, Hsinchu 300, Taiwan & Nanjing Center for Applied Mathematics, Nanjing 211135, China. E-mail:

wwlin@math.nctu.edu.tw

§Institute of Applied Mathematical Sciences, National Taiwan University, Taipei, 106, Taiwan. Email: jnwang@math.ntu.edu.tw. On behalf of all authors, the corresponding author states that there is no conflict of interest.

efficient numerical algorithm to compute as many transmission eigenvalues as
possible for time-harmonic elastic waves. Let D and Ω be open bounded domains
in R^{2} with smooth boundaries ∂D and ∂Ω, respectively. Assume that D ⊂
Ω. Let u(x) = [u1(x), u2(x)]^{>} be a two-dimensional vector representing the
displacement vector, and let its infinitesimal strain tensor be given by ε(u) =
((∇u)^{T} + ∇u)/2. We consider the linear elasticity; that is, the stress tensor
σ_{C}(u) is defined by ε(u) via Hook’s law:

σC(u) = Cε(u),

where C is the elasticity tensor. The elasticity tensor C = (Cijkl), 1 ≤ i, j, k, l ≤ 2, is a fourth-rank tensor satisfying two symmetry properties:

Cijkl= Cklij (major symmetry),

Cijkl= Cjikl (minor symmetry). (1) We require that C satisfies the strong convexity condition: there exists a κ > 0 such that for any symmetric matrix A

CA : A ≥ κ|A|^{2}, ∀ x ∈ Ω. (2)

In the following, for any two matrices A, B, we denote A : B =P

ijaijbij and

|A|^{2}= A : A. The elastic body is called isotropic if
C_{ijkl}= λδ_{ij}δ_{kl}+ µ(δ_{ik}δ_{jl}+ δ_{il}δ_{jk}),

where µ and λ are called Lam´e coefficients. In other words, for an isotropic elastic body, the stress-strain relation is given by

σC(u) = 2µε(u) + λtr(ε(u))I = 2µε(u) + λdivuI,

where I stands for the identity matrix. It is not difficult to check that the convexity condition (2) is equivalent to

µ(x) > 0, λ(x) + µ(x) > 0. (3)
The core of the ITEP is to find ω^{2} ∈ C such that there exists a nontrivial
solution (u, v) ∈ [H^{1}(Ω)]^{2}× [H^{1}(Ω \ D)]^{2}of

∇ · σC_{0}(u) + ρ0ω^{2}u = 0 in Ω, (4a)

∇ · σC1(v) + ρ_{1}ω^{2}v = 0 in Ω \ D, (4b)

Bv = 0 on ∂D, (4c)

u = v on ∂Ω, (4d)

σC_{0}(u)ν = σC_{1}(v)ν on ∂Ω, (4e)

where C_{0}, C_{1} are elasticity tensors, ρ_{0}, ρ_{1} are density functions, and ν is the
outer normal of ∂Ω. Here, we define the boundary operator B on ∂D

Bv = v|∂D or Bv = σC_{1}(v)ν|∂D, (5)

where, to abuse the notation, ν denotes the unit normal on ∂D pointing into
the interior of D. Recall that σC_{1}(v)ν represents the traction acting on ∂D or

∂Ω.

The investigation of the ITEP (4) is motivated by the following practical
problem. Let us assume that ρ0 is a positive constant and supp(C0− C1) ⊂
Ω \ D, supp(ρ_{1}− ρ0) ⊂ Ω \ D. We can regard (C_{1}− C0, ρ_{1}− ρ0) as a ”coated”

material near ∂D. We now take any solution ˆu of

∇ · σ_{C}_{0}(ˆu) + ρ_{0}ω^{2}u = 0ˆ in R^{2}

as an incident field and consider the incident field ˆu scattered by the object D
and the inhomogeneity of the material, that is, (C1, ρ1). In particular, when
ρ0= 1 and C0is isotropic with constant Lam´e coefficients λ and µ, the typical
incident field ˆu =: ˆu^{in}_{e}(x) with e = p or s is given by

ˆ

u^{in}_{p}(x) = ξe^{ik}^{p}^{x·ξ} or uˆ^{in}_{s}(x) = ξ^{⊥}e^{ik}^{s}^{x·ξ}
where ξ ∈ Ω := {kξk_{2}= 1} and k_{p} := ω/√

λ + 2µ and k_{s}:= ω/√

µ represent the compressional and shear wave numbers, respectively. Let v be the total field satisfy

∇ · σ_{C}_{1}(v) + ρ_{1}ω^{2}v = 0 in R^{2}\ D.

Then, if ω^{2} is an interior transmission eigenvalue of (4) and u|Ω = ˆu|Ω, then
the obstacle and the inhomogeneity (C1− C0, ρ1− ρ0) would be nonscattered
objects at ω^{2}when the incident field is ˆu. Consequently, if our aim is to detect
D and the inhomogeneity (C1− C0, ρ1− ρ0) by the scattering information, we
have to avoid the interior transmission eigenvalues.

On the other hand, a more interesting implication of this investigation is to

”cloak” the domain D from the elastic waves with suitable coated materials.

Now, we regard the incident field ˆu as a source of seismic waves, i.e., p-waves
or s-waves. If ω^{2} is an interior transmission eigenvalue of (4), then D with
coated material C1, ρ1 will be ”invisible” from the seismic waves propagating
at frequency ω.

The study of the ITEP originates from the validity of some qualitative ap- proaches to the inverse scattering problems in an inhomogeneous medium, such as the linear sampling method [12] and the factorization method [18]. In recent years, the ITEP has attracted much attention in the study of direct/inverse scattering problems for acoustic and electromagnetic waves in inhomogeneous media [3,4,5,6,13,14,15,19,23]. For the investigation of the ITEP for elastic waves in the case of D = ∅, there are a few theoretical results [2, 1, 9, 11,10].

Recently, theoretical results on the discreteness and the existence of interior transmission eigenvalues of (4) were proven in [7].

The purpose of this paper is to develop a numerical method to compute the transmission eigenvalues of elastic waves (4). To put this work in perspective, we mention only some results related to the ITEP for elastic waves. As far as we know, two works have studied the computation of the ITEP for elastic waves in the case of D = ∅, see [8] and [16]. In [16], a numerical method

was presented to compute a few smallest positive transmission eigenvalues of
(4). The ITEP was reformulated as locating the roots of a nonlinear function
whose values are generalized eigenvalues of a series of self-adjoint fourth-order
problems. After discretizing the fourth-order eigenvalue problems using H^{2}-
conforming finite elements, a secant-type method was employed to compute the
roots of the nonlinear function. In [8], a numerical method based on the ideas
in [21,22] was proposed to compute many interior transmission eigenvalues for
the elastic waves. In this paper, not only was the case of different densities
considered but also the case of different elasticity tensors. The strategy used in
[21, 22,8] is as follows. One first discretizes (4) by the finite element method
(FEM). Then, the ITEP is transformed into a generalized eigenvalue problem
(GEP). By some ingenious observations, this GEP can be reduced to a quadratic
eigenvalue problem (QEP) and, at the same time, unwanted eigenvalues (0 or

∞ eigenvalues) can be removed. One then applies a quadratic Jacobi-Davidson (JD) method with nonequivalence deflation to compute the eigenvalues of the resulting QEP. Contrary to the results obtained in [16], the method implemented in [8] is able to locate a large number of positive transmission eigenvalues of (4).

The paper is organized as follows. In Section2, we describe the discretization of the ITEP using the FEM. The discretization reduces the PDE problem to a GEP. We will apply the JD method to locate positive eigenvalues of the GEP.

However, the existence of zero eigenvalues of the GEP will hinder our task.

Thus, in Section3, we discuss a nonequivalence deflation technique to remove the zero eigenvalues. In Section4, we describe our numerical algorithm based on the JD method in detail. Section5 contains all numerical simulations and remarks. We conclude the paper with a summary in Section6.

### 2 Discretization and the GEP

We first review the discretization of the ITEP (4) based on the standard piece- wise linear FEM (see [14] for details). Let

Sh= The space of continuous piecewise linear functions on Ω;

S_{h}^{D}= The subspace of functions in Sh that have vanishing DOF on Ω \ D;

S_{h}^{I} = The subspace of functions in S_{h} that have vanishing DOF on ¯D ∪ ∂Ω;

S_{h}^{Σ}= The subspace of functions in S_{h}that have vanishing DOF in D ∪ (Ω \ ¯D) ∪ ∂Ω;

S_{h}^{B} = The subspace of functions in S_{h}that have vanishing DOF in Ω,

where DOF is the degrees of freedom. Let {ψ_{i}}^{n}_{i=1}^{D}, {φ_{i}}^{n}_{i=1}^{I} , {θ_{i}}^{m}_{i=1}^{Σ}, and
{ξi}^{m}_{i=1}^{B}, denote standard nodal bases for the finite element spaces of S_{h}^{D}, S_{h}^{I},

S_{h}^{Σ}, and S_{h}^{B}, respectively. We then set

Ψ1, Ψ2, · · · , Ψ2n_{D} =ψ1 ψ2 · · · ψn_{D}

ψ1 ψ2 · · · ψn_{D}

,

Φ1, Φ2, · · · , Φ2n_{I} =φ1 φ2 · · · φn_{I}

φ1 φ2 · · · φn_{I}

,

Ξ_{1}, Ξ2, · · · , Ξ2m_{B} =ξ1 ξ2 · · · ξm_{B}

ξ1 ξ2 · · · ξm_{B}

,

Θ_{1}, Θ_{2}, · · · , Θ_{2m}_{Σ} =θ_{1} θ2 · · · θmΣ

θ_{1} θ_{2} · · · θ_{m}_{Σ}

, and

u = u^{D}_{h} + u^{I}_{h}+ u^{B}_{h} + u^{Σ}_{h} =

2n_{D}

X

j=1

u^{D}_{j} Ψj+

2n_{I}

X

j=1

u^{I}_{j}Φj+

2m_{B}

X

j=1

wjΞj+

2m_{Σ}

X

j=1

u^{Σ}_{j}Θj

v = v_{h}^{I}+ v_{h}^{B}+ v_{h}^{Σ} =

2nI

X

j=1

v^{I}_{j}Φ_{j}+

2mB

X

j=1

w_{j}Ξ_{j}+

2mΣ

X

j=1

v_{j}^{Σ}Θ_{j}

Here we take into account the boundary condition u = v on ∂Ω and set u^{B}_{h} =
v^{B}_{h}. Expressed by the nodal bases, u and v have different dimensions. We will
discuss FEM for the Dirichlet and Neumann data separately.

### 2.1 Dirichlet condition: v = 0 on ∂D

For the zero Dirichlet condition, we take v^{Σ}_{j} = 0. Applying the standard piece-
wise linear FEM and using the integration by parts, we obtain

2nI

X

j=1

u^{I}_{j}(σ_{C}_{0}(Φ_{j}), ∇Φ_{i})_{Ω}+

2mB

X

j=1

w_{j}(σ_{C}_{0}(Ξ_{j}), ∇Φ_{i})_{Ω}+

2mΣ

X

j=1

u^{Σ}_{j} (σ_{C}_{0}(Θ_{j}), ∇Φ_{i})_{Ω}

= λ

2n_{I}

X

j=1

u^{I}_{j}(ρ0Φj, Φi)_{Ω}+

2m_{B}

X

j=1

wj(ρ0Ξj, Φi)_{Ω}+

2m_{Σ}

X

j=1

u^{Σ}_{j} (ρ0Θj, Φi)_{Ω}

, i = 1, · · · , 2nI,

2n_{D}

X

j=1

u^{D}_{j} (σ_{C}_{0}(Ψ_{j}), ∇Ψ_{i})_{Ω}+

2m_{Σ}

X

j=1

u^{Σ}_{j} (σ_{C}_{0}(Θ_{j}), ∇Ψ_{i})_{Ω}

= λ

2n_{D}

X

j=1

u^{D}_{j} (ρ0Ψj, Ψi)_{Ω}+

2m_{Σ}

X

j=1

u^{Σ}_{j} (ρ0Θj, Ψi)_{Ω}

, i = 1, · · · , 2nD,

2n_{I}

X

j=1

u^{I}_{j}(σ_{C}_{0}(Φ_{j}), ∇Ξ_{i})_{Ω}+

2m_{B}

X

j=1

w_{j}(σ_{C}_{0}(Ξ_{j}), ∇Ξ_{i})_{Ω}− (σ_{C}_{0}(u)ν, Ξ_{i})_{∂Ω}

= λ

2n_{I}

X

j=1

u^{I}_{j}(ρ0Φj, Ξi)_{Ω}+

2m_{B}

X

j=1

wj(ρ0Ξj, Ξi)_{Ω}

, i = 1, · · · , 2mB, (6) and

2n_{I}

X

j=1

u^{I}_{j}(σC_{0}(Φj), ∇Θi)_{Ω}+

2n_{D}

X

j=1

u^{D}_{j} (σC_{0}(Ψj), ∇Θi)_{Ω}+

2m_{Σ}

X

j=1

u^{Σ}_{j} (σC_{0}(Θj), ∇Θi)_{Ω}

= λ

2n_{I}

X

j=1

u^{I}_{j}(ρ0Φj, Θi)_{Ω}+

2n_{D}

X

j=1

u^{D}_{j} (ρ0Ψj, Θi)_{Ω}+

2m_{Σ}

X

j=1

u^{Σ}_{j} (ρ0Θj, Θi)_{Ω}

, i = 1, · · · , 2mΣ,

where λ = ω^{2}. Likewise for v but with v^{Σ}_{j} = 0, we have

2n_{I}

X

j=1

v^{I}_{j}(σ_{C}_{1}(Φ_{j}), ∇Φ_{i})_{D}c+

2m_{B}

X

j=1

w_{j}(σ_{C}_{1}(Ξ_{j}), ∇Φ_{i})_{D}c

= λ

2n_{I}

X

j=1

v_{j}^{I}(ρ1Φj, Φi)_{D}c+

2m_{B}

X

j=1

wj(ρ1Ξj, Φi)_{D}c

, i = 1, · · · , 2nI, (7) and

2n_{I}

X

j=1

v^{I}_{j}(σC_{1}(Φj), ∇Ξi)_{Ω}+

2m_{B}

X

j=1

wj(σC_{1}(Ξj), ∇Ξi)_{Ω}− (σC_{1}(v)ν, Ξi)_{∂Ω}

= λ

2n_{I}

X

j=1

v_{j}^{I}(ρ1Φj, Ξi)_{Ω}+

2m_{B}

X

j=1

wj(ρ1Ξj, Ξi)_{Ω}

, i = 1, · · · , 2mB, (8)

where ν is the unit outer normal of ∂Ω and D^{c} = Ω \ ¯D. Finally, taking
into account boundary conditions (4d), (4e), applying the linear FEM to the
difference equation between u and v in D^{c}, and performing the integration by
parts again, using (6)-(8), we can derive for i = 1, · · · , mB,

2n_{I}

X

j=1

u^{I}_{j}(σC_{0}(Φj), ∇Ξi)_{Ω}− v^{I}_{j}(σC_{1}(Φj), ∇Ξi)_{Ω}

+

2m_{B}

X

j=1

wj (σC_{0}(Ξj), ∇Ξi)_{Ω}− (σC_{1}(Ξj), ∇Ξi)_{Ω}

= λ

2n_{I}

X

j=1

u^{I}_{j}(ρ0Φj, Ξi)_{Ω}− v_{j}^{I}(ρ1Φj, Ξi)_{Ω} +

2m_{B}

X

j=1

wj (ρ0Ξj, Ξi)_{Ω}− (ρ1Ξj, Ξi)_{Ω}

, i = 1, · · · , 2mB. (9) For clarity, we define the stiffness matrices and mass matrices as in Table1.

Additionally, we set u^{I} = u^{I}_{1}, · · · , u^{I}_{2n}_{I}>

, u^{D} = u^{D}_{1}, · · · , u^{D}_{2n}_{D}>

, w =

w1, · · · , w2m_{B}>

, u^{Σ} = u^{Σ}_{1}, · · · , u^{Σ}_{2m}_{Σ}>

, and v^{I} = v^{I}_{1}, · · · , v^{I}_{2n}_{I}>

. Then, the discretization gives rise to a generalized eigenvalue problem (GEP)

Kz = λMz, (10)

K =

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB}

0 −K_{1}^{I} 0 0 K_{1}^{IB}

0 0 K_{0}^{D} K_{0}^{DΣ} 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0
(K_{0}^{IB})^{>} (K_{1}^{IB})^{>} 0 0 K_{0}^{B}− K_{1}^{B}

, (11)

Matrix Dimension Definition Stiffness Matrix

K_{`}^{I} 0 2nI× 2nI (K_{`}^{I})ij = (σC_{`}(Φj), ∇Φi)_{Ω}
K_{`}^{D} 0 2nD× 2nD (K_{`}^{D})ij = (σC_{`}(Ψj), ∇Ψi)_{Ω}
K_{`}^{B} 0 2mB× 2mB (K_{`}^{B})ij = (σC_{`}(Ξj), ∇Ξi)_{Ω}
K_{`}^{Σ} 0 2mΣ× 2mΣ (K_{`}^{Σ})ij = (σC_{`}(Θj), ∇Θi)_{Ω}

K_{`}^{IB} 2nI× 2mB (K_{`}^{IB})ij = (σC_{`}(Ξj), ∇Φi)_{Ω}
K_{`}^{IΣ} 2nI× 2mΣ (K_{`}^{IΣ})ij = (σC_{`}(Θj), ∇Φi)_{Ω}
K_{0}^{DΣ} 2n_{D}× 2m_{Σ} (K_{0}^{DΣ})_{ij} = (σ_{C}_{0}(Θ_{j}), ∇Ψ_{i})_{Ω}
Ke_{`}^{Σ} 2mΣ× 2mΣ ( eK_{`}^{Σ})ij = (σC_{`}(Θj), ∇Θi)_{D}c

Mass Matrix

M_{`}^{I} 0 2n_{I}× 2nI (M_{`}^{I})_{ij} = (ρ_{`}Φ_{j}, Φ_{i})_{Ω}
M_{`}^{D} 0 2nD× 2nD (M_{`}^{D})ij = (ρ`Ψj, Ψi)_{Ω}
M_{`}^{B} 0 2m_{B}× 2mB (M_{`}^{B})_{ij} = (ρ_{`}Ξ_{j}, Ξ_{i})_{Ω}
M_{`}^{Σ} 0 2mΣ× 2mΣ (M_{`}^{Σ})ij = (ρ`Θj, Θi)_{Ω}

M_{`}^{IB} 2nI× 2mB (M_{`}^{IB})ij = (ρ`Ξj, Φi)_{Ω}
M_{`}^{IΣ} 2nI× 2mΣ (M_{`}^{IΣ})ij = (ρ`Θj, Φi)_{Ω}
M_{0}^{DΣ} 2nD× 2mΣ (M_{0}^{DΣ})ij = (ρ0Θj, Ψi)_{Ω}
Mf_{`}^{Σ} 2mΣ× 2mΣ ( fM_{`}^{Σ})ij = (ρ`Θj, Θi)_{D}c

Table 1: Stiffness and mass matrices with ` = 0, 1

M =

M_{0}^{I} 0 0 M_{0}^{IΣ} M_{0}^{IB}

0 −M_{1}^{I} 0 0 M_{1}^{IB}

0 0 M_{0}^{D} M_{0}^{DΣ} 0

(M_{0}^{IΣ})^{>} 0 (M_{0}^{DΣ})^{>} M_{0}^{Σ} 0
(M_{0}^{IB})^{>} (M_{1}^{IB})^{>} 0 0 M_{0}^{B}− M_{1}^{B}

, z =

u^{I}

−v^{I}
u^{D}
u^{Σ}
w

.

### 2.2 Neumann condition: σ

_{C}

_{1}

### (v)ν = 0 on ∂D

We now consider the homogeneous Neumann condition. The equations for u remain the same. For v, using the integration by parts and the boundary condition σC1(v)ν = 0 on ∂D, we have

2n_{I}

X

j=1

v^{I}_{j}(σC_{1}(Φj), ∇Φi)_{Ω}+

2m_{B}

X

j=1

wj(σC_{1}(Ξj), ∇Φi)_{Ω}+

2m_{Σ}

X

j=1

v^{Σ}_{j} (σC_{1}(Θj), ∇Φi)_{D}c

= λ

2nI

X

j=1

v_{j}^{I}(ρ1Φj, Φi)_{Ω}+

2mB

X

j=1

wj(ρ1Ξj, Φi)_{Ω}+

2mΣ

X

j=1

v_{j}^{Σ}(ρ1Θj, Φi)_{D}c

,
i = 1, · · · , 2n_{I}, (12)
to replace the equation (7) and additionally we have

2n_{I}

X

j=1

v^{I}_{j}(σC_{1}(Φj), ∇Θi)_{Ω}+

2m_{Σ}

X

j=1

v_{j}^{Σ}(σC_{1}(Θj), ∇Θi)_{D}c

= λ

2n_{I}

X

j=1

v_{j}^{I}(ρ1Φj, Θi)_{Ω}+

2m_{Σ}

X

j=1

v_{j}^{Σ}(ρ1Θj, Θi)_{D}c

, i = 1, · · · , 2mΣ.

Similarly, taking into account the boundary conditions on ∂Ω, applying the
linear FEM to the difference equation between u and v, and performing the
integration by parts, we obtain equation (9). With v^{Σ} = v^{Σ}_{1}, · · · , v^{Σ}_{2m}

Σ

>

, expressing the system in the matrix form gives the following GEP

Kz = λMz, (13)

K =

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB} 0

0 −K_{1}^{I} 0 0 K_{1}^{IB} K_{1}^{IΣ}

0 0 K_{0}^{D} K_{0}^{DΣ} 0 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0 0
(K_{0}^{IB})^{>} (K_{1}^{IB})^{>} 0 0 K_{0}^{B}− K_{1}^{B} 0

0 (K_{1}^{IΣ})^{>} 0 0 0 Ke_{1}^{Σ}

, (14)

M =

M_{0}^{I} 0 0 M_{0}^{IΣ} M_{0}^{IB} 0

0 −M_{1}^{I} 0 0 M_{1}^{IB} M_{1}^{IΣ}

0 0 M_{0}^{D} M_{0}^{DΣ} 0 0

(M_{0}^{IΣ})^{>} 0 (M_{0}^{DΣ})^{>} M_{0}^{Σ} 0 0
(M_{0}^{IB})^{>} (M_{1}^{IB})^{>} 0 0 M_{0}^{B}− M_{1}^{B} 0

0 (M_{1}^{IB})^{>} 0 0 0 Mf_{1}^{Σ}

, z =

u^{I}

−v^{I}
u^{D}
u^{Σ}
w
v^{Σ}

.

### 3 Solving GEPs

In this section, we will attempt to solve GEPs (10) and (13). In each case of
the boundary condition on ∂D, we consider two situations where (i) C_{0} = C_{1}
and ρ_{0}6= ρ_{1}(Case 1) and where (ii) C_{0}6= C_{1}and ρ_{0}= ρ_{1}(Case 2). We discuss
the two cases separately.

### 3.1 Case 1 with Dirichlet condition

In this case, K_{0}^{B} = K_{1}^{B}, K_{0}^{I} = K_{1}^{I}, K_{0}^{IB} = K_{1}^{IB}, and the stiffness matrix K of
(10) becomes

K =

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB}

0 −K_{0}^{I} 0 0 K_{0}^{IB}

0 0 K_{0}^{D} K_{0}^{DΣ} 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0
(K_{0}^{IB})^{>} (K_{0}^{IB})^{>} 0 0 0

.

We would like to implement the Jacobi-Davidson (JD) method to solve (10). To make the numerical method effective, it is important to remove the null space of K as much as possible, which corresponds to zero eigenvalue (unphysical) of (10). To this end, we consider the linear system:

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB}

0 −K_{0}^{I} 0 0 K_{0}^{IB}

0 0 K_{0}^{D} K_{0}^{DΣ} 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0
(K_{0}^{IB})^{>} (K_{0}^{IB})^{>} 0 0 0

u1

−u1

u3

u4

u5

=

0 0 0 0 0

. (15)

The second equation of (15) reads

K_{0}^{I}u1+ K_{0}^{IB}u5= 0. (16)
By (16) and the first equation of (15), we have

K_{0}^{IΣ}u4= 0.

Since, in general, nI mΣ, it is reasonable to assume that K_{0}^{IΣ}is of full rank
and thus u4= 0. Using the third and fourth equations of (15), we immediately
obtain

u3= 0 and (K_{0}^{IΣ})^{>}u1= 0.

Combining this and (16) gives Au1

u5

:=

K_{0}^{I} K_{0}^{IB}
(K_{0}^{IΣ})^{>} 0

u1

u5

= 0. (17)

Note that the dimension of A is (2nI+ 2mΣ) × (2nI+ 2mB). So the nullity
of A is at most 2(mB− mΣ). Let [u^{>}_{1}, u^{>}_{5}] 6= 0 satisfy (17), then

u1

−u1

0 0 u5

(18)

is a null vector of K, i.e., a solution of (15).

Next, we need an important deflation technique to deflate eigenvalues (in- cluding zero eigenvalues) that have been computed to ∞. Suppose that we have computed

KX0= MX0Λ0, X0∈ R^{s×r}, Λ0∈ R^{r×r}, and s = 2(2nI+ nD+ mΣ+ mB),
that is, (X0, Λ0) is an eigenpair of (K, M). Let Y0be chosen so that Y_{0}^{>}MX0=
Ir. We then define

K = K − αMXe 0Y_{0}^{>}M,
M = M − MXf 0Y_{0}^{>}M,
where α 6∈ Spec(Λ0). Then we can show that
Theorem 1.

Spec( eK, fM) = (Spec(K, M) \ Spec(Λ0)) ∪ {∞}^{r}_{k=1},

where Spec( eK, fM) denotes the set of eigenvalues of the linear pencil eK − λ fM.

Proof. It suffices to compute

det( eK − λ fM) = det(K − λM − (α − λ)MX0Y_{0}^{>}M)

= det(K − λM)det(I_{s}− (α − λ)(K − λM)^{−1}MX0Y_{0}^{>}M). (19)
Recall that KX0 = MX0Λ0. We obtain (K − λM)X0 = MX0(Λ0− λIr) and
thus

(K − λM)^{−1}MX0= X0(Λ0− λIr)^{−1}.
Substituting this relation into (19) leads to

det( eK − λ fM) = det(K − λM)det(I_{s}− (α − λ)X_{0}(Λ_{0}− λI_{r})^{−1}Y_{0}^{>}M)

= det(K − λM)det(Λ0− λIr)^{−1}det(Λ0− λIr− (α − λ)Ir)

= det(K − λM)det(Λ0− λIr)^{−1}det(Λ0− αIr),

(20)
where we have used the identity det(Is+ AB) = det(Ir+ BA), where A ∈ R^{s×r},
B ∈ R^{r×s}, in the second equality above. The theorem follows easily from
(20).

In our numerical algorithm, we first use Theorem1to deflate the zero eigen- values to ∞ and compute a number of positive eigenvalues (from the smallest) of the associated matrix pair ( eK, fM). Since the zero eigenvalues have been de- flated, it will be quite effective to compute those small positive eigenvalues. We could continue deflating the computed eigenvalues to ∞. However, we will not do so due to the sparsity of matrices K, M. Note that the deflation process will destroy the sparsity of K, M. To keep the deflated matrices eK, fM as sparse as possible, after computing a number of positive eigenvalues, we first restore those eigenvalues that have been deflated to ∞ back to the matrix pair and deflate the positive eigenvalues that were just computed.

### 3.2 Case 1 with Neumann condition

Here, the stiffness matrix K of (13) becomes

K =

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB} 0

0 −K_{0}^{I} 0 0 K_{0}^{IB} K_{0}^{IΣ}

0 0 K_{0}^{D} K_{0}^{DΣ} 0 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0 0
(K_{0}^{IB})^{>} (K_{0}^{IB})^{>} 0 0 0 0
0 (K_{0}^{IΣ})^{>} 0 0 0 Ke_{0}^{Σ}

. (21)

We want to point out that eK_{1}^{Σ}here is evaluated in terms of the elasticity tensor
C0 since we have C1 = C0. Similarly, we want to find the null space of K as
much as possible. For this purpose, we want to consider the linear system

K_{0}^{I} 0 0 K_{0}^{IΣ} K_{0}^{IB} 0

0 −K_{0}^{I} 0 0 K_{0}^{IB} K_{0}^{IΣ}

0 0 K_{0}^{D} K_{0}^{DΣ} 0 0

(K_{0}^{IΣ})^{>} 0 (K_{0}^{DΣ})^{>} K_{0}^{Σ} 0 0
(K_{0}^{IB})^{>} (K_{0}^{IB})^{>} 0 0 0 0
0 (K_{0}^{IΣ})^{>} 0 0 0 Ke_{0}^{Σ}

u1

−u1

u3

u4

u5

u6

=

0 0 0 0 0 0

. (22)

By the first and second equations of (22), we have
K_{0}^{IΣ}(u4− u6) = 0.

Reasoning as above, we can assume that K_{0}^{IΣ}is of full rank and hence u4−u6=
0, i.e., u_{4}= u_{6}. It follows from the last equation of (22) that

(K_{0}^{IΣ})^{>}u_{1}= eK_{0}^{Σ}u_{4}. (23)
Substituting (23) into the fourth equation of (22) and combining its third equa-
tion gives

K_{0}^{D} K_{0}^{DΣ}
(K_{0}^{DΣ})^{>} (K_{0}^{Σ}+ eK_{0}^{Σ})

u3

u4

=0 0

. (24)

Since the matrix in (24) is symmetric and its diagonal matrices are positive-
definite, its rank is most likely small if it is singular. Therefore, we can take
u3 = u4 = 0 and (23) leads to (K_{0}^{IΣ})^{>}u1 = 0. Using the first or the second
equation again, we thus conclude that

Au_{1}
u5

:=

K_{0}^{I} K_{0}^{IB}
(K_{0}^{IΣ})^{>} 0

u_{1}
u5

= 0. (25)

Hence if [u^{>}_{1}, u^{>}_{5}] 6= 0 satisfies (25), then

u1

−u1

0
0
u_{5}

0

is a null vector of K in (21). Having discovered most of null vectors of K, we then solve the GEP (13) by the JD method with the deflation technique described as before.

### 3.3 Case 2 with Dirichlet condition

In this case, we have M_{1}^{I} = M_{0}^{I}, M_{1}^{IB} = M_{0}^{IB}, M_{1}^{B}= M_{0}^{B} and the mass matrix
M is reduced to

M =

M_{0}^{I} 0 0 M_{0}^{IΣ} M_{0}^{IB}

0 −M_{0}^{I} 0 0 M_{0}^{IB}

0 0 M_{0}^{D} M_{0}^{DΣ} 0

(M_{0}^{IΣ})^{>} 0 (M_{0}^{DΣ})^{>} M_{0}^{Σ} 0
(M_{0}^{IB})^{>} (M_{0}^{IB})^{>} 0 0 0

.

This reduced mass matrix has exactly the same form as K in (11). Arguing as above, wefindthe null space of M, which corresponds to ∞ eigenvalues of (10).

Since we are interested in locating first 500 positive eigenvalues, we do not need to carry out this step in the JD method.

### 3.4 Case 2 with Neumann condition

Similar to the case above, we have M_{1}^{I} = M_{0}^{I}, M_{1}^{IB} = M_{0}^{IB}, M_{1}^{B} = M_{0}^{B},
M_{1}^{IΣ}= M_{0}^{IΣ}, and fM_{1}^{Σ}= fM_{0}^{Σ}. Thus, the mass matrix becomes

M =

M_{0}^{I} 0 0 M_{0}^{IΣ} M_{0}^{IB} 0

0 −M_{0}^{I} 0 0 M_{0}^{IB} M_{0}^{IΣ}

0 0 M_{0}^{D} M_{0}^{DΣ} 0 0

(M_{0}^{IΣ})^{>} 0 (M_{0}^{DΣ})^{>} M_{0}^{Σ} 0 0
(M_{0}^{IB})^{>} (M_{0}^{IB})^{>} 0 0 0 0
0 (M_{0}^{IB})^{>} 0 0 0 Mf_{0}^{Σ}

;

while the stiffness matrix K is given by (14). This case can be treated similarly as in Subsection3.3.

### 4 Numerical strategies

We want to apply the JD method and the deflation technique to solve the GEP:

L(λ)z := (K − λM)z = 0. Roughly speaking, we will apply the JD method to the deflated pencil eL(λ) := eK − λ fM. In Case 1, we first deflate the zero eigenvalues to infinity and then apply the JD method to the deflated system to locate a small group of positive eigenvalues (15-20 positive eigenvalues). In Case 2, we begin with the direct implementation of the JD method to L(λ) and locate a small group of positive eigenvalues. Before continuing with the JD method, we first deflate these eigenvalues to infinity. However, we would not keep deflating found eigenvalues to infinity because that will destroy the sparsity of the system.

It is important to restore the previously deflated eigenvalues (including the zero eigenvalues in Case 1) before further deflation. In doing so, we always perturb L(λ) by matrices of lower ranks.

To make the paper self contained, we outline the JD method here (also see
[24]). Let V_{k}= [v_{1}, · · · , v_{k}] be a given orthogonal matrix and let (θ_{k}, u_{k}), θ_{k} 6=

0, be a Ritz pair (an approximate eigenpair) of eL(λ), i.e., uk= Vkskand (θk, sk)
is an eigenpair of V_{k}^{T}L(λ)Ve k, namely,

V_{k}^{T}L(θe k)Vksk = 0 (26)
with ks_{k}k = 1. Since V_{k}^{T}L(θe _{k})V_{k} is of lower rank, (26) can be solved by a usual
eigenvalue solver.

Starting from the Ritz pair (θk, uk), we aim to find a correction direction t⊥uk such that

L(λ)(ue k+ t) = 0. (27)

Let

rk= eL(θk)uk (28)

be the residual vector of eL(λ) corresponding to the Ritz pair (θk, uk). To solve t in (27), we rewrite and expand (27)

L(λ)t = − ee L(λ)uk= −r_{k}+ ( eL(θk) − eL(λ))uk

= −rk+ (λ − θk) fMuk.

(29) Using the fact that

u^{>}_{k}r_{k} = u^{>}_{k}L(θe k)u_{k} = s^{>}_{k}V_{k}^{>}L(θe k)V_{k}s_{k} = 0,
we multiply I −Muf _{k}u^{>}_{k}

u^{>}_{k}Muf k

!

on both sides of (29) to eliminate the term (λ−θ_{k})
and get

I −Muf ku^{>}_{k}
u^{>}_{k}Muf k

!

L(λ)t = −re k. (30)

Next, applying the orthogonal projection (I − uku^{>}_{k}) and approximating eL(λ)
by eL(θk), we then have the following correction equation:

I −Muf ku^{>}_{k}
u^{>}_{k}Muf k

!

L(θe _{k})(I − u_{k}u^{>}_{k})t = −r_{k} with t⊥u_{k}. (31)
In designing a stopping criterion in the JD method, we will check the resid-
ual of L(θ_{k})w_{k} because of the consideration of efficiency, where (θ_{k}, w_{k}) is a
Ritz pair of L(λ) related to (θ_{k}, u_{k}). Moreover, in the deflation process (see
Theorem1), we need to make use of the eigenpairs of the original pencil L(λ).

Therefore, it is required to transform the Ritz pair (θ_{k}, u_{k}) of eL(λ) to a Ritz
pair (θ_{k}, w_{k}) of L(λ). Recall that eL(λ) = eK − λ fM, where

K = K − αMXe 0Y_{0}^{>}M,
M = M − MXf _{0}Y_{0}^{>}M,
where

KX0= MX0Λ0, X0∈ R^{s×r}, Λ0∈ R^{r×r}, s = 2(2nI+ nD+ mΣ+ mB),
Y0 satisfies Y_{0}^{>}MX0= Ir, and α 6∈ Spec(Λ0).

Theorem 2. Let (θ, u), θ 6= 0 and θ 6∈ Spec(Λ0) ∪ {α}, be an eigenpair of eL(λ), i.e., eL(θ)u = 0. Then (θ, w) is an eigenpair of L(λ), where

w = u − X_{0}q (32)

andq = (α − θ)(Λ_{0}− θI_{r})^{−1}Y_{0}^{>}Mu.

Proof. Assume that (θ, u) satisfy

L(θ)u = ( ee K − θ fM)u = 0. (33)
Consider w = u − X_{0}q, where q ∈ R^{r×1}will be determined later. We now write
(33) as

(K − αMX0Y_{0}^{>}M)(w + X0q) = θ(M − MX0Y_{0}^{>}M)(w + X0q). (34)
In view of (34), Kw = θMw if and only if

KX0q − αMX0Y_{0}^{>}Mu = MX0Λ0q − αMX0Y_{0}^{>}Mu

= θMX0q − θMX0Y_{0}^{>}Mu,
i.e.,

MX0(Λ_{0}− θIr)q = (α − θ)MX_{0}Y_{0}^{>}Mu. (35)
Multiplying Y_{0}^{>}on both sides of (35) and using Y_{0}^{>}MX0= I_{r},one obtains that

q = (α − θ)(Λ0− θIr)^{−1}Y_{0}^{>}Mu.

By Theorem 2, if we have found a Ritz pair (θk, uk) of eL(λ) in the JD algorithm, we then transform uk into wk using (32) and check the stopping criterion to determine whether (θk, wk) is a good approximation of the eigen- pair for L(λ). We summarize our method in the following two algorithms. In Algorithm 1, we perform the deflation. Algorithm 2 describes the JD method (with partial deflation).

Algorithm 1 Deflated stiffness and mass matrices

Input: Stiffness matrix K, Mass matrix M, eigenvectors X1 with associated
eigenvalues Λ_{0}, and a α not in Λ_{0}

Output: Deflated stiffness matrix eK, deflated mass matrix fM, and an addi-
tional matrix Y_{0}

1: Compute Y_{0}= MX_{1};

2: Decompose Y_{0}= U DV^{>} by svd;

3: Update Y0= V D^{−1}U^{>};

4: Compute eK = K − αMX0Y_{0}^{>}M and fM = M − MX1Y_{0}^{>}M;

Algorithm 2 Jacobi-Davidson with partial deflation

Input: Stiffness matrix K, mass matrix M, the initial projected matrix V , expected eigenvalue number p, and deflated number l

Output: (λ_{j}, x_{j}) satisfying L(λ_{j})x_{j}≡ (K − λjM)xj= 0 for j = 1, · · · , p

1: if (Case 1) then

2: Use Algorithm 1to deflate zero eigenvalues and obtain matrices eK and M;f

3: Set X0= X1 and Θ = Λ0;

4: else

5: Set eK = K and fM = M;

6: Set X0= [ ] and Θ = [ ];

7: end if

8: for (j = 1, · · · , p) do

9: Compute W_{K}= eKV , W_{M}= fMV , H_{K}= V^{∗}W_{K}, and H_{M}= V^{∗}W_{M};

10: while (the desired eigenpair does not converge) do

11: Compute (θ, s) with ksk2= 1 such that (H_{K}− θH_{M})s = 0;

12: Compute u = V s, p = eL^{0}(θ)u, and r = eL(θ)u;

13: Solve a correction vector t ⊥ u from the correction equation;

14: Orthogonalize t against V and set v = t/ktk_{2};

15: Compute wK= eKv, wM= fMv and expand HK=

HK V^{∗}wK

v^{∗}WK v^{∗}wK

, HM=

HM V^{∗}wM

v^{∗}WM v^{∗}wM

;

16: Expand V =V v, WK=W_{K} w_{K}, and WM=W_{M} w_{M};

17: end while

18: Set λj = θ and zj= u;

19: Compute xj= zj− X0q where q = (α − λj)(Θ − λjIr)^{−1}Y_{0}^{>}Mzj;

20: if (j ≤ l) then

21: Set X0=X_{0} xj and Θ =Θ 0
0 λj

;

22: else

23: Set X0=X0(:, 2 : l) xj and Θ =Θ(2 : l, 2 : l) 0

0 λj

;

24: end if

25: Use Algorithm 1to obtain new deflated matrices eK and fM;

26: Update the initial matrix V ;

27: end for

Before presenting our numerical results in the next section, we would like to validate the GEP deduced from the FEM. In other words, the eigenvalues of the GEP are indeed approximations of the interior transmission eigenvalues.

To this end, we compare the theoretical computation for the radially symmetric transmission eigenfunctions obtained in [7] with the eigenvalues of the GEP.

For simplicity, we only consider the result of Case 1 with the zero Dirichlet
boundary condition. According to [7], Ω is a disk of radius 1 and the embedded
obstacle D is a disk of radius 0.5. The parameters are set to µ = λ = 1, ρ_{0}= 1
and ρ_{1} = 0.5. Note that the symmetry is in the sense of u = u^{2}_{1}+ u^{2}_{2} and
v = v_{1}^{2}+ v_{2}^{2}. The comparison is given in Table2. We also demonstrate that
the eigenfunctions computed from the GEP are indeed radially symmetric, as
shown in Figures1 and2.

ω Values obtained in [7] Values obtained from GEP

1 9.8599 9.8792

2 19.1916 19.3240

3 26.7594 27.1291

4 36.0246 36.8537

Table 2: Comparison of eigenvalues obtained in [7] and computed from the GEP for Case 1 with Dirichlet boundary condition.

Figure 1: Eigenfunctions u (left sub-figure) and v (right sub-figure) correspond-
ing to the eigenvalue ω = 9.8792. Here u and v are radially symmetric described
by u = u(r)e_{r}and v = v(r)e_{r}, where e_{r}is the unit vector directed at the radial
direction.