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

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² 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|², ∀ x ∈ Ω. (2)

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

ijaijbij and

|A|²= 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 ω² ∈ C such that there exists a nontrivial solution (u, v) ∈ [H¹(Ω)]²× [H¹(Ω \ D)]²of

∇ · σC₀(u) + ρ0ω²u = 0 in Ω, (4a)

∇ · σC1(v) + ρ₁ω²v = 0 in Ω \ D, (4b)

Bv = 0 on ∂D, (4c)

u = v on ∂Ω, (4d)

σC₀(u)ν = σC₁(v)ν on ∂Ω, (4e)

where C₀, C₁ are elasticity tensors, ρ₀, ρ₁ are density functions, and ν is the outer normal of ∂Ω. Here, we define the boundary operator B on ∂D

Bv = v|∂D or Bv = σC₁(v)ν|∂D, (5)

where, to abuse the notation, ν denotes the unit normal on ∂D pointing into the interior of D. Recall that σC₁(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(ρ₁− ρ0) ⊂ Ω \ D. We can regard (C₁− C0, ρ₁− ρ0) as a ”coated”

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

∇ · σ_C0(ˆu) + ρ₀ω²u = 0ˆ in R²

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ⁱⁿ_e(x) with e = p or s is given by

ˆ

uⁱⁿ_p(x) = ξe^ikpx·ξ or uˆⁱⁿ_s(x) = ξ^⊥e^iksx·ξ where ξ ∈ Ω := {kξk₂= 1} and k_p := ω/√

λ + 2µ and k_s:= ω/√

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

∇ · σ_C1(v) + ρ₁ω²v = 0 in R²\ D.

Then, if ω² 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 ω²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 ω² 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²- 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_hthat have vanishing DOF in D ∪ (Ω \ ¯D) ∪ ∂Ω;

S_h^B = The subspace of functions in S_hthat have vanishing DOF in Ω,

where DOF is the degrees of freedom. Let {ψ_i}ⁿ_i=1^D, {φ_i}ⁿ_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

 ,

Ξ₁, Ξ2, · · · , Ξ2m_B =ξ1 ξ2 · · · ξm_B

ξ1 ξ2 · · · ξm_B

 ,

Θ₁, Θ₂, · · · , Θ_2mΣ =θ₁ θ2 · · · θmΣ

θ₁ θ₂ · · · θ_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(σ_C0(Φ_j), ∇Φ_i)_Ω+

2mB

X

j=1

w_j(σ_C0(Ξ_j), ∇Φ_i)_Ω+

2mΣ

X

j=1

u^Σ_j (σ_C0(Θ_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 (σ_C0(Ψ_j), ∇Ψ_i)_Ω+

2m_Σ

X

j=1

u^Σ_j (σ_C0(Θ_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(σ_C0(Φ_j), ∇Ξ_i)_Ω+

2m_B

X

j=1

w_j(σ_C0(Ξ_j), ∇Ξ_i)_Ω− (σ_C0(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₀(Φj), ∇Θi)_Ω+

2n_D

X

j=1

u^D_j (σC₀(Ψj), ∇Θi)_Ω+

2m_Σ

X

j=1

u^Σ_j (σC₀(Θ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 λ = ω². Likewise for v but with v^Σ_j = 0, we have

2n_I

X

j=1

v^I_j(σ_C1(Φ_j), ∇Φ_i)_Dc+

2m_B

X

j=1

w_j(σ_C1(Ξ_j), ∇Φ_i)_Dc

= λ





2n_I

X

j=1

v_j^I(ρ1Φj, Φi)_Dc+

2m_B

X

j=1

wj(ρ1Ξj, Φi)_Dc



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

2n_I

X

j=1

v^I_j(σC₁(Φj), ∇Ξi)_Ω+

2m_B

X

j=1

wj(σC₁(Ξj), ∇Ξi)_Ω− (σC₁(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₀(Φj), ∇Ξi)_Ω− v^I_j(σC₁(Φj), ∇Ξi)_Ω

+

2m_B

X

j=1

wj (σC₀(Ξj), ∇Ξi)_Ω− (σC₁(Ξ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₁, · · · , u^I_2nI>

, u^D = u^D₁, · · · , u^D_2nD>

, w =

w1, · · · , w2m_B>

, u^Σ = u^Σ₁, · · · , u^Σ_2mΣ>

, and v^I = v^I₁, · · · , v^I_2nI>

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

Kz = λMz, (10)

K =













K₀^I 0 0 K₀^IΣ K₀^IB

0 −K₁^I 0 0 K₁^IB

0 0 K₀^D K₀^DΣ 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 (K₀^IB)^> (K₁^IB)^> 0 0 K₀^B− K₁^B













, (11)

Matrix Dimension Definition Stiffness Matrix

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

K_{`</sub><sup>IB</sup> 2nI× 2mB (K<sub>`}^IB)ij = (σC_{`</sub>(Ξj), ∇Φi)<sub>Ω</sub> K<sub>`}^IΣ 2nI× 2mΣ (K_{`</sub><sup>IΣ</sup>)ij = (σC<sub>`}(Θj), ∇Φi)_Ω K₀^DΣ 2n_D× 2m_Σ (K₀^DΣ)_ij = (σ_C0(Θ_j), ∇Ψ_i)_Ω Ke_{`</sub><sup>Σ</sup> 2mΣ× 2mΣ ( eK<sub>`}^Σ)ij = (σC_`(Θj), ∇Θi)_Dc

Mass Matrix

M_{`</sub><sup>I</sup>  0 2n<sub>I</sub>× 2nI (M<sub>`}^I)_ij = (ρ_{`</sub>Φ<sub>j</sub>, Φ<sub>i</sub>)<sub>Ω</sub> M<sub>`}^D 0 2nD× 2nD (M<sub>`</sub><sup>D</sup>)ij = (ρ`Ψj, Ψi)_Ω M_{`</sub><sup>B</sup> 0 2m<sub>B</sub>× 2mB (M<sub>`}^B)_ij = (ρ_{`</sub>Ξ<sub>j</sub>, Ξ<sub>i</sub>)<sub>Ω</sub> M<sub>`}^Σ 0 2mΣ× 2mΣ (M<sub>`</sub><sup>Σ</sup>)ij = (ρ`Θj, Θi)_Ω

M_{`</sub><sup>IB</sup> 2nI× 2mB (M<sub>`}^IB)ij = (ρ`Ξj, Φi)<sub>Ω</sub> M<sub>`</sub>^IΣ 2nI× 2mΣ (M<sub>`</sub><sup>IΣ</sup>)ij = (ρ`Θj, Φi)_Ω M₀^DΣ 2nD× 2mΣ (M₀^DΣ)ij = (ρ0Θj, Ψi)_Ω Mf_{`</sub><sup>Σ</sup> 2mΣ× 2mΣ ( fM<sub>`}^Σ)ij = (ρ`Θj, Θi)_Dc

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

M =













M₀^I 0 0 M₀^IΣ M₀^IB

0 −M₁^I 0 0 M₁^IB

0 0 M₀^D M₀^DΣ 0

(M₀^IΣ)^> 0 (M₀^DΣ)^> M₀^Σ 0 (M₀^IB)^> (M₁^IB)^> 0 0 M₀^B− M₁^B











 , z =











 u^I

−v^I u^D u^Σ w











 .

2.2 Neumann condition: σ

(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₁(Φj), ∇Φi)_Ω+

2m_B

X

j=1

wj(σC₁(Ξj), ∇Φi)_Ω+

2m_Σ

X

j=1

v^Σ_j (σC₁(Θj), ∇Φi)_Dc

= λ





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)_Dc



, i = 1, · · · , 2n_I, (12) to replace the equation (7) and additionally we have

2n_I

X

j=1

v^I_j(σC₁(Φj), ∇Θi)_Ω+

2m_Σ

X

j=1

v_j^Σ(σC₁(Θj), ∇Θi)_Dc

= λ





2n_I

X

j=1

v_j^I(ρ1Φj, Θi)_Ω+

2m_Σ

X

j=1

v_j^Σ(ρ1Θj, Θi)_Dc



, 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^Σ₁, · · · , v^Σ_2m

Σ

>

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

Kz = λMz, (13)

K =

















K₀^I 0 0 K₀^IΣ K₀^IB 0

0 −K₁^I 0 0 K₁^IB K₁^IΣ

0 0 K₀^D K₀^DΣ 0 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 0 (K₀^IB)^> (K₁^IB)^> 0 0 K₀^B− K₁^B 0

0 (K₁^IΣ)^> 0 0 0 Ke₁^Σ

















, (14)

M =

















M₀^I 0 0 M₀^IΣ M₀^IB 0

0 −M₁^I 0 0 M₁^IB M₁^IΣ

0 0 M₀^D M₀^DΣ 0 0

(M₀^IΣ)^> 0 (M₀^DΣ)^> M₀^Σ 0 0 (M₀^IB)^> (M₁^IB)^> 0 0 M₀^B− M₁^B 0

0 (M₁^IB)^> 0 0 0 Mf₁^Σ















 , 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₀ = C₁ and ρ₀6= ρ₁(Case 1) and where (ii) C₀6= C₁and ρ₀= ρ₁(Case 2). We discuss the two cases separately.

3.1 Case 1 with Dirichlet condition

In this case, K₀^B = K₁^B, K₀^I = K₁^I, K₀^IB = K₁^IB, and the stiffness matrix K of (10) becomes

K =













K₀^I 0 0 K₀^IΣ K₀^IB

0 −K₀^I 0 0 K₀^IB

0 0 K₀^D K₀^DΣ 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 (K₀^IB)^> (K₀^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₀^I 0 0 K₀^IΣ K₀^IB

0 −K₀^I 0 0 K₀^IB

0 0 K₀^D K₀^DΣ 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 (K₀^IB)^> (K₀^IB)^> 0 0 0























 u1

−u1

u3

u4

u5













=











 0 0 0 0 0













. (15)

The second equation of (15) reads

K₀^Iu1+ K₀^IBu5= 0. (16) By (16) and the first equation of (15), we have

K₀^IΣu4= 0.

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

u3= 0 and (K₀^IΣ)^>u1= 0.

Combining this and (16) gives Au1

u5

 :=

 K₀^I K₀^IB (K₀^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^>₁, u^>₅] 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₀^>MX0= Ir. We then define

K = K − αMXe 0Y₀^>M, M = M − MXf 0Y₀^>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₀^>M)

= det(K − λM)det(I_s− (α − λ)(K − λM)⁻¹MX0Y₀^>M). (19) Recall that KX0 = MX0Λ0. We obtain (K − λM)X0 = MX0(Λ0− λIr) and thus

(K − λM)⁻¹MX0= X0(Λ0− λIr)⁻¹. Substituting this relation into (19) leads to

det( eK − λ fM) = det(K − λM)det(I_s− (α − λ)X₀(Λ₀− λI_r)⁻¹Y₀^>M)

= det(K − λM)det(Λ0− λIr)⁻¹det(Λ0− λIr− (α − λ)Ir)

= det(K − λM)det(Λ0− λIr)⁻¹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₀^I 0 0 K₀^IΣ K₀^IB 0

0 −K₀^I 0 0 K₀^IB K₀^IΣ

0 0 K₀^D K₀^DΣ 0 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 0 (K₀^IB)^> (K₀^IB)^> 0 0 0 0 0 (K₀^IΣ)^> 0 0 0 Ke₀^Σ

















. (21)

We want to point out that eK₁^Σ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₀^I 0 0 K₀^IΣ K₀^IB 0

0 −K₀^I 0 0 K₀^IB K₀^IΣ

0 0 K₀^D K₀^DΣ 0 0

(K₀^IΣ)^> 0 (K₀^DΣ)^> K₀^Σ 0 0 (K₀^IB)^> (K₀^IB)^> 0 0 0 0 0 (K₀^IΣ)^> 0 0 0 Ke₀^Σ































 u1

−u1

u3

u4

u5

u6

















=















 0 0 0 0 0 0

















. (22)

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

Reasoning as above, we can assume that K₀^IΣis of full rank and hence u4−u6= 0, i.e., u₄= u₆. It follows from the last equation of (22) that

(K₀^IΣ)^>u₁= eK₀^Σu₄. (23) Substituting (23) into the fourth equation of (22) and combining its third equa- tion gives

 K₀^D K₀^DΣ (K₀^DΣ)^> (K₀^Σ+ eK₀^Σ)

 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₀^IΣ)^>u1 = 0. Using the first or the second equation again, we thus conclude that

Au₁ u5

 :=

 K₀^I K₀^IB (K₀^IΣ)^> 0

 u₁ u5



= 0. (25)

Hence if [u^>₁, u^>₅] 6= 0 satisfies (25), then















 u1

−u1

0 0 u₅

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₁^I = M₀^I, M₁^IB = M₀^IB, M₁^B= M₀^B and the mass matrix M is reduced to

M =













M₀^I 0 0 M₀^IΣ M₀^IB

0 −M₀^I 0 0 M₀^IB

0 0 M₀^D M₀^DΣ 0

(M₀^IΣ)^> 0 (M₀^DΣ)^> M₀^Σ 0 (M₀^IB)^> (M₀^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₁^I = M₀^I, M₁^IB = M₀^IB, M₁^B = M₀^B, M₁^IΣ= M₀^IΣ, and fM₁^Σ= fM₀^Σ. Thus, the mass matrix becomes

M =

















M₀^I 0 0 M₀^IΣ M₀^IB 0

0 −M₀^I 0 0 M₀^IB M₀^IΣ

0 0 M₀^D M₀^DΣ 0 0

(M₀^IΣ)^> 0 (M₀^DΣ)^> M₀^Σ 0 0 (M₀^IB)^> (M₀^IB)^> 0 0 0 0 0 (M₀^IB)^> 0 0 0 Mf₀^Σ

















;

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₁, · · · , 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^TL(λ)Ve k, namely,

V_k^TL(θe k)Vksk = 0 (26) with ks_kk = 1. Since V_k^TL(θ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^>_kr_k = u^>_kL(θe k)u_k = s^>_kV_k^>L(θe k)V_ks_k = 0, we multiply I −Muf _ku^>_k

u^>_kMuf k

!

on both sides of (29) to eliminate the term (λ−θ_k) and get

I −Muf ku^>_k u^>_kMuf 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^>_kMuf k

!

L(θe _k)(I − u_ku^>_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₀^>M, M = M − MXf ₀Y₀^>M, where

KX0= MX0Λ0, X0∈ R^s×r, Λ0∈ R^r×r, s = 2(2nI+ nD+ mΣ+ mB), Y0 satisfies Y₀^>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₀q (32)

andq = (α − θ)(Λ₀− θI_r)⁻¹Y₀^>Mu.

Proof. Assume that (θ, u) satisfy

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

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

KX0q − αMX0Y₀^>Mu = MX0Λ0q − αMX0Y₀^>Mu

= θMX0q − θMX0Y₀^>Mu, i.e.,

MX0(Λ₀− θIr)q = (α − θ)MX₀Y₀^>Mu. (35) Multiplying Y₀^>on both sides of (35) and using Y₀^>MX0= I_r,one obtains that

q = (α − θ)(Λ0− θIr)⁻¹Y₀^>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 Λ₀, and a α not in Λ₀

Output: Deflated stiffness matrix eK, deflated mass matrix fM, and an addi- tional matrix Y₀

1: Compute Y₀= MX₁;

2: Decompose Y₀= U DV^> by svd;

3: Update Y0= V D⁻¹U^>;

4: Compute eK = K − αMX0Y₀^>M and fM = M − MX1Y₀^>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⁰(θ)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₂;

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)⁻¹Y₀^>Mzj;

20: if (j ≤ l) then

21: Set X0=X₀ 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, ρ₀= 1 and ρ₁ = 0.5. Note that the symmetry is in the sense of u = u²₁+ u²₂ and v = v₁²+ v₂². 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_rand v = v(r)e_r, where e_ris the unit vector directed at the radial direction.