• 沒有找到結果。

An efficient numerical algorithm for computing densely distributed positive interior transmission eigenvalues

N/A
N/A
Protected

Academic year: 2022

Share "An efficient numerical algorithm for computing densely distributed positive interior transmission eigenvalues"

Copied!
23
0
0

加載中.... (立即查看全文)

全文

(1)

densely distributed positive interior transmission eigenvalues

Tiexiang Li

Department of Mathematics, Southeast University, Nanjing 211189, People’s Republic of China.

E-mail: txli@seu.edu.cn

Tsung-Ming Huang

Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan.

E-mail: min@ntnu.edu.tw

Wen-Wei Lin

Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan.

E-mail: wwlin@math.nctu.edu.tw

Jenn-Nan Wang

Institute of Applied Mathematical Sciences, NCTS, National Taiwan University, Taipei, 106, Taiwan.

E-mail: jnwang@math.ntu.edu.tw

Abstract. We propose an efficient eigensolver for computing densely distributed spectrum of the two-dimensional transmission eigenvalue problem (TEP) which is derived from Maxwell’s equations with Tellegen media and the transverse magnetic mode. The discretized governing equations by the standard piecewise linear finite element method give rise to a large-scale quadratic eigenvalue problem (QEP). Our numerical simulation shows that half of the positive eigenvalues of the QEP are densely distributed in some interval near the origin. The quadratic Jacobi-Davidson method with a so-called non-equivalence deflation technique is proposed to compute the dense spectrum of the QEP. Extensive numerical simulations show that our proposed method makes the convergence efficiently even it needs to compute more than 5000 desired eigenpairs. Numerical results also illustrate that the computed eigenvalue curves can be approximated by the nonlinear functions which can be applied to estimate the denseness of the eigenvalues for the TEP.

Keywords: Two-dimensional transmission eigenvalue problem, Tellegen model, transverse magnetic mode, quadratic Jacobi-Davidson method, non-equivalence deflation

(2)

1. Introduction

The transmission eigenvalue problems (TEP) have recently received a great deal of attention in the area of the inverse scattering which is essential in the study of direct/inverse scattering problems for nonabsorbing inhomogeneous media [3, 4, 5, 6, 8, 10, 11, 20, 27]. The existence of transmission eigenvalues is closely related to the validity of some reconstruction methods for the inverse scattering problems in an inhomogeneous medium such as the linear sampling method [1,7] and the factorization method [21]. On the other hand, it has been shown that transmission eigenvalues can be estimated by using the far-field data [2] or the near-field data (Cauchy data) [30].

In fact, both Dirichlet eigenvalues and transmission eigenvalues carry the information of the scatterer and can be estimated by either the far-field or the near-field data. For example, the Herglotz wave function [9] is applied at the first Dirichlet eigenvalue to reconstruct the shape of the sound-soft scatterer. The eigenvalue method using multiple frequency near-field data (EM2F) [32] is proposed to detect Dirichlet or transmission eigenvalues and to reconstruct the support of the scatterer by determining the indicator function. In addition, the EM2F can be used to distinguish between the sound-soft obstacle and nonabsorbing inhomogeneous medium. For further study in the theories and applications of TEP, we refer to [6] and the references therein for further details.

In this paper, we consider the scattering of acoustic waves by an inhomogeneous medium which occupies a bounded and simply connected domain D ⊆ R2. The related so-called TEP is to find λ > 0 and nontrivial functions u, v ∈ L2(D) with u − v ∈ H02(D) = { w ∈ H2| w = 0,∂w∂ν = 0 on ∂D} satisfying

∆u + λε(x, y)u = 0, in D, (1a)

∆v + λv = 0, in D, (1b)

where ν is the outer unit normal to the smooth boundary ∂D and ε(x, y) is the index of refraction. Any positive λ such that (1) has nontrivial solutions u and v is called an interior transmission eigenvalue. Note that u − v ∈ H02(D) is equivalent to

u = v on ∂D, (2a)

∂u

∂ν = ∂v

∂ν on ∂D. (2b)

Equation (1a) can be regarded as the reduced Maxwell equations with transverse magnetic (TM) modes. In this paper, we will consider the case where the parameter ε(x, y) are chosen sufficiently large. For the standard Maxwell equations, ε corresponds to the electric permittivity which is written as the product of the relative permittivity and the free space permittivity. In practice, ε can not be taken as large as we wish.

Therefore, in this paper, we will derive (1a) from Maxwell’s equations with Tellegen media (non-reciprocal and non-chiral). It turns out for this case ε is the sum of the electric permittivity and the square of Tellegen parameter. By choosing large Tellegen parameter, we can enlarge the parameter ε as well.

(3)

There are several theoretical results regarding the existence of interior transmission eigenvalues [4,5, 20,27]. Several numerical algorithms for computing the transmission eigenvalues were proposed recently. Three finite element methods (FEMs) and a coupled boundary element method were proposed for solving the two-dimensional (2D)/three-dimensional (3D) interior transmission eigenvalue problems in [10, 12, 33], and, furthermore, in a recently published book [34] for details. In view of the existence theory of the transmission eigenvalues based on the fourth order reformulation in [5], two iterative methods together with convergence analysis were studied in [31]. In [18], a mixed FEM for 2D TEP was proposed, which leads to a non-Hermitian quadratic eigenvaluue problem (QEP) and then was solved by an adaptive Arnoldi method. On the other hand, the multilevel correction method was used to reduce the solution of TEP into a series of solutions to some linear boundary value problems which could be solved by the multigrid method [19].

In many cases with general inhomogeneous medium, the desired positive transmission eigenvalues are surrounded by complex ones. Based on complex-valued contour integrals and the boundary integral equation, an accurate numerical method for computing interior transmission eigenvalues for many obstacles (different from spheres) in 3D acoustic scattering was presented in [22]. However, only constant index of refraction and smooth domains were treated there. The algorithm used in [22] was later extended to the electromagnetic case in [23,24]. The QEP above can be rewritten as a particular parametrized symmetric definite generalized eigenvalue problem (GEP).

For such GEP, the eigenvalue curves can be arranged in a monotonic order so that the desired curves are sequentially located by a new secant-type iteration (see [25] for 2D TEP and [13] for 3D TEP, respectively).

In this paper, we focus on the 2D TEP with complex media and make the following contributions.

• We derive the 2D TEP (1) with ε(x, y) = ˜ε(x, y) + γ2 from the Maxwell’s equation with non-reciprocal, non-chiral media and the transverse magnetic mode (TM).

Here ˜ε(x, y) is the electric permittivity and γ is the Tellegen parameter.

• Discretized (1) by the standard piecewise linear finite element method [10] results in a GEP. The GEP is then reduced to a QEP by deflating all nonphysical zeros.

Our numerical simulations indicate that half of the positive eigenvalues of the QEP are densely distributed in some interval.

• We adapt the quadratic Jacobi-Davidson (QJD) method with partial locking technique for computing a large number of desired eigenpairs of the QEP. In order to accelerate convergence, we also develop a so-called partial non-equivalence deflation technique combined with QJD to deflate the part of computed eigenvalues to infinity while keeping the other eigenvalues unchanged. Numerical results demonstrate that the new partial deflation technique makes the convergence efficiently for computing 5000 desired eigenpairs.

• Furthermore, we modify QJD with partial deflation technique so that it can be

(4)

applied to compute all the eigenvalues in a given interval. Therefore, we separate the densely distributed eigenvalues of the TEP into several subintervals and compute the desired eigenpairs simultaneously by the modified method. Numerical results show that this modified method can be applied to compute more than 10000 target eigenpairs in our model.

• According to the computed eigenvalues, we find that the logarithms of the eigenvalue curves vs log10(ε(x, y)) for the TEP (1) are approximated by straight lines with almost the same negative slope. Therefore, we can estimate the distribution of the eigenvalues for given ε(x, y) from these parallel straight lines.

This paper is organized as follows. Section 2 is devoted to derive a 2D TEP with TM mode in non-reciprocal, non-chiral media. A corresponding discretization TEP and its spectral analysis are given in Section 3. In Section 4, we develop a non-equivalence low-rank deflation which can be used to accelerate convergence of QJD for computing the desired positive eigenvalues. A practical QJD algorithm combined with non-equivalence deflation and numerical results are presented in Sections 5 and 6, respectively. Finally, a concluding remark is given in Section 7.

2. Description of the governing equation

We consider Maxwell’s equations in complex media [35]:

∇ × E = iω (µH + ζE) , (3a)

∇ × H = −iω (˜εE + ξH) , (3b)

where E and H are the electronic field and magnetic field, respectively, ω represents the frequency, ˜ε and µ are the permittivity and the permeability, respectively, ξ and ζ are 3-by-3 magnetoelectric parameter matrices in various forms for describing different types of materials (see [28, p.26] and [35, p.44]). In this paper, we assume µ = 1 and all other parameters are functions of x, y. For this situation, we study the TM mode for (3), i.e.,

E =h

0 0 E3(x, y)i>

, H =h

H1(x, y) H2(x, y) 0i>

. (4)

Let

ζ =

0 0 ζ1

0 0 ζ2

−ζ1 −ζ2 0

, ξ =

0 0 ξ1

0 0 ξ2

−ξ1 −ξ2 0

. (5)

Then, Eqs. (3a) implies that

yE3

−∂xE3 0

= iω

 H1 H2 0

+

 ζ1E3 ζ2E3 0

. (6)

(5)

Substituting (6) into (3b) yields

(iω)−1

0 0

−

2

∂x2 +∂y22

 E3

−

0 0

∂x2E3) −∂y1E3)

= − iω

ε˜

 0 0 E3

−

0 0 ξ1H1+ ξ2H2

, which implies that

− ∆E3

= ω2



˜

εE3− ξ1



(iω)−1

∂yE3− ζ1E3

 + ξ2



(iω)−1

∂xE3+ ζ2E3



+ iω ∂

∂x(ζ2E3) − ∂

∂y(ζ1E3)



= ω2(˜ε + ξ1ζ1+ ξ2ζ2) E3+ iω ∂

∂x(ζ2E3) − ∂

∂y (ζ1E3) + ξ1

∂yE3− ξ2

∂xE3

 . If we choose ζ1 = ξ1 = γ1 and ζ2 = ξ2 = γ2, where γ1, γ2 are real constants, then we have

−∆E3 = ω2 ˜ε + γ12+ γ22 E3 ≡ ω2ε(x, y)E3 (7) with

ε(x, y) = ˜ε(x, y) + γ12+ γ22. (8) Maxwell’s equations with real parameters ζ = ξ are called Tellegen model.

We now discuss the boundary conditions. Recall that ν is the outer unit normal to

∂D. Then it is readily seen that

E × ν =

 0 0 E3

×

 ν1 ν2 0

=

−ν2E3 ν1E3

0

and

(∇ × E) × ν =

∂yE3

∂x E3 0

×

 ν1 ν2 0

=

∂E3

∂y ν2+ ∂E∂x3ν1 0 0

=

∂E3

∂ν

0 0

.

In other words, boundary conditions E × ν|∂D and (∇ × E) × ν|∂D are equivalent to E3|∂D and ∂E3

∂ν |∂D.

We then arrive at the TEP for (1a)-(1b) and (2a)-(2b) with ε(x, y) given in (8).

(6)

3. Discretization of TEP and its spectral analysis

We briefly review the discretization of the TEP (1) based on the standard piecewise linear FEM (see [10] for details). Let

Sh = The space of continuous piecewise linear functions on D, ShI = The subspace of functions in Sh with vanishing DoF on ∂D, ShB = The subspace of functions in Sh with vanishing DoF in D.

Here DoF stands for degrees of freedom. Let {φi}ni=1 and {ψi}mi=1 denote standard nodal bases for the finite element spaces of ShI and ShB, respectively, then

u = uIh+ uBh =

n

X

i=1

uiφi+

m

X

i=1

wiψi, (9a)

v = vhI+ vhB =

n

X

i=1

viφi +

m

X

i=1

wiψi. (9b)

Applying the standard piecewise linear finite element method to (1a) and using integration by parts, we get

n

X

i=1

ui(∇φi, ∇φj) +

m

X

j=1

wi(∇ψi, ∇φj)

= ω2

n

X

i=1

ui(εφi, φj) +

m

X

i=1

wi(εψi, φj)

!

. (10)

Similarly, applying the standard piecewise linear finite element method to (1b), we have

n

X

i=1

vi(∇φi, ∇φj) +

m

X

j=1

wi(∇ψi, ∇φj) = ω2

n

X

i=1

vii, φj) +

m

X

i=1

wii, φj)

!

. (11)

Finally, applying the linear finite element method with boundary conditions (2a), (2b) and the integration by parts to the difference equation between (1a) and (1b), we obtain

n

X

i=1

(ui− vi)∇φi, ∇ψj

!

= ω2

n

X

i=1

ui(εφi, ψj) +

m

X

i=1

wi(εφi, ψj) −

n

X

i=1

vii, ψj) −

m

X

i=1

wii, ψj)

!

. (12) Hereafter, we define the stiffness matrices K, E, and mass matrices M1, Mε, F1, Fε, G1 and Gε as in Table 1. In addition, we set u = [u1, . . . , un]>, v = [v1, . . . , vn]>, and w = [w1, . . . , wm]>. Then, the discretizations of (10), (11) and (12) give rise to a generalized eigenvalue problem (GEP)

Az = λBz (13a)

(7)

stiffness matrix for interior meshes K = [(∇φi, ∇φj)]  0 ∈ Rn×n stiffness matrix for

E = [(∇φi, ∇ψj)] ∈ Rn×m interior/boundary meshes

mass matrices for interior M1 = [(φi, φj)]  0 ∈ Rn×n meshes with ε = 1 or ε > 1 Mε = [(εφi, φj)]  0 ∈ Rn×n mass matrices for interior/boundary F1 = [(φi, ψj)] ∈ Rn×m meshes with ε = 1 or ε > 1 Fε= [(εφi, ψj)] ∈ Rn×m mass matrices for boundary G1 = [(ψi, ψj)]  0 ∈ Rm×m meshes with ε = 1 or ε > 1 Gε = [(εψi, ψj)]  0 ∈ Rm×m

Table 1. Stiffness and mass matrices with ε(x, y) > 1 for (x, y) ∈ ¯D.

with λ = ω2,

A =

K 0 E

0 −K E

E> E> 0

, B =

Mε 0 Fε

0 −M1 F1 Fε> F1> Gε− G1

, z =

 u v w

. (13b) For simplicity, we denote

G = Gε− G1, M = Mε− M1, F = Fε− F1, (14a) Mc1 = M1− F1G−1F>, M = M − F Gc −1F>, K = K − EGb −1F>, (14b) and

S =h

K Ei

, T1 =h

M1 F1i

, M =

"

M F

F> G

#

. (14c)

Suppose that M  0 symmetric positive definite. Then it follows that G  0, M  0 and cM  0. The quadratic eigenvalue problem (QEP) is defined as

Q(λ)x ≡ λ2A2+ λA1 + A0 x = 0, (15) where A2, A1 and A0 are all n × n symmetric matrices given by

A2 = M1+ cM1Mc−1Mc1>+ F1G−1F1> (16a)

= M1+ T1M−1T1>,

A1 = −K − bK cM−1Mc1>− cM1Mc−1Kb>− EG−1F1>− F1G−1E> (16b)

= −K − SM−1T1>− T1M−1S>,

A0 = bK cM−1Kb>+ EG−1E> (16c)

= SM−1S>.

It has been shown [13] that the GEP (13) can be reduced to the QEP as in (15) and (16) with x = u − v in which all nonphysical zero are removed.

(8)

Theorem 1 ([13]). Let L(λ) and Q(λ) be matrix pencils defined as in (13) and (15), respectively. Then

σ(L(λ)) = σ(Q(λ)) ∪ {0, · · · , 0}

| {z }

m

.

Here, σ(·) denotes the spectrum of the associated matrix pencil.

Let (λ, x) be an eigenpair of (15), then

λ2(xA2x) + λ(xA1x) + (xA0x) = 0. (17) Suppose that A1 is symmetric negative definite and A2, A0  0. We then obtain that

b ≡ −xA1x > 0, a ≡ xA2x > 0, c ≡ xA0x > 0 which imply that the roots of the quadratic equation (17) are

λ+= b +√

b2− 4ac

2a > 0, λ= 2c b +√

b2− 4ac > 0 (18) provided that b2− 4ac > 0. Consequently, there are 2n positive eigenvalues of (15) and the associated eigenvectors are real.

Theorem 2. Let

W0 =

"

M F

F> G

#−1/2"

K E>

#

, W1 =

"

M F

F> G

#−1/2"

M1 F1>

#

, (19)

d0 = kW0k2, d1 = kW1k2, (20)

and

( a0 = λmin(A0), ¯a0 = λmax(A0) = d20,

a2 = λmin(A2) , ¯a2 = λmax(A2). (21) Suppose that

a1 = λmin(K) − 2d0d1 > 0, ¯a1 = λmax(K) + 2d0d1, (22) δ = λmin(K)2− 4d0(d1λmin(K) + d0λmax(M1)) > 0. (23) Then there are n positive eigenvalues of (15) in the interval (r, r), where

r = 2d20 a1+√

δ, r = 2λmin(A0)

¯ a1+√

¯

a1− 4a2a0 > 0. (24) Proof. By the definitions of W0 and W1, A1 in (16b) can be expressed by

A1 = −(K + W0>W1+ W1>W0). (25)

(9)

Let x be the unit eigenvector of (15). Note that x is real. Eq. (25) implies that b = −x>A1x = x>Kx + x>W0>W1x + x>W1>W0x

≥ λmin(K) − 2d0d1 = a1 > 0.

Recall that a = x>A2x and c = x>A0x. Then we have that b2− 4ac ≥ (λmin(K) − 2d0d1)2− 4d20min(M1) + d21)

= λmin(K)2− 4d0(d1λmin(K) + d0λmax(M1)) = δ > 0 and thus

λ = 2c

b +√

b2− 4ac ≤ 2d20 a1+√

δ = r. On the other hand, we can see that

b = −x>A1x ≤ λmax(K) + 2d0d1 = ¯a1, b2− 4ac ≤ ¯a21− 4a2a0,

which implies

λ = 2c

b +√

b2− 4ac ≥ 2a0

¯

a1+p¯a21− 4a2a0 = r.

It follows from (18) that there are n smallest positive eigenvalues on (r, r).

Remark 3. By (16c), the value x>A0x is dominated by x>K cbM−1Kb>x provided that λmin(K) ≡ O(κ)  1. From (16a), x>A2x ≈ O(1) holds. If we can choose the coefficient ε(x, y) in (1a) so that max |(M )ij| ≈ max |(K)ij|, then from (24) it follows that

λ ≈ 2O(κ)

O(κ) +pO(κ)2+ O(κ) ≈ O(1).

In other words, there are n positive eigenvalues of (1) which are densely distributed in the interval (0, O(1)). This motivates us to develop efficient numerical algorithms to compute all smallest clustering positive eigenvalues.

4. Non-equivalence low-rank deflation

In this section, we introduce the non-equivalence low-rank deflation [14] to locate the successive eigenpairs of the QEP in (15). Once the smallest positive eigenvalue is retrieved, we then transform it to infinity by the deflation scheme, while keeping all other eigenvalues unchanged. The next successive eigenvalue thus becomes the target positive eigenvalue of the new transformed problem, which is then repeatly solved by the proposed method.

(10)

Definition 4 ([14]). Let (Λ1, X1) ∈ Rr×r × Rn×r with rank(X1) = r be a given pair.

1, X1) is called an eigenmatrix pair of Q(λ) in (15) if it satisfies

A2X1Λ21+ A1X1Λ1+ A0X1 = 0. (26) In particular, (diag(∞, · · · , ∞), X1) is called an ∞-eigenmatrix pair for Q(λ) if A2X1 = 0.

Given an eigenmatrix pair (Λ1, X1) ∈ Rr×r× Rn×r(r ≤ n) of Q(λ) in (15) with Λ1 being nonsingular, we define a new deflated QEP as

Qd(λ) := λ2Ae2+ λ eA1+ eA0, (27) where

Ae2 := A2− A2X1Θ1X1>A2, (28a) Ae1 := A1+ A2X1Θ1Λ−>1 X1>A0+ A0X1Λ−11 Θ1X1>A2, (28b) Ae0 := A0− A0X1Λ−11 Θ1Λ−>1 X1>A0, (28c) and

Θ1 := (X1>A2X1)−1. (28d) The nonequivalence deflation (28) allows us to transform Q(λ) into a new QEP Qd(λ) with the same eigenvalues, except that the eigenvalues of Λ1 are replaced by r infinities.

On the other hand, let (Λ2, X2) ∈ Rs×s× Rn×s be another eigenmatrix pair of Q(λ).

Suppose that σ(Λ1) ∩ σ(Λ2) = ∅. Then the following orthogonality relation holds [14]

X2>A0X1− Λ>2(X2>A2X11 = 0. (29) Using this orthogonality relation, we can see that (Λ2, X2) is also an eigenmatrix pair of Qd(λ).

5. Jacobi-Davidson method for quadratic eigenvalue problems

In this section, we propose a quadratic Jacobi-Davidson (QJD) method [15, 29]

combined with non-equivalence deflation scheme to solve the QEP (15). Suppose Vk

is a k-dimensional subspace with an orthogonal unitary basis {v1, v2, . . . , vk}. Let Vk = [v1, · · · , vk] and (θk, sk) be an eigenpair of VkQ(λ)Vks = 0 and (θk, uk ≡ Vksk) be a Ritz pair of Q(λ) with kskk2 = 1. To expand the subspace Vk successively, the QJD method seaks an approximate solution for the correction equation:



I − pkuk ukpk



Q(θk) (I − ukuk) t = −rk, t⊥uk, (30) where rk= Q(θk)uk and pk = (2θkA2+ A1)uk. This is a crucial step in the QJD method that may affect the overall performance significantly.

(11)

There are different schemes [15] for solving (30). Here, based on the coefficient matrices A0, A1, A2 in (16), we adopt the scheme proposed in [15] to solve (30). Using the condition t⊥uk, Eq. (30) can be rewritten as

Q(θk)t = ukQ(θk)t

ukpk pk− rk ≡ ηkpk− rk. (31) We can then solve the two linear systems

Q(θk)t1 = pk, Q(θk)t2 = rk (32a) and compute the solution t of (31) as

t = ηkt1 − t2 with ηk= ukt2

ukt1. (32b)

Based on the above discussions, we summarize the quadratic Jacobi-Davidson method for the computation of the desired eigenvalue of the QEP (15) in Algorithm 1.

Algorithm 1 QJD method for Q(λ)x ≡ (λ2A2 + λA1+ A0)x = 0.

Input: Coefficient matrices A0, A1, A2 and an initial orthonormal matrix V . Output: The desired eigenpair (λ, x).

1: Compute Wi = AiV and Mi = VWi for i = 0, 1, 2.

2: while (the desired eigenpair is not convergent) do

3: Compute the eigenpairs (θ, s) of (θ2M2+ θM1+ M0)s = 0.

4: Select the desired eigenpair (θ, s) with ksk2 = 1.

5: Compute u = V s, p = (2θA2 + A1)u, r = Q(θ)u.

6: Solve the correction vector t in (32).

7: Orthogonalize t against V ; set v = t/ktk2.

8: Compute

wi = Aiv, Mi =

"

Mi Vwi vWi vwi

#

for i = 0, 1, 2.

9: Expand V = [V, v] and Wi = [Wi, wi] for i = 0, 1, 2.

10: end while

11: Set λ = θ and x = u.

Note that the solutions t1 and t2 in (32a) can be efficiently computed by the following way. Substituting A2, A1 and A0 in (16) into (32a), Eq. (32a) can be represented as

n

θk2M1− θkK +

θkMc1− bK

Mc−1

θkMc1>− bK>

+ (θkF1− E) G−1 θkF1>− E> t = b, (33)

(12)

where b = pk or b = rk. Let ˆt = cM−1

θkMc1>− bK>

t, ˜t = G−1 θkF1>− E> t, which is equivalent to



Kb>− θkMc1>

t + cMˆt = 0, E>− θkF1> t + G˜t = 0. (34a) Then (33) can be expressed by

θkK − θk2M1 t +

K − θb kMc1

ˆt + (E − θkF1) ˜t = −b. (34b) Combining (34a) with (34b), the solution t in (33) can be solved from the augmented linear system

Mc Kb>− θkMc1>

G E>− θkF1>

K − θb kMc1 E − θkF1 θkK − θk2M1

 ˆt

˜t t

=

 0 0

−b

. (35)

5.1. Partial locking scheme

To compute successively all other desired eigenvalues, deflation [14, 16] or locking [15, 16, 17, 26] scheme is necessary. The Jacobi-Davidson method incorporated with locking scheme is capable of calculating the smallest positive eigenvalue first and then of computing all other desired eigenvalues successively by suitably choosing the orthonormal searching space span(V ≡ [Vc, V0]), where the columns of Vc form an orthonormal basis for the subspace generated by the convergent eigenvectors and V0 is any matrix satisfying V0>V0 = I. Therefore, in each iteration of Algorithm 1, the convergent eigenvalues λ1, . . . , λj will be included in the set of eigenvalues of the projective QEP (θ2M2+ θM1+ M0)s = 0 in Line3 of Algorithm 1. So, the target Ritz value θ in Line 4of Algorithm 1is chosen as θ /∈ {λ1, . . . , λj}.

Let {λ1, . . . , λm} be the desired eigenvalues. If m is small, then the locking scheme can be applied to compute all desired eigenvalues successively. However, when m is large, locking all convergent eigenvectors in the searching subspace span(V ) will reduce the efficiency significantly because it increases the costs of computing eigenpairs of (θ2M2 + θM1 + M0)s = 0, the Ritz vectors u, and the reorthogonalization of the correction vector t against V in Lines 1,5 and 7, respectively, of Algorithm1. In order to remedy this drawback, we propose a partial locking scheme with at most locking ` convergent eigenvectors in each iteration. Namely, for the computation of the (j + 1)- th eigenpair (λj+1, xj+1) with j + 1 ≤ `, all the convergent eigenvectors x1, . . . , xj are locked in V which means that the columns of Vc is an orthonormal basis of the subspace span{x1, . . . , xj}. If j + 1 > `, then only the convergent eigenvectors xj+1−`, . . . , xj are locked. The algorithm is summarized in Algorithm 2.

(13)

Algorithm 2 Quadratic Jacobi-Davidson method with partial locking scheme.

Input: Coefficient matrices A0, A1, A2, number p of desired eigenvalues, locking number

` (` < p) and an initial orthonormal matrix V . Output: The desired eigenpair (λj, xj) for j = 1, . . . , p.

1: Set Vc = [ ];

2: for j = 1, . . . , p do

3: Use Algorithm 1with initial matrix V to compute the desired eigenpair (λj, xj);

4: if j ≤ ` then

5: Orthogonalize xj against Vc; set Vc= [Vc, xj/kxjk2];

6: else

7: Orthogonalize xj against Vc(:, 2 : `); set Vc= [Vc(:, 2 : `), xj/kxjk2];

8: end if

9: Find an initial matrix V0 such that V0>V0 = I with V ≡ [Vc, V0].

10: end for

5.2. Partial deflation scheme

In Section 4, an explicit nonequivalence low-rank deflation method is proposed to transform the convergent eigenvalues to the infinity, while all other eigenvalues remain unchanged. The next successive eigenvalue thus becomes the smallest positive one of the transformed problem. In this subsection, we will discuss how to efficiently apply QJD to solve the deflated QEP Qd(λ)x = 0 in (27).

Let Y0 = A0X1Λ−11 ∈ Rn×r and Y2 = A2X1 ∈ Rn×r. Then the matrices eA2, eA1 and Ae0 defined in (28c) can be written as

Ae2 = A2− Y2Θ1Y2>, (36a) Ae1 = A1+ Y2Θ1Y0>+ Y0Θ1Y2>, (36b) Ae0 = A0− Y0Θ1Y0>. (36c) As stated in (32a), to solve the correction vector td, we need to solve the linear system



θ2kAe2 + θkAe1+ eA0

t = b. (37)

In view of (36), (37) can be rewritten as

Q(θk) − (θkY2− Y0) Θ1 θkY2>− Y0> t = b. (38) Denote

U = θkY2− Y0.

Applying the Sherman-Morrison-Woodbury formula, the solution of (38) can be

(14)

computed as

t = Q(θk) − U Θ1U> −1

b

= Q(θk)−1b + Q(θk)−1U I − Θ1U>Q(θk)−1U−1

Θ1U>Q(θk)−1b

= Q(θk)−1b + Q(θk)−1U Θ−11 − U>Q(θk)−1U−1

U>Q(θk)−1b.

In other words, in each iteration, we need to solve the linear systems Q(θk)Z =h

˜

pk ˜rk Ui

, (39a)

where ˜rk = Qdk)˜uk and ˜pk = (2θkAe2+ eA1)˜uk, and compute the correction vector td as

td = ηk˜t1− ˜t2 with ηk= u˜k˜t2k˜t1

(39b) where

˜t1 = Z(:, 1) + Z(:, 3 : r + 2)(Θ−11 − U>Z(:, 3 : r + 2))−1U>Z(:, 1), (39c)

˜t2 = Z(:, 2) + Z(:, 3 : r + 2)(Θ−11 − U>Z(:, 3 : r + 2))−1U>Z(:, 2). (39d) From above, we can see that the computational cost of solving Qd(λ)x = 0 by QJD will be increasing as r increases. In order to reduce the computational cost, similar to the concept of partial locking scheme, we propose a partial deflation scheme with at most deflating ` convergent eigenvectors in each iteration. That is, for computing the (j + 1)- th eigenpair (λj+1, xj+1) with j + 1 ≤ `, all the convergent eigenvectors x1, . . . , xj are deflated. If j + 1 > `, then only the convergent eigenvectors xj+1−`, . . . , xj are deflated.

We summarize the steps in Algorithm 3.

6. Numerical results

In what follows, we will compare the efficiency of Algorithm 2 with ` = 20 and Algorithm 3 with ` = 10 for computing desired positive real transmission eigenvalues λi > 0, i = 1, . . . , p, on four different domains [25] as shown in Figure 1. The meshes for these domains are constructed by triangular meshes. The associated matrix dimensions n and m of matrices in Table 1are listed in Table 2. The distributions of {λ1, . . . , λp} with ε(x, y) = 50, 100, 500, 1000 and matrix dimensions n and m in Table 2 are shown in Figures 1(e), 1(f), 1(g) and 1(h), respectively. All the eigenvalues almost have an uniform distribution.

All computations for numerical test examples are carried out in MATLAB 2015b.

The linear system in (35) is solved by Gaussian elimination (i.e., left matrix divide in MATLAB). Tic and toc functions are used to measure elapsed time in second. For the hardware configuration, we use an HP server equipped with the RedHat Linux operating system, two Intel Quad-Core Xeon E5-2643 3.33 GHz CPUs and 96 GB of main memory.

(15)

Algorithm 3 Quadratic Jacobi-Davidson method with partial deflation scheme.

Input: Coefficient matrices A0, A1, A2, number p of desired eigenvalues, number ` (` < p) of deflation and an initial orthonormal matrix V .

Output: The desired eigenpair (λj, xj) for j = 1, . . . , p.

1: Set X1 = [ ], Y0 = [ ], Y2 = [ ] and Θ = [ ];

2: for j = 1, . . . , p do

3: Use Algorithm 1 with the initial matrix V and solve the correction vector td by (39) to compute the first desired eigenpair (λj, xj) of Qd(λ)x = 0;

4: Compute y0 = λ−1j A0xj and y2 = A2xj;

5: if j ≤ ` then

6: Set Θ =

"

Θ X1>y2

x>jY2 x>jy2

#

, Θ1 = Θ−1, X1 = [X1, xj], Y0 = [Y0, y0] and Y2 = [Y2, y2];

7: else

8: Set Θ =

"

Θ(2 : `, 2 : `) X1(:, 2 : `)>y2 x>j Y2(:, 2 : `) x>j y2

#

and Θ1 = Θ−1;

9: Set X1 = [X1(:, 2 : `), xj], Y0 = [Y0(:, 2 : `), y0] and Y2 = [Y2(:, 2 : `), y2];

10: end if

11: Update the initial matrix V .

12: end for

(a) Disk (b) Ellipse (c) Dumbbell (d) Peanut

102 103

ε(x, y) = ε0 0

50 100 150 200 250 300 350 400 450

eigenvalueλi

(e) λi for disk

102 103

ε(x, y) = ε0 0

50 100 150 200 250

eigenvalueλi

(f) λi for ellipse

102 103

ε(x, y) = ε0 0

100 200 300 400 500 600 700

eigenvalueλi

(g) λi for dumb- bell

102 103

ε(x, y) = ε0 0

100 200 300 400 500 600

eigenvalueλi

(h) λifor peanut

Figure 1. Model domains and the associated distributions of the 5000 positive target eigenvalues for ε(x, y) = 50, 100, 500, 1000.

Table 2. The matrix dimension n, m (K ∈ Rn×n, E ∈ Rn×m) of the matrices for the benchmark problems. Here the mesh size h ≈ 0.04.

Domain disk ellipse dumbbell peanut

(n, m) (124631, 1150) (71546, 976) (149051, 1871) (168548, 1492)

(16)

Table 3. The first and pth eigenvalues λ1 and λp computed by Algorithm3 and the associated (e10), ep0)) defined in (40), where p = 5000 and ε0= 5000.

Domain disk ellipse

(e10), ep0)) (0.00290, 4.25) (1.837 × 10−3, 2.129) (λ1, λp) (0.00294, 4.28) (2.087 × 10−3, 2.210)

(r1, rp) (0.014, 0.007) (0.119, 0.037)

Domain dumbbell peanut

(e10), ep0)) (0.0113, 6.42) (6.392 × 10−3, 5.621) (λ1, λp) (0.0114, 6.49) (6.440 × 10−3, 5.673)

(r1, rp) (0.009, 0.011) (0.007, 0.009)

6.1. Numerical validation for the clustering eigenvalues

In this section, we shall numerically validate that the TEP has densely distributed eigenvalues in the interval (0, O(1)) if the coefficient ε(x, y) in (1a) is sufficient large as shown in Remark3. Furthermore, we shall demonstrate that each eigencurve vs ε(x, y) can be numerically approximated by a nonlinear function.

In order to observe the variety of the distribution of the eigenvalues with various coefficients ε(x, y), we compute the fifty smallest positive real eigenvalues for each given constant ε(x, y) = ε0 and show the computed eigenvalues in Figures 2(a), 2(c), 2(e) and 2(g). From these results, we see that the distributions of these fifty eigenvalues are clustered to the interval (0.2,1) as ε0 approaches to 103.

On the other hand, these results also show that each eigencurve λi0) can be approximated by a nonlinear function

λi0) ≈ ei0) ≡ bi× ε−a0 i (40) with constant ai and bi for i = 1, . . . , 50. We show the nonlinear functions e1(ε) and e50(ε) in Figures 2(a), 2(c), 2(e) and 2(g) with red lines. These approximations can be extended to other eigencurves. Using the eigenvalues shown in Figures 1(e), 1(f), 1(g) and 1(h), we get the coefficients ai and bi for i = 1, . . . , 5000 as shown in Figures 2(b), 2(d), 2(f) and 2(h), respectively. The approximation in (40) can be used to estimate the eigenvalues for a given ε0. In Table3, we demonstrate the computed eigenvalues λ1 and λp by Algorithm 3 and (e10), ep0)) in (40) for ε0 = 5000 and p = 5000 with the domains in Figure 1. The results show that the relative errors ri ≡ |λi− ei|/|λi| can be achieved about 0.01 for i = 1 and p.

Furthermore, the curves of the coefficients ai and bi for i = 1, . . . , 5000 in Figures 2(b), 2(d), 2(f) and 2(h) can be approximated by a linear function

ai ≈ `a(i) ≡ α1× i + α0 (41a) and a nonlinear function

bi ≈ eb(i) ≡ β0× iβ1 × 10β2(log10(i))2, (41b)

(17)

ε

102 103

eigenvalues

10-2 10-1 100 101

e1(ε) = 15.16 × ε−1.0047

e50(ε) = 264.2 × ε−0.99995

(a) Eigenvalues for disk

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

i 0.5

1 1.5 2 2.5 3 3.5 4 4.5 5

values

eb(i) = 10.99 × i0.7615× 100.03604(log10(i))2

a(i) = 1.004 + 5.51 × 10−7× i ai

log10(bi) log10(eb(i))

(b) Coefficients ai and bi for disk

ε

102 103

eigenvalues

10-2 10-1 100 101

e1(ε) = 13.24 × ε−1.0031

e50(ε) = 141.4 × ε−1.0006

(c) Eigenvalues for ellipse

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

i 0.5

1 1.5 2 2.5 3 3.5 4 4.5 5

values eb(i) = 3.894 × i0.8401× 100.02652(log10(i))2

a(i) = 0.9962 + 2.707 × 10−6× i ai

log10(bi) log10(eb(i))

(d) Coefficients ai and bi for ellipse

ε

102 103

eigenvalues

10-2 10-1 100 101

e1(ε) = 59.03 × ε−1.0048

e50(ε) = 489.4 × ε−1.0038

(e) Eigenvalues for dumbbell

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

i 0.5

1 1.5 2 2.5 3 3.5 4 4.5 5

values

eb(i) = 33.3 × i0.5876× 100.06127(log10(i))2

a(i) = 1.005 + 3.838 × 10−7× i ai

log10(bi) log10(eb(i))

(f) Coefficients ai and bi for dumbbell

ε

102 103

eigenvalues

10-2 10-1 100 101

e1(ε) = 32.89 × ε−1.0031

e50(ε) = 368.8 × ε−1.0007

(g) Eigenvalues for peanut

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

i 0.5

1 1.5 2 2.5 3 3.5 4 4.5 5

values

eb(i) = 15.08 × i0.7587× 100.03568(log10(i))2

a(i) = 1.003 + 8.248 × 10−7× i ai

log10(bi) log10(eb(i))

(h) Coefficients ai and bi for peanut

Figure 2. Eigenvalues with various ε(x, y) and the coefficients ai and bi in the nonlinear functions.

參考文獻

相關文件

Results for such increasing stability phenomena in the inverse source problems for the acoustic, electromagnetic, and elastic waves can be found in [ABF02, BLT10, BHKY18, BLZ20,

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

Keywords: Interior transmission eigenvalues, elastic waves, generalized eigen- value problems, quadratic eigenvalue problems, quadratic Jacobi-Davidson method..

In the inverse boundary value problems of isotropic elasticity and complex conductivity, we derive estimates for the volume fraction of an inclusion whose physical parameters

In this paper, we propose a practical numerical method based on the LSM and the truncated SVD to reconstruct the support of the inhomogeneity in the acoustic equation with

Like the proximal point algorithm using D-function [5, 8], we under some mild assumptions es- tablish the global convergence of the algorithm expressed in terms of function values,

Optim. Humes, The symmetric eigenvalue complementarity problem, Math. Rohn, An algorithm for solving the absolute value equation, Eletron. Seeger and Torki, On eigenvalues induced by

Due to the important role played by σ -term in the analysis of second-order cone, before developing the sufficient conditions for the existence of local saddle points, we shall