Efficient Arnoldi-type algorithms for rational eigenvalue problems arising in fluid-solid systems

Efficient Arnoldi-type algorithms for rational eigenvalue problems

arising in fluid–solid systems

So-Hsiang Chou

a,⇑

, Tsung-Ming Huang

b

, Wei-Qiang Huang

c

, Wen-Wei Lin

d

a

Department of Mathematics and Statistics, Bowling Green State University, Bowling Green, OH 43403-0221, USA

b

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

c

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

d

Department of Mathematics, National Taiwan University, Taipei 106, Taiwan

a r t i c l e

i n f o

Article history:

Received 31 August 2010

Received in revised form 14 December 2010 Accepted 14 December 2010

Available online 21 December 2010 Keywords:

Fluid–structure interaction Finite elements

Rational eigenvalue problem Trimmed linearization Arnoldi algorithm

a b s t r a c t

We develop and analyze efficient methods for computing damped vibration modes of an acoustic fluid confined in a cavity with absorbing walls capable of dissipating acoustic energy. The discretization in terms of pressure nodal finite elements gives rise to a rational eigenvalue problem. Numerical evidence shows that there are no spurious eigenmodes for such discretization and also confirms that the discretization based on nodal pressures is much more efficient than that based on Raviart–Thomas finite elements for the displace-ment field. The trimmed linearization method is used to linearize the associated rational eigenvalue problem into a generalized eigenvalue problem (GEP) of the form Ax ¼ kBx. For solving the GEP we apply Arnoldi algorithm to two different types of single matrices B1_{A and AB}1_{. Numerical accuracy shows that the application of Arnoldi on AB}1_{is better} than that on B1_A.

 2010 Elsevier Inc. All rights reserved.

1. Introduction

Efficient and correct computation of the damped vibration modes generated by an inviscid, compressible, barotropic fluid in a cavity with absorbing walls is an important issue when for example one is interested in decreasing the level of noise in aircraft or cars. In general, one needs first a mathematical model consisted of partial differential equations with proper boundary and initial conditions. After this first phase of mathematical formulation, the next phase is to find efficient meth-ods to compute the modes. This phase involves correct discretization of the mathematical formulation and computation of large scale nonlinear eigenvalue problems, be it quadratic, cubic, or even rational. Choosing correct discretization schemes to avoid spurious modes and finding efficient methods to locate eigenvalues that lie in the interior of the spectrum are among important issues to deal with. In the mathematical formulation phase, we have interaction between the fluid and structure (cavity walls), and the displacement variable natural for the solid could be chosen for the fluid as well so that compatibility and equilibrium (cf.(2.3) and (2.7)below) through the fluid–solid interface can be satisfied automatically. A drawback lurk-ing behind the displacement formulation is the possible presence of nonphysical zero-frequency spurious circulation modes, if one is not careful in choosing the discretization scheme associated with the underlying partial differential system. For example discretization by standard finite elements or finite differences often exhibit such a phenomenon. Approaches cir-cumventing this drawback can be found in[2,8,11,12,29], among others.

0021-9991/$ - see front matter  2010 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2010.12.022

⇑Corresponding author. Tel.: +1 419 372 8225; fax: +1 419 372 6092.

E-mail addresses:[email protected](S.-H. Chou),[email protected](T.-M. Huang),[email protected](W.-Q. Huang),wwlin@math. ntu.edu.tw(W.-W. Lin).

Contents lists available atScienceDirect

Journal of Computational Physics

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / j c p

One of the discretizations we will be using in this paper is the edge-based or Raviart–Thomas finite elements for the dis-placement field, following[3,5]. The main concerns in[3,19]are pure mathematical issues of proving that their numerical approximation is free of spurious modes and has second order convergence rate. Efficient computation of the modes is not a concern, as they solved the associated quadratic eigenvalue problem by the standard eigensolver eigs from Matlab that employs Arnoldi iterations.

In this paper our primary concern is to develop and study efficient eigensolvers for the spectral approximation of the damped vibration modes. Two approximations are investigated, one constructed from the edge-based displacement space (cf. Eq.(2.11)below), which results in quadratic eigenvalue problems (QEPs) and one from the node-based pressure space (cf. Eq.(2.12)), which results in rational eigenvalues problems (REPs). Our first approximation is identical to that in[3,5], but we further develop efficient methods for solving the associated QEP. However, we show in Section2that this problem has a large zero-frequency or null space and this fact may influence the efficiency of Arnoldi-type algorithms. Motivated by this, we extensively explore the second approximation of using the pressure space, which has a much smaller eigenvalue system to solve and which has a one dimensional null space. Instead of a QEP, the associated eigenvalue problem is a rational one having the form R(k)p = 0, where RðkÞ :¼ k2_{M þ} k2

aþbkA þ K is the rational k-matrix with coefficient matrices M, A and K.

While there is an extensive literature on QEPs problems[26], REPs are much less studied[25,27,28]. Although on the surface the rational R(k)p = 0 could be turned into a cubic one by multiplying out the denominator, we will preserve its rational structure and design efficient methods to numerically solve it in Section3.

The organization of this paper is as follows. We describe the underlying model fluid–solid problem of this paper in Section

2, where the edge-based displacement approximation and the node-based pressure approximation are derived. We pay par-ticular attention to identifying the dimension of the associated null space, which may influence performance of the numer-ical method introduced later. In Section3, we use the general strategy of turning a nonlinear eigenvalue problem into a standard one by some sort of linearization techniques. We then apply the Arnoldi type algorithms to solve it. For the two nonlinear eigenvalue problems, the QEP is as usual turned into a generalized eigenvalue problem (GEP), from which two types of standard eigenvalue problems (SEP)(3.7.1) and (3.7.2)are derived. The REP is trimmed-linearized into two types of three by three block SEPs(3.19.1) and (3.19.2). The important issue of residual error bound analysis is addressed here. We then apply Arnoldi method with Schur-restarting described in Section4to the resulting SEPs. The important issues of stopping criteria and computational costs for applying Arnoldi method to the QEP and REP are also derived in this section. In Section5, we present numerical results and evaluate the merits of the schemes involved where we also demonstrate the role of normwise scaling in preprocessing the eigenvalue problems. Conclusions are included in Section6.

2. Model problem

Let us consider a simple model of a rigid container filled with an inviscid compressible barotropic fluid and its acoustic energy is absorbed through a thin layer of a viscoelastic material applied to some or all of its walls. For simplicity we assume the fluid domainX Rn(n = 2 or 3) to be polyhedral, and the boundary oX=CA[CR, where the absorbing boundaryCAis

the union of all the different faces ofXand is covered by damping material. The rigid boundaryCRis the remaining part ofC.

An example of the setup is inFig. 1(i) on Section5, where the top boundary is absorbing and the remaining boundary is rigid. The dynamic variables of our model problem are the fluid pressure P and the displacement field U, which satisfy ([6,15])

q

@

2_U

@t2þ

r

P ¼ 0 in

X

; ð2:1Þ

P ¼ 

q

c2_{div U in}

_X

_;

ð2:2Þ

P ¼

a

U  n þ b@U @t  n

 

on

C

A; ð2:3Þ

U  n ¼ 0 on

C

R: ð2:4Þ

Here

q

is the fluid density, c, the acoustic speed, and n, the unit outer normal vector along oX. At the absorbing boundary

(2.3)indicates that the pressure is balanced by the effects of the viscous damping (the b term) and the elastic behavior (the

a

term). We assume the coefficients

a

and b are given positive constants.

To look for the damped vibration modes we assume(2.1)–(2.4)has complex solution of the form U(x, t) = ektu(x) and P(x, t) = ekt_{p(x). This leads to a problem of finding k 2 C; u :}

X! Cnand p :X! C; ðu; pÞ – ð0; 0Þ such that

q

k2u þ

r

p ¼ 0 in

X

; ð2:5Þ

p ¼ 

q

c2_{div u in}

_X

_{; ð2:6Þ}

p ¼ ð

a

þ kbÞu  n on

C

A; ð2:7Þ

u  n ¼ 0 on

C

R: ð2:8Þ

The boundary condition(2.7)makes this eigenvalue problem nonlinear. For each damped vibration mode,

x

:¼ Imk is the vibration angular frequency and Rek the decay rate. In practice, we select a range of

x

values and are interested in the least decaying modes in this range. We next describe the natural variational formulation of the above problem on which the numerical approximation will be based.

Let

V :¼ f

v

2 Hðdiv;

X

Þ :

v

 n 2 L2ð@

X

Þ and

v

 n ¼ 0 on

C

Rg:

Here we employ standard Sobolev spaces notation. For example, H( div,X) stands for the space of all L2_{vector functions v on}

Xwith L2integrable divergence.

Testing(2.5)by 

v

2 V and integrating by parts, we obtain a variational formulation of problem(2.5)–(2.8)involving only the displacement variable: Find k 2 C and u 2 V; u – 0, such that

Z X

q

c2_{div u div }

_v

_þZ CA

a

u  n

v

 n þ k Z CA bu  n

v

 n þ k2 Z X

q

u  

v

¼ 0

8

v

2 V: ð2:9Þ

(We use the symbol " to mean ‘for all’ throughout the paper.) This is a quadratic eigenvalue problem. Note that k = 0 is an eigenvalue and the dimension of its eigenspace

K :¼ fu 2 V : div u ¼ 0 in

X

and u  n ¼ 0 on @

X

g

is infinity. All nonzero eigenvalues have finite multiplicity (the dimension of the eigenspace is finite)[4]. It is shown in[4]

that all the other solutions of(2.9), the decay rate is strictly negative. That is, if an eigenpair 0 – k 2 C and 0 – u 2 V is a solution of problem(2.9)then Rek < 0.

Alternatively we can derive a variational formulation in terms of the pressure: Find k 2 C and p 2 H1(X) such that

Z X

r

p 

r

_{q þ} k 2

a

þ kb Z CA

q

pq þk 2 c2 Z X pq ¼ 0

8

q 2 H1ð

X

Þ: ð2:10Þ

However, in this case the eigenvalue problem is rational, which is rarely studied compared with linear and quadratic eigen-value problems. Note that in contrast to the displacement formulation, the eigenspace corresponding to k = 0 is now one dimensional. Thus this formulation has a much smaller null space or kernel, which may be more stable and efficient when used in conjunction with projection-like spectral approximation methods.

2.1. Spectral approximation

We now turn to the finite element methods for approximating the solutions of the quadratic eigenvalue problem(2.9)

and the rational eigenvalue problem(2.10). Spurious modes are usually present when standard finite elements are used in a displacement formulation. However Bermúdez et al.[4]successfully demonstrated that the spurious modes can be avoided by using the lowest order Raviart–Thomas elements in Rn_{;n ¼ 2; 3 (see, for instance,[7,18]). For simplicity we will}

consider only the two dimensional case. Let fThg be a regular family of triangulations ofXindexed by h, the maximum

diameter of the elements. Let

Vh:¼ f

v

h2 Hðdiv;

X

Þ :

v

hjT2 P n

0 P0x

8

T 2 Th and

v

h n ¼ 0 on

C

Rg  V

where n = 2 and Pkdenotes the set of polynomials of degree at most k. Thus locally vhtakes the form (a + sx, b + sy)>(>

denotes transpose and x = (x, y)>_{). The discrete problem associated with(2.9)is: Find k 2 C and u}

h2 Vh;uh–0, such that Z X

q

c2_{div u} hdiv 

v

hþ Z CA

a

uh n

v

h n þ k Z CA buh n

v

h n þ k2 Z X

q

uh 

v

h¼ 0;

8

v

h2 Vh: ð2:11Þ

Theorem 1. The dimension of the zero eigenspace E0 associated with (2.11) equals the number of interior nodes in the

triangulation.

Proof. Setting vh= uhand k = 0 in(2.11), we see that

div uh¼ 0 on

X

and uh n ¼ 0 on @

X

:

Since uh= (a + sx, b + sy)>on T 2 Th, the divergence free condition implies that uhis a constant vector (a, b)>on T. By direct

computation, we see that there exists a linear polynomial

w

Tsuch that

@wT

@x ¼ b and @wT

@y ¼ a:

Let n = (n1, n2)>be a unit normal to an edge e of T, so t = (n2, n1)>is a unit tangent vector to e. We see that

uh n ¼

r

wT t ¼ @wT

@t :

So if an edge e is common to T1and T2then in general wT1and wT2differ by a constant only by the continuity of uh n across e.

At an interior node Nj, we can assign a common value for all

w

Tat that node. Here T are all triangles sharing Njas the common

node. We then spread this defining process outward to allXusing the induced values on other nodes. Consequently,Wis continuous piecewise linear overX. Letr?_{:¼ ð}@

@y;@x@Þ >

and define

r

?_S

h:¼ f

r

?

W

h:

W

h is continuous piecewise linear and vanishes on the boundaryg:

Thus we have just shown the zero eigenspace E0is containedr\Shand the opposite inclusion is also easily checked. Hence

E0¼

r

?Sh:

We now find the dimension ofr\

Sh. Let N be the number of interior nodes and letWj, j = 1, . . . , N, be the nodal basis functions

such that Wj(Nk) = dkj. The linear independence of Wj’s is preserved by the perp-gradient operation. In fact, suppose

PN

j¼1cjr?Wj¼ 0. Then this impliesPNj¼1cjWj¼ c for some constant c. Hence cj= c by the conditionWj(Nk) = dkj. Consequently,

cðP_jWj 1Þ ¼ 0. But we knowPNj¼1Wj–1 due to the vanishing boundary condition. Thus cj= c = 0 and we conclude that the

dimension of the zero eigenspace dimE0¼ dimr?Shequals the number of interior nodes in the mesh. h

Define the conforming P1finite element space

Hh:¼ fph2 H 1

ð

X

Þ : phjT2 P1

8

T 2 Thg:

This is the subspace of H1₍

X) consisted of continuous piecewise linears. The alternative discrete problem in terms of the approximate pressure field is: Find k 2 C and ph2 Hhsuch that

Z X

r

ph

r

qhþ k2

a

þ kb Z CA

q

phqhþ k2 c2 Z X phqh¼ 0

8

qh2 Hh: ð2:12Þ

Letting qh= phand k = 0 in(2.12)we can easily see that the dimension of the zero eigenspace in this case is one, which is the

same as the original problem(2.10).

Again we see that the pressure formulation has a much smaller null space than the displacement formulation. Also the number of unknowns is much smaller. Thus the pressure formulation turns out to be a very good alternative, once in addi-tion we show in the remaining secaddi-tions that its associated eigenvalue problem can be efficiently solved. A minor remark is in order here.

Remark 2.1. Suppose an eigenpair (k, ph),k – 0 has been computed, what if, in addition, one wants to know a corresponding

displacement approximation uh? One must not find uh by solving an additional system linear equations again so as to

maintain the advantage of the pressure formulation. It should be given by a simple formula. A naive way is to use the relation

(2.5)to evaluate a uh, but this would be ill conceived since the computed displacement would be piecewise constant.

Consequently,r uh= 0, which certainly does not approximate(2.6). Fortunately, a general principle for such a problem

(recovery of uhfrom the pressure approximation ph) has been provided in[9]where one can obtain an accurate uhin the

Raviart–Thomas space by a simple evaluation formula which is a modification of the above naive formula.

3. Linearization of nonlinear eigenvalue problems

In this section we start to address the computational issues related to the displacement approximation(2.11)and the pressure approximation(2.12).

3.1. Linearization of quadratic eigenvalue problem

Suppose the total number of interior and absorbing edges is n1. Let f/jg n1

j¼1denote the cardinal basis of Vh, so that on the

edge ej, /jhas the unit normal flux and zero normal flux on the remaining n1 1 edges. That is,

R

ei/j nid

r

¼ dij. For uh2 Vh,

we write uh¼Pnj¼11 uj/jand denote u ¼ ½u>1;   ; u>n1

>_{. Note that the unknown vector u contains normal fluxes in its}

compo-nents. Then, the discrete problem(2.11)can be expressed as the following QEP:

Q ðkÞu  k2Muu þ ð

a

þ kbÞAuu þ Kuu ¼ 0; ð3:1Þ

where Mu ½Muij and Ku ½Kuij are mass and stiffness matrices, respectively, and Au ½Auij is used to describe the effect of the

absorbing wall. Here

Mu ij¼ Z X

q

/i /j; K u ij¼ Z X

q

c2_{div /} idiv /j; A u ij¼ Z CA /i n /j n; ð3:2Þ

for i, j = 1, . . . , n1. For this problem, we are only interested in eigenvalues that are located in the interior of the spectrum.

Sup-pose that the eigenvalues near

r

are of interest. Accordingly, the QEP(3.1)is shifted into

l

2_Me uþ

l

Deuþ eKu   u ¼ 0 ð3:3Þ with

l

= k 

r

and e Mu¼ Mu; eDu¼ 2

r

Muþ bAu; eKu¼

r

2Muþ ð

a

þ

r

bÞAuþ Ku: ð3:4Þ

On the one hand, one can numerically solve(3.3)without transforming it further. Among such direct methods we mention the SOAR (second-order Arnoldi) algorithm[1]and the Jacobi–Davidson algorithm applied to polynomial eigenvalue prob-lems[21]. On the other hand, it is more common to transform or linearize(3.3)into a SEP[26]. In this paper, we let

Au¼ 0  eMu I  eDu " # ; Bu¼ I 0 0 Keu   ð3:5Þ

and linearize(3.3)into the GEP

Aux ¼ 1

l

Bux with x  

l

Meuu u " # 

v

u   : ð3:6Þ

The matrix eKuin(3.5)is nonsingular because the shift value

r

is not an eigenvalue of(3.1). Furthermore, the GEP(3.6)can

then be transformed into two types of SEPs of the forms B1

u Aux ¼

l

1x and AuB1u y ¼

l

1y, respectively, where y ¼ Bux.

Therefore, from(3.5) and (3.6)we have

ðQ  SEP1Þ B1 u Au

v

u   ¼ 0  eMu e K1 u  eK1u Deu " #

v

u   ¼1

l

v

u   ð3:7:1Þ and ðQ  SEP2Þ AuB1u

v

w   ¼ 0  eMuKe 1 u I  eDuKe1u " #

v

w   ¼1

l

v

w   ; w ¼ eKuu: ð3:7:2Þ

Note that the SEPs of(3.7.1) and (3.7.2)derived by the QEP in(3.3), are called Q-SEP1 and Q-SEP2, respectively. The standard Arnoldi method can then be applied to solve Q-SEPs, and the details will be given in Section4.

3.2. Trimmed linearization method for rational eigenvalue problem Let fwjg

n2

j¼1be a nodal basis of Hh. For ph 2 Hh, we write ph¼

Pn2

j¼1pjwjand denote p ¼ ½p1;   ; pn2

>

. Then, the discrete problem(2.12)can be written as the following REP:

RðkÞp  k 2 c2Mpþ Kpþ k2 kbþ

a

Ap ! p ¼ 0; ð3:8Þ

where Mp ½Mpij and Kp ½Kpij are mass and stiffness matrices, respectively, and Ap ½Apij describes the effect of the

absorb-ing wall. Here,

Mp ij¼ Z X wiwj; K p ij¼ Z X

r

wi

r

wj; A p ij¼ Z CA

q

wiwj ð3:9Þ for i, j = 1, . . . , n2.

To solve REP(3.8), one approach is to multiply equation(3.8)by the scalar kb +

a

and expand it into a cubic polynomial eigenvalue problem, and then solve it by Jacobi–Davidson method[14]. An alternative approach is to treat(3.8)as nonlinear eigenvalue problem and solve it by a nonlinear eigensolver, such as Newton’s method, nonlinear Arnoldi method, or nonlin-ear Jacobi–Davidson method[20,27,28]. Recently, a trimmed linearization is proposed in[25]which linearizes(3.8)into a GEP so that the standard Arnoldi method can be applied. We introduce the trimmed linearization below.

Given a shift value

r

. With

l

= k 

r

, the rational k-matrix R(k) in(3.8)can be rewritten as

RðkÞ ¼ðk 

r

þ

r

Þ 2 c2 Mpþ Kpþ ðk 

r

þ

r

Þ2 ðk 

r

þ

r

Þb þ

a

Ap ¼ðk 

r

Þ 2 þ 2ðk 

r

Þ

r

þ

r

2 c2 Mpþ Kpþ ðk 

r

Þ2þ 2ðk 

r

Þ

r

þ

r

2 ðk 

r

Þb þ

r

bþ

a

Ap ¼

l

2 1 c2Mp   þ

l

2

r

c2Mp   þ

r

2 c2Mpþ Kp   þ

l

2_{þ 2}

_lr

_þ

_r

2

l

bþ

r

bþ

a

Ap ¼

l

2 1 c2Mp   þ1

l

2

r

c2Mp   þ 1

l

2

r

2 c2Mpþ Kp   þ1 þ 2

r

=

l

þ

r

2₌

_l

2

l

bþ ð

r

bþ

a

Þ Ap   : ð3:10Þ

Setting

m

= 1/

l

, the rational term in(3.10)can be simplified into the following form by applying the long division

r

2

_m

2_{þ 2}

_rm

_{þ 1} b=

m

þ ð

r

bþ

a

Þ¼

r

2

r

bþ

a

m

2_þ

r

2bþ 2

ra

ð

r

bþ

a

Þ2

m

þ

a

2 ð

r

bþ

a

Þ3 ð

r

bþ

a

Þ3

a

2 þ ð

r

bþ

a

Þ4

a

2_b

m

!1 :

This implies that

RðkÞ ¼ 1

m

2

m

2

r

2 c2Mpþ Kpþ

r

2

r

bþ

a

Ap   þ

m

2

r

c2Mpþ

r

2_{bþ 2}

_ra

ð

r

bþ

a

Þ2Ap ! þ 1 c2Mpþ

a

2 ð

r

bþ

a

Þ3Ap ! "  ð

r

bþ

a

Þ 3

a

2 þ ð

r

bþ

a

Þ4

a

2_b

m

!1 Ap 3 5 ¼

l

2_Me pþ

l

Depþ eKp

l

2#

.l

1 1 LpR>p; ð3:11Þ where e Mp¼ 1 c2Mpþ

a

2 ð

r

bþ

a

Þ3Ap; ð3:12Þ e Dp¼ 2

r

c2Mpþ

r

2 bþ 2

ra

ð

r

bþ

a

Þ2 Ap; ð3:13Þ e Kp¼

r

2 c2Mpþ Kpþ

r

2

r

bþ

a

Ap; ð3:14Þ #¼ð

r

bþ

a

Þ 3

a

2 ;

.

¼  ð

r

bþ

a

Þ4

a

2_b ð3:15Þ

and LpR>p ¼ Ap is the full-rank decomposition of Apwith Lp;Rp2 Rn2m. Introducing an auxiliary vector

q ¼ # 

.l

1 1_R>

pp; ð3:16Þ

the REP in(3.8)can be reformulated as

l

2_Me

pþ

l

Depþ eKp

 

p 

l

2_L

pq ¼ 0: ð3:17Þ

Using(3.16) and (3.17), we get the GEP

Apx  0  eMp Lp In2  eDp 0 0 R>p #Im 2 6 6 4 3 7 7 5x ¼

l

1 In2 0 0 0 Kep 0 0 0

.

Im 2 6 4 3 7 5x 

_l

1Bpx; ð3:18Þ where x ¼ ½ðð

l

1_Ke

pþ eDpÞpÞ>;p>;q>>. As before, the matrix eKpin(3.14)is nonsingular because the shift value

r

is not an

eigenvalue of(3.8). As in(3.7.1) and (3.7.2), the GEP(3.18)can then be, respectively, transformed into the following two types of the SEPs of the forms B1

ðR  SEP1Þ B1 p Apx ¼ 0  eMp Lp e K1 p  eK1p eDp 0 0 

.

1_R> p

.

1#Im 2 6 6 4 3 7 7 5x ¼

l

1x; ð3:19:1Þ and ðR  SEP2Þ ApB1p y ¼ 0  eMpKe1p

.

1Lp In2  eDpKe 1 p 0 0 R>pKe1p

.

1#Im 2 6 6 4 3 7 7 5y ¼ 1

l

y; y ¼ Bpx: ð3:19:2Þ

Note that the SEPs of(3.19.1) and (3.19.2)derived by the REP in(3.17)are called R-SEP1 and R-SEP2, respectively.

3.3. Error analysis

In this subsection, we will discuss residuals of QEP(3.1) and REP (3.8) by using linearizations(3.7.1) and (3.19.1), respectively.

We first derive residual bounds of approximate eigenpairs for QEP(3.1)by by using linearizations Q-SEP1 and Q-SEP2, respectively. Let ð

l

1

1 ;½v>1;u>1 >

Þ be an approximate eigenpair of(3.7.1)and ½f>₁₁;f>₁₂>be the associated residual vector. That is, f11 f12   ¼ _e0  eMu K1 u  eK1u Deu " #

v

1 u1   

_l

1 1

v

1 u1   ¼

_l

1 1 

v

1

l

1Meuu1 e K1 u ð

l

1

v

1

l

1Deuu1 eKuu1Þ " # : It follows that

l

2 1Meuu1þ

l

1Deuu1þ eKuu1¼

l

1ð

v

1

l

1f11Þ þ

l

1

v

1

l

1Keuf12¼ 

l

21f11

l

1Keuf12:

Therefore, with k1=

l

1+

r

and from(3.1)we have

kQ ðk1Þu1k ku1k ¼k

l

2 1Meuu1þ

l

1Deuu1þ eKuu1k ku1k 6j

l

1j 2 kf11k þ j

l

1jk eKukkf12k ku1k : ð3:20Þ

On the other hand, let ð

l

1 2 ;½v>2;w>2

>

Þ be an approximate eigenpair of(3.7.2)and ½f> 21;f

> 22

>_{be the associated residual vector.}

That is, f21 f22   ¼ 0  eMuKe 1 u I  eDuKe1u " #

v

2 w2   1

l

2

v

2 w2   ¼  eMu e K1 u w2l12

v

2

v

2 eDuKe1u w2_l1 2w2 2 4 3 5: It follows that

l

2 2MeuKe1u w2þ

l

2DeuKe1u w2þ w2¼

l

2ð

v

2

l

2f21Þ þ

l

2

v

2

l

2f22¼ 

l

22f21

l

2f22:

Letting u2¼ eK1u w2, with k2 =

l

2+

r

and from(3.1)we have,

kQ ðk2Þu2k ku2k ¼k

l

2 2Meuu2þ

l

2Deuu2þ eKuu2k ku2k 6j

l

2j 2 kf21k þ j

l

2jkf22k ku2k : ð3:21Þ

Now, we derive residual bounds of approximate eigenpairs for REP(3.8)by using linearizations R-SEP1 and R-SEP2, respec-tively. Let ð

l

1

1 ;½s>1;p>1;q>1 >

Þ be an approximate eigenpair of(3.19.1)and ½g>

11;g>12;g>13

>_{be the associated residual vector. That}

is, g11 g12 g13 2 6 4 3 7 5 ¼ 0  eMp Lp e K1 p  eK1p Dep 0 0 

.

1_R> p

.

1#Im 2 6 6 4 3 7 7 5 s1 p1 q1 2 6 4 3 7 5 

_l

1 1 s1 p1 q1 2 6 4 3 7 5:

This implies that

s1¼ 

l

1Mepp1þ

l

1Lpq1

l

1g11; ð3:22Þ g12¼ eK1p s1 eK1p Depp1 1

l

₁p1; ð3:23Þ q1¼

l

1

.

1_{# 1}  1

l

₁ g13þ

.

1R > pp1   : ð3:24Þ

s1¼ 

l

1Mepp1þ

l

21

l

1

.

1# 1

 1

Lpg13þ

.

1LpR>pp1

 



l

1g11: ð3:25Þ

Substituting(3.25)into(3.23), from(3.11) and (3.15)we get

Rð

r

þ

l

1Þp1¼

l

21Mepp1þ

l

1Depp1þ eKpp1

l

21 ð

r

bþ

a

Þ3

a

2 þ ð

r

bþ

a

Þ4

a

2_b

_l

1 " #1 LpR>pp1 ¼ 

l

2 1g11

l

1Kepg12

l

2 1 b

r

bþ

a

þ 1

l

₁  1 Lpg13

which implies that

kRð

r

þ

l

1Þp1k kp1k 6 1 kp1k j

l

1j 2 kg11k þ j

l

1jk eKpk kg12k þ

l

21 b

r

bþ

a

þ 1

l

1  1 kLpk kg13k ( ) : ð3:26Þ

On the other hand, let ð

l

1

2 ;½s>2;t>2;q>2 >

Þ be an approximate eigenpair of(3.19.2)and ½g>

21;g>22;g>23 >

be the associated residual vector. That is,

g21 g22 g23 2 6 4 3 7 5 ¼ 0  eMpKe1p

.

1Lp In2  eDpeK 1 p 0 0 R> pKe1p

.

1#Im 2 6 6 4 3 7 7 5 s2 t2 q2 2 6 4 3 7 5 

_l

1 2 s2 t2 q2 2 6 4 3 7 5:

This implies that

g21¼  eMpeK1p t2þ

.

1Lpq2 1

l

2 s2; ð3:27Þ s2¼ eDpKe1p t2þ 1

l

2 t2þ g22; ð3:28Þ q2¼

.

1# 1

l

₂  1 R> peK1p t2þ g23   : ð3:29Þ

Substituting(3.28) and (3.29)into(3.27), we have

l

2 2MepKe1p t2þ

l

2eDpKe1p t2þ t2

l

22 #

.l

12  1 LpR>pKe1p t2¼ 

l

22g21

l

2g22þ

l

2 2 #

.l

12  1 Lpg23:

Letting p2¼ eK1p t2, from(3.11)we get

Rð

r

þ

l

2Þp2¼

l

22Mepp2þ

l

2Depp2þ eKpp2

l

22 ð

r

bþ

a

Þ3

a

2 þ ð

r

bþ

a

Þ4

a

2_b

_l

2 " #1 LpR>pp2 ¼ 

l

2 2g21

l

2g22þ

l

22

a

2_b ð

r

bþ

a

Þ4 b

r

bþ

a

þ 1

l

2  1 Lpg23: Hence, kRð

r

þ

l

2Þp2k kp2k 6 1 kp2k j

l

₂j2kg21k þ j

l

2jkg22k þ

l

22

a

2_b ð

r

bþ

a

Þ4 b

r

bþ

a

þ 1

l

2  1 kLpkkg23k ( ) : ð3:30Þ

Remark 3.1. In order to check the tightness of upper bounds in(3.20) and (3.21), as well as,(3.26) and (3.30)for residuals, respectively, we refer to the coefficient matrices generated inExample 1of Section5. For(2.9)we adopt the data as in[4]by setting

q

= 1 kg/m3_{, c = 340 m/s,}

_a

_{= 5  10}4_N/m3_{, and b = 200 Ns/m}3_{. In addition, we choose}

_r

_{¼ 25 þ 600}

_p

_{ıðı }pffiffiffiffiffiffiffi_{1Þ as}

the shift value. Then

(i) from(3.2), the element mass and stiffness matrices are

h2 6

q

2 1 0 1 2 0 0 0 2 2 6 4 3 7 5;

q

c2 2 2 2pffiffiffi2 2 2 2pffiffiffi2 2pffiffiffi2 2pffiffiffi2 4 2 6 4 3 7 5;

respectively. Hence, by(3.4)the infinity norm of eKucan be estimated by k eKuk1 kKuk1¼ Oð

q

c2Þ ¼ Oð10 5

Þ. From(3.20)and

(3.21), we conclude that the upper bound for the residual of the approximate eigenpair (

l

1+

r

, u1) of(3.1)by solving Q-SEP1

(ii) From(3.9), the element mass and stiffness matrices are h2 24 2 1 1 1 2 1 1 1 2 2 6 4 3 7 5; 1 1=2 1=2 1=2 1=2 0 1=2 0 1=2 2 6 4 3 7 5;

respectively. Hence, by(3.14)we have that k eKpk1 kKpk1¼ Oð1Þ. If the eigenvalue k is one of the desired eigenvalues in

Fig. 2, then with

l

= k 

r

we have

4  107_<

_l

2 b

r

bþ

a

þ 1

l

 1 <3:1  10 10_; 0:001 <

l

2

a

2b ð

r

bþ

a

Þ4 b

r

bþ

a

þ 1

l

 1 <0:8:

Clearly, from(3.26) and (3.30)we conclude that the upper bound for the residual of the approximate eigenpair (

r

+

l

1, p1) of

REP(3.8)by solving R-SEP1 is larger than that of (

r

+

l

2, p2) of(3.8)by solving R-SEP2.

4. Arnoldi method with schur-restarting

The Arnoldi method is the most popular method for solving large sparse SEPs:Lx ¼ kx. In Arnoldi process, an orthonormal matrix Vk+1is generated to satisfy

LVk¼ VkHkþ hkþ1;k

v

kþ1e>k; ð4:1Þ

where Hk2 Rkkis upper Hessenberg. If the dimension of the Krylov subspace span {Vk} is larger than a certain value, then

the process of Arnoldi decomposition will be restarted.

For the restarting process, we can use an implicit restart scheme[17,22]. The package ARPACK[16]includes a very suc-cessful implementation of the implicitly restarted Arnoldi algorithm. It has been used by many engineering fields and re-mains a popular choice for solving eigenvalue problems. However, these implicitly restart type schemes may suffer from numerical instability due to rounding errors. Stewart proposed the Krylov–Schur method[13,23,24]that relaxes the need to preserve the structure of the Arnoldi decomposition and therefore ease the complications of the purging and deflating.

We state the Schur-restarting scheme as follows. Let

Hk¼ U½ m U‘ Tm Tf 0 T‘   _U m U ‘   ð4:2Þ

be a Schur decomposition of Hkwhere Tmand T‘are upper triangular, and the eigenvalues of Tmare of interest. Here and

hereafter U⁄_{denotes the conjugate transpose of the matrix U. Substituting(4.2)into(4.1), we see that}

−3500 −300 −250 −200 −150 −100 −50 0 500 1000 1500 2000 2500 3000 3500 4000 real part imaginary part

λ

₁

λ

₂

λ

₃

λ

₄

λ

5

λ

₆

λ

₇

λ

8

λ

₉

λ

₁₀

L Vð k½Um U‘Þ ¼ Vð k½Um U‘Þ Tm Tf 0 T‘   þ hkþ1;k

v

kþ1 e>k½Um U‘  ;

which implies that

L eVm¼ eVmTmþ ~

v

mþ1tm; ð4:3Þ

where eVm VkUm; ~

v

mþ1¼

v

kþ1and tm hkþ1;ke>kUm. Let Q1be a Householder matrix with

t

mQ1¼

s

e>m:

Then(4.3)can be rewritten as

Lð eVmQ1Þ ¼ ð eVmQ1ÞðQ

1TmQ1Þ þ

s

v

~mþ1e>m: ð4:4Þ

The matrix Q

1TmQ1can be reduced to a new Hessenberg matrix eHmby using Householder matrices Qifor i = 2, . . . , m  1 with

Q m1   Q 2ðQ 1TmQ1ÞQ2   Qm1¼ eHm; e>jQ2   Qm1¼ e>m:

Multiplying(4.4)by Qi, i = 2, . . . , m  1, a new Arnoldi decomposition of order m

L eVm¼ eVmHemþ

s

v

~mþ1e>m

is obtained where eVm:¼ eVmQ1   Qm1and the Arnoldi process can be applied to generate it to order k in(4.1). One repeats

the above process until the desired eigenvalues are convergent. The process is summarized in Algorithm 1. Algorithm 1. Arnoldi method with Schur-restarting for solving Lx ¼ kx

Input: L: coefficient matrix, tolL: tolerance for convergence, rmax: maximum number of Schur-restartings.

Output: The desired m eigenpairs.

1: Build an initial Arnoldi decomposition of order m as in(4.1)and set r = 0. 2: restart

3: Extend Arnoldi decomposition of order m to order k = m + ‘ and set r = r + 1. 4: Compute all Ritz pairs ð

l

1

i ;ziÞ with Hkzi¼

l

1i zi,i = 1, . . . , k and sorting Ritz values so that{(

l

1, z1), . . . , (

l

m, zm)}

are wanted. 5: for i = 1, . . . , m

6: Check convergence by jhkþ1;kjje>kzij < tolL.

7: end for

8: if (Not all m desired eigenvalues are convergent andr < rmax) then

9: Compute the Schur decomposition of Hkas in(4.2), where the eigenvalues of Tmare of interest.

10: Set Vm:¼ VkUm, vm+1:¼ vk+1and tm:¼ hkþ1;ke>kUm.

11: Compute Householder transformation Q1such that tmQ1¼

s

e>m.

12: Reduce Q1TmQ1to a new Hessenberg matrix Hmby using Householder transformations Qifor i = 2, . . . , m  1.

13: Set Vm:¼ VmQ1  Qm1and hm+1,m=

s

to get the new Arnoldi decomposition with order m:

LVm¼ VmHmþ hmþ1;m

v

mþ1e>m: ð4:5Þ

14: end if

15: until (desired m eigenpairs are convergent or r P rmax)

Now, we will apply the Algorithm 1 to solve QEP(3.1)and REP(3.8), respectively, by setting L to be the coefficient matri-ces in(3.7.1) and (3.19.1), respectively.

4.1. Stopping criteria for QEPs and REPs Let (

l

1_{, z) be a Ritz pair and satisfy H}

kz =

l

1z. From(4.1)and Q-SEP1 in(3.7.1)we have

0  eMu e K1 u  eK1u Deu " # Vk1 Vk2   z ¼1

l

Vk1 Vk2   z þ hkþ1;k

v

kþ1;1

v

kþ1;2   e> kz; ð4:6Þ where Vk¼ V_Vk1 k2   and

v

kþ1¼

v

_v

kþ1;1 kþ1;2  

are partitioned with compatible sizes. Using the first equation of(4.6), we can elim-inate Vk1z in the second equation and get

Q ðkÞu1 k k u1 k k ¼ kð

l

2_Me uþ

l

eDuþ eKuÞu1k u1 k k ¼ j

l

jjhkþ1;kj e>kz f1 ku1k  q1ð

l

Þ; ð4:7Þ

where u1= Vk2z, k =

l

+

r

and f1¼ k

l

v

kþ1;1þ eKu

v

kþ1;2k. Without ambiguity by using the same notations as above in

Algo-rithm 1, from(4.1)and Q-SEP2 in(3.7.2)we also have

0  eMueK1u I  eDuKe1u " # Vk1 Vk2   z ¼1

l

Vk1 Vk2   z þ hkþ1;k

v

kþ1;1

v

kþ1;2   e> kz and Q ðkÞu2 k k ku2k ¼ kð

l

2_Me uþ

l

Deuþ eKuÞu2k ku2k ¼ j

l

jjhkþ1;kj e>kz f2 ku2k  q2ð

l

Þ; ð4:8Þ

where u2¼ eK1u Vk2z; k ¼

l

þ

r

and f2¼ k

l

v

kþ1;1þ

v

kþ1;2k. Therefore, q1(

l

) in(4.7)and q2(

l

) in (4.8), respectively, can be

used as stopping criteria for residuals while Algorithm 1 is applied to solved QEPs(3.1).

Similarly, we can apply Algorithm 1 to solve REPs(3.8). As above, we let (

l

1_{, z) be a Ritz pair and satisfy H}

kz =

l

1z. From

(4.1), and R-SEP1, R-SEP2 in(3.19.1)we have

0  eMp Lp e K1 p  eK1p Dep 0 0 

.

1_R> p

.

1#Im 2 6 6 4 3 7 7 5 Vk1 Vk2 Vk3 2 6 4 3 7 5z ¼

_l

1 Vk1 Vk2 Vk3 2 6 4 3 7 5z þ hkþ1;k

v

kþ1;1

v

kþ1;2

v

kþ1;3 2 6 4 3 7 5e> kz ð4:9Þ and 0  eMpKe1p

.

1Lp In2  eDpKe 1 p 0 0 R> pKe1p

.

1#Im 2 6 6 4 3 7 7 5 Vk1 Vk2 Vk3 2 6 4 3 7 5z ¼

_l

1 Vk1 Vk2 Vk3 2 6 4 3 7 5z þ hkþ1;k

v

kþ1;1

v

kþ1;2

v

kþ1;3 2 6 4 3 7 5e> kz; ð4:10Þ where Vk¼ ½V>k1;V > k2;V > k3 > and

v

kþ1¼ ½

v

>kþ1;1;

v

>kþ1;2;

v

>kþ1;3 >

are partitioned with compatible sizes. Using the first and the third equations of(4.9)and(4.10), we can eliminate V1z and V3z in the second equation of(4.9) and (4.10), respectively,

and get RðkÞp1 k k kp1k ¼k½

l

2_Me pþ

l

eDpþ eKp

l

2ð# 

.l

1Þ1App1k kp1k ¼j

l

jjhkþ1;kj e > kz n1 kp1k  r1ð

l

Þ; ð4:11Þ where p1= Vk2z, k =

l

+

r

and n1¼ k

l

v

kþ1;1þ eKp

v

kþ1;2 .l 2 #l.Lp

v

kþ1;3k, and kRðkÞp2k kp2k ¼k

l

2_Me pþ

l

Depþ eKp

l

2ð# 

.l

1Þ1Ap h i p2k kp2k ¼j

l

jjhkþ1;kj e > kz n2 kp2k  r2ð

l

Þ; ð4:12Þ where p2¼ eK1p Vk2z; k ¼

l

þ

r

and n2¼ k

l

v

kþ1;1þ

v

kþ1;2 l 2

#l.Lp

v

kþ1;3k. Therefore, r1(

l

) in(4.11)and r2(

l

) in(4.12)can be

used as stopping criteria for residuals while Algorithm 1 is applied to solve REPs(3.8).

Applying Algorithm 1 to solve QEPs(3.1)and REPs(3.8)are summarized in Algorithms 2 and 3, respectively.

Algorithm 2. Arnoldi method with Schur-restarting for solving QEP in(3.1)

Input: Coefficient matrices Mu, Duand Ku, parameters c,

a

and b,

r

: shift value, tolQ: tolerance for convergence, rmax:

maximum number of Schur-restartings.

Output: The desired eigenpairs (ki, ui) fori = 1, . . . , m.

1: Construct matrices eMu; eDuand eKudefined in(3.4)and set r = 0.

2: Compute initial Arnoldi decomposition in Line 1 of Algorithm 1 with L in Q-SEP1 or Q-SEP2. 3: restart

4: Do the steps in Lines 3 and 4 of Algorithm 1. 5: for i = 1, . . . , m do

6: Compute

u

ð

l

iÞ ¼ ðj

r

þ

l

1i j 2

kMuk þ j

a

þ ð

r

þ

l

1i ÞbjkAuk þ kKukÞ.

7: Check convergence of QEP by q‘(

l

i)/

u

(

l

i) < tolQwith q‘(

l

i) in(4.7)or(4.8),‘ = 1,2.

8: end for

9: if (Not all m desired eigenvalues are convergent and r < rmax) then

10: Do the Schur-restarting in Lines 9–13 of Algorithm 1. 11: end if

12: until (desire m eigenpairs are convergent or r P rmax)

13: Set ki¼

r

þ

l

1i and ui= Vk2zifori = 1, . . . , m.

14: if Q-SEP2 is solved then 15: ui eK1u ui;i ¼ 1; . . . ; m.

Algorithm 3. Arnoldi method with Schur-restarting for solving REP in(3.8)

Input: Coefficient matrices Mp, Kpand Ap, parameters c,

a

and b,

r

: shift value, tolR: tolerance for convergence, rmax:

maximum number of Schur-restartings.

Output: The desired eigenpairs (ki, pi) for i = 1, . . . , m.

1: Construct matrices eMp; eDpand eKpdefined in(3.12), (3.13) and (3.14), respectively, and set r = 0.

2: Compute the full-rank decomposition of Ap:LpR>p ¼ Ap.

3: Compute initial Arnoldi decomposition in Line 1 of Algorithm 1 with L in R-SEP1 or R-SEP2. 4: restart

5: Do the steps in Lines 3 and 4 of Algorithm 1. 6: for i = 1, . . . , m do 7: Compute wð

l

iÞ ¼ j ðrþl1 i Þ 2 c2 jkMpk þ kKpk þ j ðrþl 1 i Þ 2 aþðrþl1 i Þb jkApk.

8: Check convergence by r‘(

l

i)/

w

(

l

i) < tolRwith r‘(

l

i) in(4.11)or(4.12),‘ = 1,2.

9: end for

10: if (Not all m desired eigenvalues are convergent and r < rmax) then

11: Do the Schur-restarting in Lines 9–13 of Algorithm 1. 12: end if

13: until (desire m eigenpairs are convergent or r P rmax)

14: Set ki¼

r

þ

l

1i and pi= Vk2zifor i = 1, . . . , m.

15: if R-SEP2 is solved then 16: pi eK1p pi;i ¼ 1; . . . ; m.

17: end if

4.2. Computational costs

In this subsection, we compare the computational costs of the jth Arnoldi step of Algorithm 1 for solving Q-SEPs(3.7.1)

and R-SEPs(3.19.1), respectively. This is of general interest, because a comparison of the CPU times is sensible only if the number of outer iterations of Algorithm 2 or Algorithm 3 is the same for each algorithm. From(4.5), the unit vector vj+1

is generated by

L

v

j¼ Xj

i¼1

hj;i

v

iþ hjþ1;j

v

jþ1

where hj;i¼

v

iLvj for i = 1, . . . , j and hjþ1;j¼ kLvjPji¼1hj;i

v

ik2. For convenience, we let

v

j¼ ½v>j1;

v

>j2 >_with

_v

ji2 Cn, i = 1,2.

The matrix–vector product L

v

jin Algorithm 2 for solving QEP(3.1)by Q-SEP1(3.7.1)and Q-SEP2(3.7.2)can be, respectively,

represented by B1 u Au

v

j¼  eMu

v

j2 e K1 u ð

v

j1 eDu

v

j2Þ " # ; AuB1u

v

j¼  eMugu

v

j1 eDugu " #

with gu¼ eK1u

v

j2. This implies that Algorithm 2 for Q-SEP1 and Q-SEP2 needs the same computational costs for generating

the unit vector vj+1for each j.

On the other hand, by letting

v

j¼ ½

v

>j1;

v

>j2;

v

>j3 >

with

v

ji2 Cn;i ¼ 1; 2 and

v

j32 Cm, the matrix–vector product L

v

jin

Algo-rithm 3 for solving REPs by R-SEP1(3.19.1)and R-SEP2(3.19.2)can be, respectively, represented by

B1p Ap

v

j¼ Lp

v

j3 eMp

v

j2 e K1 p ð

v

j1 eDp

v

j2Þ

.

1_#

_v

j3

.

1R>p

v

j2 2 6 6 4 3 7 7 5; ApB1p

v

j¼

.

1_L p

v

j3 eMpgp

v

j1 eDpgp

.

1_#

_v

j3 R>pgp 2 6 6 4 3 7 7 5 with g_p¼ eK1

p

v

j2. Consequently, the computational cost of ApB1p

v

jneeds an extra cost for the computation of

.

1(Lpvj3)

compared to that B1

p Ap

v

j. The cost for generating the unit vector vj+1by R-SEP1 is slightly cheaper than that by R-SEP2.

Remark 4.1. In the numerical implementation, the vectors gu¼ eK1u

v

j2 and gp¼ eK1p

v

j2 forj = 1, . . . , k can be saved in

Gu ½ eK1u

v

12   eK1u

v

k2 and Gp ½ eK1p

v

12   eK1p

v

k2, respectively, so that the vectors u2, p2 in (4.8) and (4.12) can be

computed byu2= Guz and p2= Gpz directly. Hence, it requires the same computational costs for computingu1, u2in(4.7) and

(4.8), as well as,p1, p2in(4.11) and (4.12), respectively. Consequently, the computational costs of Q-SEP1 for the convergence

test in Algorithm 2 need one extra matrix–vector product eKu

v

kþ1;2than those of Q-SEP2 in computing f1andf2. Similarly, the

computational costs of R-SEP1 for the convergence test in Algorithm 3 need one extra matrix–vector product eKp

v

kþ1;2than

those of R-SEP2 in computing n1 and n2. Therefore, we conclude that Algorithm 2 for Q-SEP1 and Q-SEP2, as well as,

Algorithm 3 for R-SEP1 and R-SEP2, respectively, almost have the same computational costs provided that they have the same outer iterations.

5. Numerical results

We conduct numerical experiments to evaluate performance and accuracy of the eigenvalue solvers described in Section

4. To distinguish between various eigenvalue problems, we use notations Q1, Q2, R1 and R2 defined as follows:

 Q1: Applying Algorithm 2 to solve the QEP(3.1)with Q-SEP1 in(3.7.1).  Q2: Applying Algorithm 2 to solve the QEP(3.1)with Q-SEP2 in(3.7.2).  R1: Applying Algorithm 3 to solve the REP(3.8)with R-SEP1 in(3.19.1).  R2: Applying Algorithm 3 to solve the REP(3.8)with R-SEP2 in(3.19.2).

All computations are carried out in MATLAB 2009a on a HP workstation with an Intel Quad-Core Xeon X5570 2.93 GHz and 72 GB main memory, using IEEE double-precision floating-point arithmetic. We apply Algorithms 2 and 3 to solve the following examples arising in fluid–solid systems. The order k of Arnoldi decomposition in Line 3 of Algorithm 1 is set k = 40, the maximum number rmaxof Schur-restartings is set rmax= 15 and the number of desired eigenpairs is m = 10. The relative

residuals of approximate eigenpairs (ki, ui) and (ki, pi) computed by Q1 and Q2, as well as, R1 and R2 are, respectively, defined

by

kQðkiÞuik

u

ðkiÞkuik

and kRðkiÞpik wðkiÞkpik

where

u

(ki) and

w

(ki) are given in Algorithms 2 and 3, respectively. Tolerances for relative residuals of QEPs and REPs are

chosen by tolQ= tolR= 5  1015. The linear systems in Algorithms 2 and 3 are solved by LU-factorization with the shift value

r

=25 + 600

p

ı. Fronbenius norm for matrices and2-norm for vectors are used.

Example 1 [4]. We take the geometrical data: the domainX= [0m,1m]  [  0.75m, 0m],CA= [0m, 1m]  {0m} given in

Fig. 1(i) and the following physical data:

q

= 1 kg/m3_{, c = 340 m/s,}

_a

_{= 5  10}4_N/m3_{, and b = 200 Ns/m}3_. Table 1

Computational costs of the jth Arnoldi step of Algorithm 1 for Q-SEP2 and R-SEP2, where e

Mu; eDu; eKu2 Rn1n1, eMp; eDp; eKp2 Rn2n2and Lp;Rp2 Rn2mwith m  n2. The length of the

vectors in the inner products for Q-SEP2 and R-SEP2 are 2n1and 2n2+ m, respectively.

Q-SEP2(3.7.2) R-SEP2(3.19.2)

Solving linear system eKuxu¼ bu Kepxp¼ bp

Matrix–vector products Meubu; eDubu Mepbp; eDpbp;Lpcp;R>pc>p

Inner products j + 1 j + 1

Saxpy operators j + 1 j + 2

Scale-vector product 1 1

Table 2

Dimension information and convergence rates of k1.

(M, N) Matrix size (QEP) Matrix size (REP) Conv. rate

(3M  1)  N (M + 1)  (N + 1) k1 Q2 R2 (48, 36) 5,148 1,813 (96, 72) 20,664 7,081 (192, 144) 82,800 27,985 rate[1,1] 1.9979 2.0010 (384, 288) 331,488 111,265 rate[1,2] 1.9995 2.0003 (768, 576) 1,326,528 443,713 rate[1,3] 1.9999 2.0001

The rectangular domainXis uniformly partitioned into M byN rectangles and each rectangle is further refined into two triangles, seeFig. 1(ii). The dimensions of coefficient matrices in QEP(3.1)and REP(3.8)are (3M  1)  N and (M + 1)  (N + 1), respectively.Fig. 2plots the analytic solutions of the desired eigenvalues k1, . . . , k10of(2.5)–(2.8)(see[4]) with the

lowest positive vibration frequencies satisfying 0 <ImðkiÞ

2p <600 Hz.

Convergence test: We first demonstrate convergence rates of Q2 and R2 while computing the desired eigenvalues in

Fig. 2. To measure the convergence rate, we run the test over the five successively refined meshes (See the first column ofTable 2) and then calculate the rates by

rate½i;j¼ log2

jk½i;j k½i;jþ1j

jk½i;jþ1 k½i;jþ2j

 

; for i ¼ 1; . . . ; 10; j ¼ 1; 2; 3;

where k[i,j]for j = 1, . . . , 5 denote the approximate eigenvalues computed by Q2 and R2 corresponding tokiobtained from the

meshes described inTable 2. The 5th and the 6th columns ofTable 2illustrate the quadratic convergence of rate[1,j]{j = 1,2,3

for k1of QEP(3.1)and REP(3.8), respectively. In our numerical experiment, the convergence rate are always close to 2 for all

desired eigenvalues, ki, i = 1,. . .,10, computed by Q2 and R2 as well as Q1 and R1.

Normwise scaling of QEP: Balancing norms of coefficient matrices is an important issue[26]before solving a QEP of the form:

Pð

l

Þx  ð

l

2_P

2þ

l

P1þ P0Þx ¼ 0: ð5:1Þ

In[10]authors give an elegant way to scale the norms of coefficient matrices of(5.1)as follows. Define

bPð

m

Þx  ð

m

2_bP

2þ

m

bP1þ bP0Þx ¼ 0

with

m

¼

l

=f; bP2¼ f2

g

P2; bP1¼ f

g

P1and bP0¼

g

P0, where f and

g

are scaling factors. Taking f and

g

as f¼

ffiffiffiffiffiffiffiffiffiffiffiffi

c

0=

c

2

p

and

g

¼ 2=ð

c

0þ

c

1fÞ with

c

2:¼ kP2k2;

c

1:¼ kP1k2;

c

0:¼ kP0k2, it is proved in[10]that the problem

min

f;g max jkbP2k2 1j; jkbP1k2 1j; jkbP0k2 1j

n o

achieves the optimum at f⁄ and

g

⁄. In our implementation, the values of

c

i, for i = 0, 1, 2 are computed by

c

₂¼ k eMukF;

c

1¼ k eDukF;

c

0¼ k eKukFand

c

2¼ k eMpkF;

c

1¼ k eDpkF;

c

0¼ k eKpkFfor QEP(3.3)and REP(3.17), respectively. We

de-note ‘‘#It’’ the number of Schur-restartings (outer iterations). InTable 3, we show # Its for computing 10 desired eigenvalues ofExample 1with (M, N) = (768, 576) by Q1, Q2, R1 and R2 with/without scaling. The tolerances tolQand tolRfor relative

residuals are chosen to be 5  1015_{. We see that the convergence rate of scaled Q-SEPs or R-SEPs is faster than that of}

un-scaled Q-SEPs or R-SEPs. The performance of Q2 and R2 is also better than that of Q1 and R1, respectively. In the case of unscaled REP, the norms of eMp; eDp and eKp in(3.12)–(3.14)are Oð1010Þ; Oð105Þ and Oð1Þ, respectively. Since the norms

of coefficient matrices vary too much, R1 can even fail to converge to 10 eigenpairs after 15 outer iterations.

No spurious eigenmodes: In[4], it has been proved that there are no spurious eigenmodes for the discretization based on Raviart–Thomas finite elements. We compute twenty desired eigenvalues of QEP(3.1)and REP(3.8)by Q2 and R2, respec-tively, with scaling and various mesh sizes as shown inTable 2(we computed 20 instead of 10 eigenvalues to be better con-firmed). The desired eigenvalues of REP are in one-to-one correspondence to those of QEP which match well with relative error less than 106_{, that is, no spurious eigenmodes ever appear. We numerically conclude that there are no spurious}

eigen-modes for the discretization in terms of pressure nodal finite elements.

Null space considerations:Theorem 1shows that the dimension of the null space of QEP(3.1)is equal to the number of interior nodes, i.e., (M  1)(N  1). In order to observe the interference of such a large null space in the convergence of Q1 and Q2, we give six different shift values denoted by the ‘‘+’’ inFig. 3to observe variation in the #Its for Q1 and Q2. The inte-ger pair (i, j) under each shift value ‘‘+’’ denotes the#Its for Q1 and Q2, respectively. The results inFig. 3shows that the #It needed decreases, as the shift value

r

is chosen relatively far away from zero.

Comparison of pressure and displacement formulations: In this paragraph, we shall discuss the advantages of using the nodal pressure finite elements with various mesh sizes described inTable 2. The notations ‘‘TQ2’’ and ‘‘TR2’’ denote the

total CPU times for Q2 and R2, respectively. We summarize the results as follows:

 Accuracy of eigenpairs: FromRemark 3.1, the upper bound for relative residual of the approximate eigepairs of QEP(3.1)

(or REP(3.8)) by using Q-SEP2(3.7.2)(or R-SEP2(3.19.2)) is much smaller than that by using Q-SEP1(3.7.1)(or R-SEP1

(3.19.1)). On applying Q1 and Q2 to solve QEP(3.1)with#It = 2, inFig. 4, we see that the relative residuals of eigenpairs

Table 3

#Its for k1, . . . , k10of Q-SEPs and R-SEPs with/without scaling.

Q1 Q2 R1 R2

#It (scaled) 3 2 4 3

Fig. 4. The relative residuals of computed eigenpairs, obtained by Q1, Q2 for QEP(3.1)and R1, R2 for REP(3.8)with (M, N) = (768, 576). −3500 −300 −250 −200 −150 −100 −50 0 500 1000 1500 2000 2500 3000 3500 4000 Real part Imaginary part ( 6, 5 ) ( 3, 2 ) ( 6, 4 ) ( 2, 2 ) ( 4, 4 ) ( 2, 2 ) Eigenvalues Shift values

Fig. 3. The #Its of Q1 and Q2 with different shift values. ‘‘o’’ denotes desired eigenvalues k1, . . . , k10. ‘‘(i, j)’’ denotes the #Its for Q1 and Q2, respectively.

Table 4

Iteration numbers and CPU times for Q2 and R2.

(M, N) Q2 R2 TR2 TQ 2 #It TQ2 #It TR2 (48, 36) 2 1.316 2 0.471 0.36 ( 96, 72) 2 7.717 2 2.387 0.31 (192, 144) 2 55.27 2 14.95 0.27 (384, 288) 2 567.8 2 134.0 0.24 (768, 576) 2 8152 2 1645 0.20

corresponding to k4and k5computed by Q2 are improved by about 1 significant digit than those by Q1. The other

eigen-pairs almost have the same accuracy. On applying R1 and R2 to solve REP(3.8)with#It = 2, inFig. 4, we see that the rel-ative residuals of eigenpairs computed by R2 are improved by about 2 to 4 significant digits than those by R1.  Comparison R2 with Q2: From Sub Section4.2we see that Q1 and Q2, as well as, R1 and R2 have the same computational

costs, respectively. FromFig. 4, we favor applying Q2 and R2 to solve QEP(3.1)and REP(3.8), respectively. From column 12 ofTable 4, we see that the CPU times for solving the REP(3.8)by R2 is only 1/5 to 1/3 of that for solving the QEP(3.1)

by Q2. The accuracy of the computed eigenpairs for REP(3.8)is also better than that of QEP(3.1). These results tell us that using R2 to solve nodal pressure finite elements for the discrete problem(2.12)is better than that using Q2 to solve Rav-iart–Thomas displacement finite elements for the discrete problem(2.11).

−6000 −500 −400 −300 −200 −100 0 500 1000 1500 2000 2500 3000 3500 4000

λ

₁

λ

₂

λ

₃

λ

₄

λ

₅

λ

₆

λ

₇

λ

₈

λ

₉

λ

₁₀ Real part Imaginary part

Fig. 5. The distribution of the ten desired eigenvalues k1, . . . , k10forExample 2.

We now want to apply our methods to a more complicated configuration in which the absorbing walls are located on three sides.

Example 2. We use the same geometric data and physical data inExample 1except that the absorbing wall is extended to one half of the rigid walls in the left and right boundaries, that isC= [0, 1]  {0} [ {0}  [  0.375, 0] [ {1}  [  0.375, 0].

InExample 1, we numerically show that there are no spurious eigenmodes for the discretization in terms of pressure nodal finite elements. Moreover, the computational cost for solving the associated REP(3.8)is obviously less than that of solving QEP (3.1) which is obtained from using Raviart–Thomas displacement finite elements to the discrete problem

(2.11). Therefore, in this example we only use nodal finite elements to discretize the model and compare the accuracy of R1 and R2 for solving the associated REP. The computed eigenvalues k1, . . . , k10with lowest positive vibration frequencies

satisfying 0 <ImðkiÞ

2p <600 Hz are shown inFig. 5. The convergence rates for k1, . . . , k10obtained from various the mesh sizes

described inTable 2are also close to 2. The relative residuals computed by R1 and R2 are presented inFig. 6which shows that the accuracy of the eigenpairs produced by R2 is better than R1.

6. Conclusions

We propose efficient methods for computing damped vibration modes of an acoustic fluid confined in a cavity with absorbing walls capable of dissipating acoustic energy. Two approximations are investigated, one constructed from the edge-based displacement space, which results in QEPs(3.1)and one from the node-based pressure space, which results in REPs(3.8). Our numerical results show that both nodal and edge-based finite elements have second-order convergence rate. We theoretically prove that the nullity of the QEP(3.1)equals the number of the interior grid points. Numerical results show that if the shift value is close to zero, then such a large null space interfere with the convergence of the eigensolver. Furthermore, numerical evidence also shows that there are no spurious eigenmodes for the discretization in terms of pres-sure nodal finite elements and the CPU times for solving the corresponding REP(3.8)are only 1/5 to 1/3 of the CPU times for solving the QEP(3.1).

For solving the nonlinear eigenvalue problems, a linearization and a trimmed-linearization methods are used to linearize QEP(3.1)and REP(3.8)into four different types of SEPs which can be solved by Q1 and Q2, as well as, R1 and R2. Numerical accuracy shows that Q2 and R2 algorithms are better than Q1 and R1, respectively.

Acknowledgments

The authors are grateful to the anonymous referees for their useful comments and suggestions. This work is partially sup-ported by the National Science Council and the National Center for Theoretical Sciences of Taiwan.

References

[1] Z. Bai, Y. Su, A second-order Arnoldi method for the solution of the quadratic eigenvalue problem, SIAM J. Matrix Anal. 26 (3) (2005) 640–659. [2] K.J. Bathe, C. Nitikitpaiboon, X. Wang, A mixed displacement-based finite element formulation for acoustic fluid–structure interaction, Comput. Struct.

56 (1995) 225–237.

[3] A. Bermúdez, R. Durán, M.A. Muschietti, R. Rodrı´guez, J. Solomin, Finite element vibration analysis of fluid–solid systems without spurious modes, SIAM J. Numer. Anal. 32 (1995) 1280–1295.

[4] A. Bermúdez, R.G. Durán, R. Rodrı´guez, J. Solomin, Finite Element anlysis of a quadratic eigenvalue problem arising in disspative acoustic, SIAM J. Numer. Anal. 38 (2000) 267–291.

[5] A. Bermúdez, R. Rodrı´guez, Finite element computation of the vibration modes of a fluid–solid system, Comput. Methods Appl. Mech. Eng. 119 (1994) 355–370.

[6] A. Bermúdez, R. Rodrı´guez, Modeling and numerical solution of elastoacoustic vibrations with interface damping, Int. J. Numer. Methods Eng. 46 (1999) 1763–1779.

[7] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991.

[8] H.C. Chen, R.L. Taylor, Vibration analysis of fluid–solid systems using a finite element displacement formulation, Int. J. Numer. Methods Eng. 29 (1990) 683–698.

[9] S.H. Chou, S. Tang, ConservativeP1 conforming and nonconforming Galerkin FEMs: effective flux evaluation via a non-mixed method approach, SIAM J. Numer. Anal. 38 (2000) 660–680.

[10] H.Y. Fan, W.W. Lin, P. Van Dooren, Normwise scaling of second order polynomial matrices, SIAM J. Matrix Anal. Appl. 26 (2004) 252–256. [11] L. Gastaldi, Mixed finite element methods in fluid structure systems, Numer. Math. 74 (1996) 153–176.

[12] M. Hamdi, Y. Ouset, G. Verchery, A displacement method for the analysis of vibrations of coupled fluid–structure systems, Int. J. Numer. Methods Eng. 13 (1978) 139–150.

[13] V. Hernandez, JE Roman, A. Tomas, V. Vidal. Krylov–Schur Methods in SLEPc, Technical report, CSE-2008-13, University of California, Davis, USA, 2008. Available online<http://www.grycap.upv.es/slepc>.

[14] T.-M. Hwang, W.-W. Lin, J.-L. Liu, W. Wang, Jacobi–Davidson methods for cubic eigenvalue problems, Numer. Linear Algebr. Appl. 12 (2005) 605–624.

[15] V. Kehr-Kandille, R. Ohayon, Elastoacoustic damped vibrations. Finite element and modal reduction methods, in: O.C. Zienkiewicz, P. Ladev‘eze (Eds.), New Advances in Computational Structural Mechanics, Elsevier, Amsterdam, 1992, pp. 321–334.

[16] R.B. Lehoucq, D.C. Sorensen, C. Yang, ARPACK USERS GUIDE: Solution of Large Scale Eigenvalue Problems With Implicitely Restarted Arnoldi Methods, SIAM, Philadelphia, 1998.

[18] P.A. Raviart, J.M. Thomas, A mixed finite element method for second order elliptic problems, Mathematical Aspects of Finite Element Methods, Lecture Notes in Math, vol. 606, Springer-Verlag, Berlin, Heidelberg, 1977, pp. 292–315.

[19] R. Rodrı´guez, J. Solomin, The order of convergence of eigenfrequencies in finite element approximations of fluid-structure interaction problems, Math. Comput. 65 (1996) 1463–1475.

[20] A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674–689.

[21] G.L.G. Sleijpen, A.G.L. Booten, D.R. Fokkema, H.A. van der Vorst, Jacobi-Davidson type methods for generalized eigenproblems and polynomial eigenproblems, BIT 36 (1996) 595–633.

[22] D.C. Sorensen, Implicit application of polynomial filters in ak-step Arnoldi method, SIAM J. Matrix Anal. Appl. 13 (1992) 357–385. [23] G.W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM J. Matrix Anal. Appl. 23 (2001) 601–614.

[24] G.W. Stewart, Addendum to A Krylov–Schur algorithm for large eigenproblems, SIAM J. Matrix Anal. Appl. 24 (2002) 599–601.

[25] Y. Su, Z. Bai. Solving Rational Eigenvalue Problems via Linearization, Technical report, 2008-13, University of California, Davis, USA, 2008. Available online<http://www.cs.ucdavis.edu/research/tech-reports/2008/CSE-2008-13.pdf>.

[26] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev. 43 (2001) 235–286. [27] H. Voss, An Arnoldi method for nonlinear eigenvalue problems, BIT 44 (2004) 387–401.

[28] H. Voss, Iterative projection methods for computing relevant energy states of a quantum dot, J. Comput. Phys. 217 (2006) 824–833.

[29] X. Wang, K.J. Bathe, Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems, Int. J. Numer. Methods Eng. 40 (1997) 2001–2017.