Efficient Arnoldi-Type Algorithms for Rational Eigenvalue Problems
Arising in Fluid-Solid Systems
So-Hsiang Choua,∗, Tsung-Ming Huangb, Wei-Qiang Huangc, Wen-Wei Lind
aDepartment of Mathematics and Statistics, Bowling Green State University, Bowling Green, OH, 43403-0221.
email: [email protected].
bDepartment of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan. E-mail: [email protected].
cDepartment of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan. E-mail: [email protected].
dDepartment of Mathematics, National Taiwan University, Taipei 106, Taiwan. E-mail: [email protected].
Abstract
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 displacement field. The trimmed linearization method is used to linearize the associated rational eigenvalue problem into a generalized eigenvalue problem (GEP) of the form Ax = λBx. For solving the GEP we apply Arnoldi algorithm to two different types of single matrices B−1A and AB−1. Numerical accuracy shows that the application of Arnoldi on AB−1 is better than that on B−1A.
Key words: fluid-structure interaction, finite elements, rational eigenvalue problem, trimmed linearization, Arnoldi algorithm
1. Introduction
Efficient and correct computation of the damped vibration modes generated by an inviscid, com-pressible, 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 mathemati-cal 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 methods 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 dis-placement 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 automati-cally. A drawback lurking 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 circumventing this drawback can be found in [2, 8, 11, 12, 30], among others.
One of the discretizations we will be using in this paper is the edge-based or Raviart-Thomas finite elements for the displacement field, following [3, 5]. The main concerns in [3, 20] 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 Section 2 that 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(λ)p = 0, where R(λ) := λ2M + λ2
α+βλA + K is the rational λ–matrix with coefficient matrices M, A and K. While there is an extensive literature on QEPs problems [27], REPs are much less studied [26, 28, 29]. Although on the surface the rational R(λ)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 Section 3.
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 particular attention to identifying the dimension of the associated null space, which may influence performance of the numerical method introduced later. In Section 3, 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.6.1) and (3.6.2) are derived. The REP is trimmed-linearized into two types of three by three block SEPs (3.17.1) and (3.17.2). We then apply Arnoldi method with Krylov-Schur restarting described in Section 4 to the resulting SEPs. The important issue of residual error bound analysis is addressed here. In Section 5, 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 Section 6. Basically, we conclude that the pressure formulation has much fewer variables and so has a small eigenvalue problem to solve. Once a good method is designed to solve its associated rational eigenvalue problem, the pressure formulation should be the preferred choice.
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 domain Ω ⊂ Rn (n = 2 or 3) to be polyhedral, and the boundary ∂Ω = ΓA∪ ΓR, where the absorbing boundary ΓA is the union of all the different faces of Ω and is covered by damping material. The rigid boundary ΓRis the remaining part of Γ. An example of the setup is in Figure 1, 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]) ρ∂ 2U ∂t2 + ∇P = 0 in Ω (2.1) P = −ρc2divU in Ω (2.2) P = αU · n + β∂U ∂t · n on ΓA (2.3) U · n = 0 on ΓR. (2.4)
Here ρ is the fluid density, c for the acoustic speed, and the unit outer normal vector along ∂Ω is denoted by n. At the absorbing boundary (2.3) indicates that the pressure is balanced by the effects of the viscous damping (the β term) and the elastic behavior (the α term). We assume the coefficients α and β 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) = eλtu(x) and P (x, t) = eλt
p(x). This leads to a problem of finding λ ∈ C, u : Ω → Cn and p : Ω → C, (u, p) 6= (0, 0) such that
ρλ2u + ∇p = 0 in Ω (2.5)
p = −ρc2 divu in Ω (2.6)
p = (α + λβ)u · n on ΓA (2.7)
u · n = 0 on ΓR. (2.8)
The boundary condition (2.7) makes this eigenvalue problem nonlinear. For each damped vibration mode, ω := Imλ is the vibration angular frequency and Reλ the decay rate. In practice, we select a range of ω 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 := {v ∈ H(div, Ω) : v · n ∈ L2(∂Ω) and v · n = 0 on Γ R}.
Here we employ standard Sobolev spaces notation. For example, H(div, Ω) stands for the space of all L2vector functions v on Ω with L2integrable divergence.
Testing (2.5) by ¯v ∈ V and integrating by parts, we obtain a variational formulation of problem (2.5)–(2.8) involving only the displacement variable: Find λ ∈ C and u ∈ V, u 6= 0, such that
Z Ω ρc2 divu div¯v + Z ΓA αu · n¯v · n + λ Z ΓA βu · n¯v · n + λ2 Z Ω ρu · ¯v = 0 ∀ v ∈ V. (2.9) (We use the symbol ∀ to mean ‘for all’ throughout the paper.) This is a quadratic eigenvalue problem. Note that λ = 0 is an eigenvalue and the dimension of its eigenspace
K := {u ∈ V : div u = 0 in Ω and u = 0 on ∂Ω}
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 6= λ ∈ C and 0 6= u ∈ V is a solution of problem (2.9) then Reλ < 0.
Alternatively we can derive a variational formulation in terms of the pressure: Find λ ∈ C and p ∈ H1(Ω) such that Z Ω ∇p · ∇¯q + λ 2 α + λβ Z ΓA ρp¯q + λ 2 c2 Z Ω p¯q = 0 ∀ q ∈ H1(Ω). (2.10)
However, in this case the eigenvalue problem is rational, which is rarely studied compared with linear and quadratic eigenvalue problems. Note that in contrast to the displacement formulation, the eigenspace corresponding to λ = 0 is now one dimensional. Thus this formulation has a much smaller null space or kernel, which may be more stable and efficient to 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´udez 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, 19]). For simplicity we will consider only the two dimensional case. Let {Th} be a regular family of triangulations of Ω indexed by h, the maximum diameter of the elements. Let
Vh:= {vh∈ H(div, Ω) : vh|T ∈ P0n⊕ P0x ∀ T ∈ Th and vh· n = 0 on ΓR} ⊂ V
where n = 2 and Pk denotes the set of polynomials of degree at most k. Thus locally vh takes the form (a + sx, b + sy)T (T denotes transpose and x = (x, y)T). The discrete problem associated with (2.9) is : Find λ ∈ C and uh∈ Vh, uh6= 0, such that
Z Ω ρc2 divuhdiv¯vh+ Z ΓA αuh· n ¯vh· n + λ Z ΓA βuh· n ¯vh· n + λ2 Z Ω ρuh· ¯vh= 0, ∀ vh∈ 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= uh and λ = 0 in (2.11), we see that
divuh= 0 on Ω and uh· n = 0 on ∂Ω.
Since uh = (a + sx, b + sy)T on T ∈ Th, the divergence free condition implies that uh is a constant vector (a, b)T on T . By direct computation, we see that there exists a linear polynomial ψ
T such that ∂ψT
∂x = −b and ∂ψT
∂y = a.
Let n = (n1, n2)T be a unit normal to an edge e of T , so t = (−n2, n1)T is a unit tangent vector to e. We see that
uh· n = ∇ψT · t = ∂ψT
∂t .
So if an edge e is common to T1and T2 then in general ψT1 and ψT2 differ 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 ψT at that node. Here T are all triangles sharing Nj as the common node. We then spread this defining process outward to all Ω using the induced values on other nodes. Consequently, Ψ is continuous piecewise linear over Ω. Let ∇⊥:= (−∂y∂ ,∂x∂ )T and define
∇⊥Sh:= {∇⊥Ψh: Ψh is continuous piecewise linear and vanishes on the boundary}.
Thus we have just shown the zero eigenspace E0 is contained ∇⊥Shand the opposite inclusion is also easily checked. Hence
We now find the dimension of ∇⊥Sh. Let N be the number of interior nodes and let Ψj, j = 1, . . . , N , be the nodal basis functions such that Ψj(Nk) = δkj. The linear independence of Ψj’s is preserved by the perp-gradient operation. In fact, supposePN
j=1cj∇ ⊥Ψ
j= 0. Then this impliesPNj=1cjΨj = c for some constant c. Hence cj = c by the condition Ψj(Nk) = δkj. Consequently, c(PjΨj− 1) = 0. But we knowPN
j=1Ψj 6= 1 due to the vanishing boundary condition. Thus cj = c = 0 and we conclude that the dimension of the zero eigenspace dim E0= dim ∇⊥Shequals the number of interior nodes in the mesh.
Define
Hh:= {ph∈ L2(Ω) : ph|T ∈ P0 ∀ T ∈ Th} ⊂ H1(Ω).
The alternative discrete problem in terms of the approximate pressure field is: Find λ ∈ C and ph∈ Hhsuch that Z Ω ∇ph· ∇¯qh+ λ2 α + λβ Z ΓA ρphq¯h+ λ2 c2 Z Ω phq¯h= 0 ∀ qh∈ Hh. (2.12) Letting qh = ph and λ = 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.
Again we see that the pressure formulation has a much smaller null space than the displacement formulation. Also the number of unknowns is also much smaller. Thus the pressure formulation is a very good alternative, once we show that its associated eigenvalue problem can be efficiently solved, as we show in the remaining sections.
Remark 2.1. Suppose an eigenpair (λ, ph), λ 6= 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 again so as to main 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, ∇·uh= 0, which certainly does not approximate (2.6). Fortunately, a general principle for such a problem (recovery of uh from the pressure approximation ph) has been provided in [9] where one can obtain an accurate uh in the Raviart–Thomas space by a simple evaluation formula which is a modification of the above naive formula. The detail of this approach will be given elsewhere.
3. Linearization of nonlinear eigenvalue problems
In this section we 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 {φj}nj=11 denote the cardinal basis of Vh, so that on the edge ej, φj has the unit normal flux and zero normal flux on the remaining n1− 1 edges. That is,
R
eiφj · nidσ = δij. For uh ∈ Vh, we write uh = Pn1
j=1ujφj and denote u = [uT1, · · · , uTn1]
T. Note that the unknown vector u contains normal fluxes in its components. Then, the discrete problem (2.11) can be expressed as the following QEP:
λ2Muu + (α + λβ)Auu + Kuu = 0, (3.1)
where Mu ≡ [Miju] and Ku ≡ [Kiju] are mass and stiffness matrices, respectively, and Au ≡ [Auij] is used to describe the effect of the absorbing wall. Here
Miju = Z Ω ρφi· ¯φj, Kiju = Z Ω ρc2 divφi div ¯φj, Auij = Z ΓA φi· n ¯φj· n, (3.2)
for i, j = 1, . . . , n1. For this problem, we are only interested in the eigenvalues that are located in the interior of the spectrum. Suppose that the eigenvalues near σ are of interest. Accordingly, the QEP (3.1) is shifted into µ2Mfu+ µ eDu+ eKu u = 0 (3.3) with µ = λ − σ and f Mu= Mu, eDu= 2σMu+ βAu, eKu= σ2Mu+ (α + σβ)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 problems [22]. On the other hand, it is more common to transform or linearize (3.3) into a standard eigenvalue problem [17]. In this paper we linearize (3.3) into the GEP
" 0 − fMu I − eDu # v u = 1 µ I 0 0 Keu v u . (3.5)
The matrix eKu in (3.5) is nonsingular because the shift value σ is not an eigenvalue of (3.1). Thus, the GEP (3.5) can be transformed into the following two types of the SEP:
" 0 − fMu e K−1 u − eKu−1Deu # v u = 1 µ v u (3.6.1) and " 0 − fMuKeu−1 I − eDuKeu−1 # v w = 1 µ v w , v = eKu−1w. (3.6.2)
The standard Arnoldi method can then be applied to solve (3.6.1) and (3.6.2). We will give details in Section 4.
3.2. Trimmed linearization method for rational eigenvalue problem
Let {ψj}nj=12 be a nodal basis of Hh. For ph ∈ Hh, we write ph = P n2
j=1pjψj and denote p = [p1, · · · , pn2]
T. Then, the discrete problem (2.12) can be written as the following REP:
R(λ)p ≡ λ 2 c2Mp+ Kp+ λ2 λβ + αAp p = 0, (3.7)
where Mp ≡ [Mijp] and Kp ≡ [Kijp] are mass and stiffness matrices, respectively, and Ap ≡ [Apij] describes the effect of the absorbing wall. Here,
Mijp = Z Ω ψiψ¯j, Kijp = Z Ω ∇ψi· ∇ ¯ψj, Apij = Z ΓA ρψiψ¯j (3.8) for i, j = 1, . . . , n2.
To solve REP (3.7), one approach is to multiply equation (3.7) by the scalar λβ + α 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.7) as nonlinear eigenvalue problem and solve it by a nonlinear eigen-solver, such as Newton’s method, nonlinear Arnoldi method, or nonlinear Jacobi-Davidson method
[21, 28, 29]. Recently, a trimmed inearization is proposed in [26] which linearizes (3.7) into a GEP so that the standard Arnoldi method can be applied. We introduce the trimmed linearization below.
The rational λ-matrix R(λ) in (3.7) can be rewritten as R(λ) = (λ − σ + σ) 2 c2 Mp+ Kp+ (λ − σ + σ)2 (λ − σ + σ)β + αAp = (λ − σ) 2+ 2(λ − σ)σ + σ2 c2 Mp+ Kp+ (λ − σ)2+ 2(λ − σ)σ + σ2 (λ − σ)β + σβ + α Ap = µ2 1 c2Mp + µ 2σ c2Mp + σ 2 c2Mp+ Kp +µ 2+ 2µσ + σ2 µβ + σβ + α Ap = µ2 1 c2Mp +1 µ 2σ c2Mp + 1 µ2 σ2 c2Mp+ Kp +1 + 2σ/µ + σ 2/µ2 µβ + (σβ + α) Ap (3.9) with µ = λ − σ, where σ is a given shift value. Setting ν = 1/µ, the rational term in (3.9) can be simplified into the following form by applying the long division
σ2ν2+ 2σν + 1 β/ν + (σβ + α) = σ 2 σβ + αν 2+σ 2β + 2σα (σβ + α)2ν + α2 (σβ + α)3 − (σβ + α)4 α2β ν + (σβ + α)3 α2 −1 . This implies that
R(λ) = 1 ν2 ν2 σ 2 c2Mp+ Kp+ σ2 σβ + αAp + ν 2σ c2Mp+ σ2β + 2σα (σβ + α)2Ap + 1 c2Mp+ α2 (σβ + α)3Ap − (σβ + α) 4 α2β ν + (σβ + α)3 α2 −1 Ap # = µ2Mfp+ µ eDp+ eKp− µ2 (σβ + α)4 α2βµ + (σβ + α)3 α2 −1 LpRpT where f Mp = 1 c2Mp+ α2 (σβ + α)3Ap, (3.10) e Dp = 2σ c2Mp+ σ2β + 2σα (σβ + α)2Ap, (3.11) e Kp = σ2 c2Mp+ Kp+ σ2 σβ + αAp, (3.12)
where LpRTp = Apis the full-rank decomposition of Apwith Lp, Rp∈ Rn2×m. Introducing an auxiliary vector q = (σβ + α) 4 α2βµ + (σβ + α)3 α2 −1 RTpp (3.13)
the REP in (3.7) can be reformulated as
µ2Mfp+ µ eDp+ eKp
Using (3.13) and (3.14), we get the GEP 0 − fMp Lp In2 − eDp 0 0 −RT p ϑIm x = 1 µ In2 0 0 0 Kep 0 0 0 %Im x, (3.15) where x = (1µKep+ eDp)p p q , % = − (σβ + α)4 α2β , ϑ = (σβ + α)3 α2 . (3.16)
The matrix eKp in (3.12) is nonsingular because the shift value σ is not an eigenvalue of (3.7). Con-sequently, as in (3.6.1) and (3.6.2), the GEP (3.15) can then be transformed into the following two types of the SEP:
0 − fMp Lp e Kp−1 − eKp−1Dep 0 0 −%−1RT p %−1ϑIm x = 1 µx, (3.17.1) and 0 − fMpKep−1 %−1Lp In2 − eDpKe −1 p 0 0 −RT pKep−1 %−1ϑIm y = 1 µy, y = In2 0 0 0 Kep 0 0 0 %Im x. (3.17.2)
For convenience, the SEPs in (3.6.1) (or (3.17.1)) and in (3.6.2) (or (3.17.2)) can be referred as to a type-1 and type-2 SEPs, respectively.
3.3. Error analysis
In this subsection, we will discuss residuals of QEP (3.1) and REP (3.7) by using linearizations of (3.6) and (3.17), respectively.
We first derive residual bounds of approximate eigenpairs for QEP (3.1) by using linearizations (3.6.1) and (3.6.2), respectively. Let (µ1, [vT1, uT1]T) be an approximate eigenpair of (3.6.1) and [fT
11, f12T]T be the associated residual vector. That is, f11 f12 = " 0 − fMu e Ku−1 − eKu−1Deu # v1 u1 − 1 µ1 v1 u1 = 1 µ1 " −v1− µ1Mfuu1 e Ku−1(µ1v1− µ1Deuu1− eKuu1) # . It follows that µ21Mfuu1+ µ eDuu1+ eKuu1 = µ1(−v1− µ1f11) + µ1v1− µ1Keuf12 = −µ1Keuf12− µ21f11. Therefore, kµ2 1Mfuu1+ µ eDuu1+ eKuu1k ku1k ≤ |µ1|k eKukkf12k + |µ1| 2kf 11k ku1k . (3.18)
On the other hand, let (µ2, [vT2, wT2]T) be an approximate eigenpair of (3.6.2) and [f21T, f22T]T be the associated residual vector. That is,
f21 f22 = " 0 − fMuKeu−1 I − eDuKeu−1 # v2 w2 − 1 µ2 v2 w2 = " − fMuKeu−1w2−µ12v2 v2− eDuKeu−1w2−µ1 2w2 # . It follows that v2 = −µ2MfuKeu−1w2− µ2f21, −µ2f22− µ22f21 = µ22MfuKeu−1w2+ µ2DeuKeu−1w2+ w2. (3.19) Letting u2= eKu−1w2, we have kµ2 2Mfuu2+ µ2Deuu2+ eKuu2k ku2k ≤ |µ2|kf22k + |µ2| 2kf 21k ku2k . (3.20)
Now, we derive residual bounds of approximate eigenpairs for REP (3.7) by using linearizations of (3.17.1) and (3.17.2), respectively. Let (µ1, [sT1, pT1, qT1]T) be an approximate eigenpair of (3.17.1) and [gT
11, g12T, gT13]T be the associated residual vector. That is, g11 g12 g13 = 0 − fMp Lp e Kp−1 − eKp−1Dep 0 0 −%−1RT p %−1ϑIm s1 p1 q1 − 1 µ1 s1 p1 q1 .
This implies that
s1 = −µ1Mfpp1+ µ1Lpq1− µ1g11, (3.21) g12 = Kep−1s1− eKp−1Depp1− 1 µ1 p1, (3.22) q1 = µ1%−1ϑ − 1 −1 µ1 g13+ %−1RTpp1 . (3.23) Substituting (3.23) into (3.21) and using the definitions of % and ϑ in (3.16), s1can be represented by s1 = −µ1Mfpp1+ µ21 µ1%−1ϑ − 1 −1 Lpg13+ %−1LpRTpp1 − µ1g11 = −µ1Mfpp1+ µ1 (σβ + α)3 α2 + (σβ + α)4 α2βµ 1 −1 LpRTpp1− µ1 β σβ + α + 1 µ1 −1 Lpg13− µ1g11. (3.24) Substituting (3.24) into (3.22), we get
R(σ + µ1)p1 = µ21Mfpp1+ µ1Depp1+ eKpp1− µ21 (σβ + α)3 α2 + (σβ + α)4 α2βµ 1 −1 LpRTpp1 = −µ2 1g11− µ1Kepg12− µ21 β σβ + α+ 1 µ1 −1 Lpg13 which implies that
kR(σ + µ1)p1k kp1k ≤ 1 kp1k ( |µ1|2kg11k + |µ1|k eKpk kg12k + µ21 β σβ + α+ 1 µ1 −1 kLpk kg13k ) . (3.25)
On the other hand, let (µ2, [sT2, tT2, qT2]T) be an approximate eigenpair of (3.17.2) and [gT21, gT22, g23T]T be the associated residual vector. That is,
g21 g22 g23 = 0 − fMpKep−1 %−1Lp In2 − eDpKe −1 p 0 0 −RT e Kp−1 %−1ϑIm s2 t2 q2 − 1 µ2 s2 t2 q2 .
This implies that
g21 = − fMpKep−1t2+ %−1Lpq2− 1 µ2 s2, (3.26) s2 = DepKep−1t2+ 1 µ2 t2+ g22, (3.27) q2 = %−1ϑ − 1 µ2 −1 RTpKep−1t2+ g23 . (3.28)
Substituting (3.27) and (3.28) into (3.26), we have
µ22g21 = −µ22MfpKep−1t2+ µ22 (σβ + α)3 α2 + (σβ + α)4 α2βµ 2 −1 LpRTpKep−1t2+ Lpg23 −µ2DepKep−1t2− t2− µ2g22. Letting p2= eKp−1t2, we get R(σ + µ2)p2 = µ22Mfpp2+ µ2Depp2+ eKpp2− µ22 (σβ + α)3 α2 + (σβ + α)4 α2βµ 2 −1 LpRTpp2 = −µ2 2g21− µ2g22+ µ22 α2β (σβ + α)4 β σβ + α+ 1 µ2 −1 Lpg23. Hence, kR(σ + µ2)p2k kp2k ≤ 1 kp2k ( |µ2|2kg21k + |µ2|kg22k + µ22 α 2β (σβ + α)4 β σβ + α+ 1 µ2 −1 kLpkkg23k ) . (3.29) Remark 3.1. Suppose in (2.9) we adopt the following data as in [4]. That is we set ρ = 1kg/m3, c = 340 m/s, α = 5 × 104 N/m3, and β = 200 Ns/m3. In addition let us choose σ = −25 + 600πı as the shift value. Then
(i) From (3.4), we see that k eKuk∞≈ (8 + 6 √
2)ρc2≈ 1.9 × 106. From (3.18) and (3.20), we conclude that the upper bound of the residual of the approximate eigenpair (σ + µ1, u1) for (3.1) by using (3.6.1) is larger than that of the approximate eigenpair (σ + µ2, u2) for (3.1) by using (3.6.2). (ii) From (3.12) we see that k eKpk∞ ≈ 8. If the eigenvalue λ is one of the desired eigenvalues in
Figure 3, then with µ = λ − σ we have
4 × 107< µ2 β σβ + α+ 1 µ −1 < 3.1 × 1010, 0.001 < µ2 α 2β (σβ + α)4 β σβ + α+ 1 µ −1 < 0.8.
Clearly, from (3.25) and (3.29) we conclude that the upper bound of the residual of the approxi-mate eigenpair (σ + µ1, p1) for REP (3.7) by using (3.17.1) is larger than that of (σ + µ2, p2) for (3.7) by using (3.17.2).
4. Arnoldi method with Krylov-Schur restarting
Arnoldi method is the most popular method for solving large sparse linear eigenvalue problems such as (3.6) and (3.17). In Arnoldi process, an orthonormal matrix Vk+1 is generated to satisfy
AVk = VkHk+ hk+1,kvk+1eTk, (4.1)
where Hk ∈ Rk×k is upper Hessenberg and A is the coefficient matrix of the eigenvalue problem in (3.6) or (3.17). The Arnoldi process is restarted once the index k becomes large. 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 [18, 23] in which the Arnoldi process is combined with the implicitly shifted QR algorithm. The package ARPACK [16] includes a very successful implementation of the implicitly restarted Arnoldi algorithm. It has been used by many engineering fields and remains 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, 24, 25] 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 Krylov-Schur restarting scheme as follows. Let Hk= Uj U` Tj Tf 0 T` Uj∗ U`∗ (4.2) be a Schur decomposition of Hk where Tj and T` are upper triangular, and the eigenvalues of Tj are of interest. Substituting (4.2) into (4.1), we see that
A(Vk Uj U` ) = (Vk Uj U` ) Tj Tf 0 T` + hk+1,kvk+1(eTk Uj U` ), which implies that
A ˜Vj= ˜VjTj+ ˜vj+1t∗j, (4.3)
where ˜Vj≡ VkUj, ˜vj+1= vk+1 and t∗j ≡ hk+1,keTkUj. Let Q1be a Householder matrix with t∗jQ1= τ eTj.
Then (4.3) can be rewritten as
A( ˜VjQ1) = ( ˜VjQ1)(Q1∗TjQ1) + τ ˜vj+1eTj. (4.4) The matrix Q∗1TjQ1 can be reduced to a new Hessenberg matrix ˜Hj by using Householder matrices Qi for i = 2, . . . , j − 1 with
Multiplying (4.4) by Qi for i = 2, . . . , j − 1, a new Arnoldi decomposition of order j A ˜Vj= ˜VjH˜j+ τ ˜vj+1eTj
is obtained where ˜Vj:= ˜VjQ1· · · Qj−1 and 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.
Applying Arnoldi method to solve QEP (3.1) and REP (3.7) are summarized in Algorithms 2 and 3, respectively.
Algorithm 1 Arnoldi method with Krylov-Schur restarting for solving Ax = λx Input: Coefficient matrix A.
Output: The desired j eigenpairs.
1: Build an initial Arnoldi decomposition of order j + `:
AVj+`= Vj+`Hj+`+ hj+`+1,j+`vj+`+1eTj+`.
2: repeat
3: Test for convergence.
4: Compute Schur decomposition of Hj+`:
Hj+`= Uj U` Tj Tf 0 T` Uj∗ U`∗ , where the eigenvalues of Tj are of interest.
5: Set Vj := Vj+`Uj, vj+1:= vj+`+1and t∗j := hj+`+1,j+`eTj+`Uj.
6: Compute Householder transformation Q1such that t∗jQ1= τ eTj.
7: Reduce Q∗1TjQ1 to a new Hessenberg matrix Hj by using Householder transformations Qi for i = 2, . . . , j − 1.
8: Set Vj := VjQ1· · · Qj−1and hj+1,j= τ to get the new Arnoldi decomposition with order j: AVj = VjHj+ hj+1,jvj+1eTj. (4.5)
9: Extend (4.5) to a new Arnoldi decomposition of order j + `.
10: until desired j eigenpairs are convergent
Algorithm 2 Arnoldi method for solving QEP in (3.1)
Input: Coefficient matrices Mu, Du and Ku, parameters c, α and β and a shift value σ. Output: The desired eigenpairs (λi, ui) for i = 1, . . . , r.
1: Construct matrices fMu, eDuand eKu defined in (3.4).
2: Apply Algorithm 1 to compute eigenpairs (µi, [uTi, vTi ]T), for i = 1, . . . , r, of (3.6.1) or (3.6.2).
3: Set λi = σ + µ−1i for i = 1, . . . , r.
4: if (3.6.2) is solved then
5: Solve linear system eKuxi= ui and reset ui:= xi for i = 1, . . . , r.
Algorithm 3 Arnoldi method for solving REP in (3.7)
Input: Coefficient matrices Mp, Kp and Ap, parameters c, α and β and shift value σ. Output: The desired eigenpairs (λi, pi) for i = 1, . . . , r.
1: Construct matrices fMp, eDp and eKp defined in (3.12), (3.11) and (3.10), respectively.
2: Compute the full-rank decomposition of Ap: LpRTp = Ap.
3: Apply Algorithm 1 to compute eigenpairs (µi, [sTi, p T i , q T i] T), for i = 1, . . . , r, of (3.17.1) or (3.17.2). 4: Set λi = σ + µ−1i for i = 1, . . . , r. 5: if (3.17.2) is solved then
6: Solve linear system eKpxi= pi and reset pi:= xi for i = 1, . . . , r.
7: end if
5. Numerical Results
We conduct numerical experiments to evaluate performance and accuracy of 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 type-1 SEP as in (3.6.1). • Q2: Applying Algorithm 2 to solve the QEP (3.1) with type-2 SEP as in (3.6.2). • R1: Applying Algorithm 3 to solve the REP (3.7) with type-1 SEP as in (3.17.1). • R2: Applying Algorithm 3 to solve the REP (3.7) with type-2 SEP as in (3.17.2).
All computations are carried out in MATLAB 2009a on a HP workstation with an Intel Quad-Core Xeon X5570 2.93GHz and 72 GB main memory, using IEEE double-precision floating-point arithmetic. In Algorithms 2 and 3, the linear system is solved by LU factorization with the shift value σ = −25 + 600πı and λ1, . . . , λ10in Figure 3 are the desired eigenvalues. In Algorithm 1, j + ` is set to be 40 and the residual tolerance ε in the standard stopping criterion is set to 10−5.
Let (λ, x) be an approximate eigenpair of QEP in (3.1) or REP in (3.7). We define the associated relative residual of (λ, x) for (3.1) and (3.7) by
kλ2M ux + λDux + KuxkF (|λ2|kM ukF+ |λ|kDukF + kKukF)kxkF and kλ2 c2Mpx + Kpx + λ 2 λβ+αApxkF (|λc22|kMpkF + kKpkF+ | λ 2 λβ+α|kApkF)kxkF , respectively.
Example 1. [4] We take the geometrical data: the domain Ω = [0m, 1m] × [−0.75m, 0m], ΓA = [0m, 1m] × {0m} given in Figure 1 and the following physical data: ρ = 1kg/m3, c = 340 m/s, α = 5 × 104 N/m3, and β = 200 Ns/m3.
Figure 3 shows that the desired eigenvalues λ1, . . . , λ10with the lowest positive vibration frequen-cies satisfying 0 < Im(λi)
2π < 600Hz. The rectangular domain Ω is uniformly partitioned into M by N rectangles and each rectangle is further refined into two triangles, see Figure 2. The dimensions of coefficient matrices in QEP (3.1) and REP (3.7) are (3M − 1) × N and (M + 1) × (N + 1), respectively.
0.75 m b= 1.00 m a= Ω absorbing wall ΓA R Γ
n
rigid wallsFigure 1: Fluid in a cavity with one absorbing wall.
(M, N ) Matrix size Matrix size Eigenvalue Rate of convergence N × (3M − 1) (M + 1) × (N + 1)
( 48, 36) 5, 148 1, 813 λ[i,1]
-( 96, 72) 20, 664 7, 081 λ[i,2]
-(192, 144) 82, 800 27, 985 λ[i,3] rate[i,1]
(384, 288) 331, 488 111, 265 λ[i,4] rate[i,2]
(768, 576) 1, 326, 528 443, 713 λ[i,5] rate[i,3]
Table 1: Dimension information regarding the computation of convergence rates.
−3500 −300 −250 −200 −150 −100 −50 0 500 1000 1500 2000 2500 3000 3500 4000 real part imaginary part λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10
Figure 3: The distribution of the ten target eigenvalues λ1, . . . , λ10.
We first demonstrate convergence rates of the two schemes while computing all desired eigenvalues λi for i = 1, . . . , 10 in Figure 3. To measure the convergence rate, we run the test over the five successively refined meshes (first column of Table 1) and then calculate the rates by
rate[i,j]= log2 |λ
[i,j]− λ[i,j+1]| |λ[i,j+1]− λ[i,j+2]|
, for i = 1, . . . , 10, j = 1, . . . , 3
where λ[i,j] for j = 1, . . . , 5 denote the approximate eigenvalues corresponding to λi obtained from the meshes described in Table 1. Figure 4 (i) and (ii) illustrate the convergence rates, rate[i,j] for i = 1, . . . , 10, j = 1, . . . , 3, for all 10 desired eigenvalues of QEP (3.1) and REP (3.7), respectively. From Figure 4, we see that all convergence rates are close to 2.
In [4], it has been proved that there are no spurious eigenmodes for the discretization based on Raviart-Thomas finite elements. In Figure 5, we numerically show that there are no spurious eigen-modes for the discretization in terms of pressure nodal finite elements. Twenty desired eigenvalues of QEP (3.1) and REP (3.7), respectively, are computed and shown in Figure 5. The desired eigenvalues of REP are in one-to-one correspondence to those of QEP. No spurious eigenmodes ever appear. (We computed 20 instead of 10 eigenvalues to be better confirmed.)
1.95 2 R a te o f C o n ve rg e n ce rate[1] rate [2] rate [3] l3 l1 l2 l4 l5 l6 l7 l8 l9 l10 2.05
(i) Raviart-Thomas finite element
1.95 2 R a te o f C o n ve rg e n ce rate[1] rate [2] rate [3] l l l3 l1 l2 l4 l5 l6 l7 l8 l9 l10 2.05
(ii) nodal finite element
Figure 4: The convergence rates of eigenvalues for Raviart-Thomas and nodal finite elements are shown in part (i) and (ii), respectively. −3500 −300 −250 −200 −150 −100 −50 0 1000 2000 3000 4000 5000 6000 Real part Imaginary part QEP REP
−3500 −300 −250 −200 −150 −100 −50 0 500 1000 1500 2000 2500 3000 3500 4000 real part imaginary part 2.6e−15 1.1e−16 4.1e−15 2.4e−16 8.2e−11 1.3e−14 1.8e−9 3.9e−11 1.6e−11 3.8e−5 Exact Unscaling (i) QEP −3500 −300 −250 −200 −150 −100 −50 0 500 1000 1500 2000 2500 3000 3500 4000 real part imaginary part 2.5e−11 1.9e−9 3.0e−12 1.0e−11 7.5e−8 7.0e−11 1.5e−11 4.6e−8 1.9e−10 1.1e−7 Exact Unscaling (ii) REP
Figure 6: The computed eigenpairs and the corresponding residual for QEP in (3.1) and REP in (3.7) without scaling are shown in part (i) and (ii), respectively. The notations “o” and “×” denote the exact and computed eigenvalues, respectively. The value near the notation “×” represents the relative residual of the computed eigenpair.
Normwise scaling of QEP: Balancing norms of coefficient matrices is an important issue [27] before solving a QEP of the form
P (µ)x ≡ (µ2P2+ µ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 b
P (ν)x ≡ (ν2Pb2+ ν bP1+ bP0)x = 0
with ν = µ/ζ, bP2= ζ2ηP2, bP1= ζηP1and bP0= ηP0, where ζ and η are scaling factors. Taking ζ and η as
ζ∗=pγ0/γ2, η∗= 2/(γ0+ γ1 p
γ0/γ2)
with γ2:= kP2k2, γ1:= kP1k2, γ0:= kP0k2, it is proved in [10] that the problem min
ζ,η max n
|k bP2k2− 1|, |k bP1k2− 1|, |k bP0k2− 1| o
achieves the optimum at ζ∗ and η∗. In our implementation, the values of γ
i, for i = 0, 1, 2 are computed by γ2= k fMukF, γ1 = k eDukF, γ0 = k eKukF and γ2= k fMpkF, γ1= k eDpkF, γ0 = k eKpkF for QEP (3.3) and REP (3.14), respectively. The computed eigenvalues and the corresponding relative residuals, obtained by Q2 and R2, respectively, for (3.1) and (3.7) without scaling are shown in (i) and (ii) of Figures 6, respectively. The results for (3.1) and (3.7) with scaling are shown in Figure 7. They demonstrate that the accuracy of eigenpairs is significantly improved by using the scaling technique. Theorem 1 shows that the dimension of the null space of QEP (3.1) is (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 “+” in Figure 8 to observe variation in the iteration number for Q1 and Q2 with different shift values. The integer pair (i, j) under the notation “+” in Figure 8 denotes the iteration numbers i and j for Q1 and Q2, respectively, with shift value at “+”. The results in Figure 8 shows that the iteration number needed decreases, as the shift value σ is chosen relatively far away from zero.
10-18 10-17 10-16 10-15 10-14 10-13 10-12 R e la tiv e R e si d u a l R 1 Q1 Q2 R 2 l3 l1 l2 l4 l5 l6 l7 l8 l9 l10
Figure 7: The relative residuals of computed eigenpairs, obtained by Q1, Q2 for QEP (3.1) and R1, R2 for REP (3.7) 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 ) ( 2, 2 ) ( 6, 4 ) ( 2, 2 ) ( 4, 4 ) ( 2, 2 ) Eigenvalues Shift values
Figure 8: The number of outer iterations of Q1 and Q2 with different shift values. “o” denotes desired eigenvalues λ1, . . . , λ10. “(i, j)” denotes the iteration numbers for Q1 and Q2 by i and j, respectively.
In this paragraph, we shall discuss the advantages of using the nodal pressure finite elements. We compare the results of solving (3.1) and (3.7) with scaling technique by using Algorithms 2 and 3, respectively. The relative residuals, produced by Q1, Q2, R1 and R2, of (3.1) and (3.7) with (M, N ) = (768, 576) are shown in Figure 7. Listed in Table 2 are the CPU times and the number of outer iterations (number of performing Krylov Schur restarting in Algorithm 1) for Q1, Q2, R1 and R2 with the meshes described in Table 1. The notation “#it” in Table 2 denotes the number of outer iterations. “TQ1”, “TQ2”, “TR1” and “TR2” denote the total CPU time for Q1, Q2, R1 and R2, respectively. We summarize the results as follows:
• CPU times for computing desired eigenpairs: In Algorithms 2 and 3, the computational cost for solving linear systems dominates that of matrix-product-vector and Krylov-Schur restarting. When the eigenvectors of SEPs (3.6.2) and (3.17.2) are computed, Q2 and R2 need to solve one extra linear system to get the eigenvectors of QEP (3.1) and REP (3.7), respectively. Therefore, the computational costs of Q2 and R2 are slightly more expensive than that of Q1 and R1, respectively, when they have the same number of outer iterations. This is verified by the results in Table 2. From columns six and eleven of Table 2, the extra CPU times TQ2− TQ1and TR2− TR1 of Q2 and R2 are only 1% to 7% of the total CPU times for Q1 and R1, respectively, when the mesh length is small enough.
• Accuracy of eigenpairs: From Remark 3.1, the upper bound of the residual for the approximate eigepairs of QEP (3.1) (or REP (3.7)) by using type-2 SEP in (3.6.2) (or (3.17.2)) is less than that by using type-1 SEP in (3.6.1) (or (3.17.1)). In the following, we compare the accuracy of approximate eigenpairs. On applying Q1 and Q2 to solve QEP (3.1), in Figure 7, we see that the relative residuals of eigenpairs corresponding to λ4 and λ5 computed by Q2 are improved about 1 significant digit than those by Q1. The other eigenpairs almost have the same accuracy. On applying R1 and R2 to solve REP (3.7), in Figure 7, the relative residuals of eigenpairs computed by R2 are improved about 2 to 4 significant digits than those by R1.
• Comparison R2 with Q2: From Table 2, the computational cost of Q2 (orR2) is slightly more expensive than that of Q1 (or R1). However, from Figure 7, we see that the relative residuals of the eigenpairs computed by Q2 (or R2) are improved about 1 (or 2 ∼ 4) significant digit than those by Q1 (or R1). Therefore, we favor applying Q2 and R2 to solve QEP (3.1) and REP (3.7), respectively. From column 12 of Table 2, we see that the computational cost for solving the REP (3.7) by R2 has 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.7) 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 Raviart-Thomas displacement finite elements for the discrete problem (2.11).
Based on the above discussion, the computational cost of solving REP (3.7) is obviously less than that of solving QEP (3.7). For solving REP (3.7), the accuracy of R2 is better than that of R1 when the stopping tolerance ε is equal to 10−5. In order to more precisely compare the accuracy of eigenpairs computed by R1 and R2, we show iteration numbers and relative residuals computed by R1 with different ε in Table 3. The matrix dimension is fixed to 443, 713. We do not improve the accuracy of λ2, . . . , λ9 when ε is changed from 10−5 to 10−12 and 10−15. The improvements only occur in the accuracy of λ1 and λ10. The cost for such improvements is the increase of the iteration number. Therefore, from Tables 2, 3 and Figure 7, it’s better to choose R2 over R1 when solving REP (3.7).
We now want to apply our methods to a more complicated configuration in which the absorbing wall are located on three sides.
Q1 Q2 TQ2 TQ1 R1 R2 TR2 TR1 TR2 TQ2
(M, N ) #it TQ1 #it TQ2 #it TR1 #it TR2
( 48, 36) 2 1.162 2 1.316 1.13 2 0.452 2 0.471 1.04 0.36
( 96, 72) 2 6.658 2 7.717 1.16 2 2.073 2 2.387 1.15 0.31
(192, 144) 2 51.85 2 55.27 1.07 2 14.16 2 14.95 1.06 0.27
(384, 288) 2 553.5 2 567.8 1.03 2 130.5 2 134.0 1.03 0.24
(768, 576) 2 7924 2 8152 1.03 2 1630 2 1645 1.01 0.20
Table 2: Iteration number and CPU times for Q1, Q2, R1 and R2.
R1 ε = 10−5 ε = 10−12 ε = 10−15
#it 2 3 4
RRes1 6.4 × 10−14 1.5 × 10−14 1.5 × 10−14 RRes10 3.8 × 10−13 1.5 × 10−14 1.5 × 10−14
Table 3: Iteration number and relative residual for different stopping tolerance ε. RRes1and RRes10denote the relative residuals of λ1and λ10, respectively.
Example 2. Here we use the same geometric data and physical data in Example 1 except that the absorbing wall is extended to one half of the rigid walls in the left and right boundaries. Γ = [0, 1] × {0} ∪ {0} × [−0.375, 0] ∪ {1} × [−0.375, 0].
In Example 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.7) 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 and efficiency of R1 and R2 for solving the associated REP. The computed eigenvalues λ1, . . . , λ10 with lowest positive vibration frequencies satisfying 0 < Im(λi)
2π < 600Hz are shown in Figure 9. The convergence rates for λ1, . . . , λ10 obtained from the meshes described in Table 1 is shown in (i) of Figure 10. All the convergence rates are also close to 2. As shown in Table 4, the computational cost of R2 is almost equal to that of R1. However, as shown in (ii) of Figure 10, the accuracy of the eigenpairs produced by R2 is better than that by 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
R1 R2 TR2
TR1
(M, N ) #it Time #it Time
( 48, 36) 2 0.503 2 0.522 1.04
( 96, 72) 2 2.256 2 2.552 1.13
(192, 144) 3 18.12 3 19.96 1.10
(384, 288) 3 161.2 3 172.6 1.07
(768, 576) 3 1830 3 1896 1.04
−6000 −500 −400 −300 −200 −100 0 500 1000 1500 2000 2500 3000 3500 4000 λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 Real part Imaginary part
Figure 9: The distribution of the ten target eigenvalues λ1, . . . , λ10for Example 2.
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 R a te o f C o n ve rg e n ce rate [1] rate [2] rate [3] l3 l1 l2 l4 l5 l6 l7 l8 l9 l10 (i) convergence rate
10-20 10-18 10-16 10-14 10-12 10-10 R e la tiv e R e si d u a l R 1 R 2 l3 l1 l2 l4 l5 l6 l7 l8 l9 l10 (ii) relative residual
Figure 10: The convergence rates of eigenvalues and relative residual of the computed eigenpairs for rational eigenvalue problem with (M, N ) = (768, 576) in (i) and (ii), respectively.
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.7). 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 pressure nodal finite elements and the CPU times for solving the corresponding REP (3.7) 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.7) into four different types of SEPs which can be solved by Q1/Q2 as well as R1/R2. Numerical accuracy shows that Q2 and R2 algorithms are better than Q1 and R1, respectively.
Acknowledgments
This work is partially supported by the National Science Council and the National Center for Theoretical Sciences of Taiwan.
References
[1] Z. Bai and Y. Su. A second-order Arnoldi method for the solution of the quadratic eigenvalue problem, SIAM J. Matrix Anal., 26(3):640–659, 2005.
[2] K. J. Bathe, C. Nitikitpaiboon, and X. Wang, A mixed displacement-based finite element for-mulation for acoustic fluid-structure interaction, Comput. Struct., 56, pp. 225–237, 1995. [3] A. Berm´udez, R. Dur´an, M. A. Muschietti, R. Rodr´ıguez, and J. Solomin, Finite element
vi-bration analysis of fluid-solid systems without spurious modes, SIAM J. Numer. Anal., 32, pp. 1280–1295, 1995.
[4] A. Berm´udez, R. G. Dur´an, R. Rodr´ıguez, and J. Solomin, Finite Element anlysis of a quadratic eigenvalue problem arising in disspative acoustic, SIAM J. Numer. Anal. 38, pp. 267–291, 2000. [5] A. Berm´udez and R. Rodr´ıguez, Finite element computation of the vibration modes of a
fluid-solid system, Comput. Meth. Appl. Mech. Eng., 119, pp. 355–370, 1994.
[6] A. Berm´udez and R. Rodr´ıguez, Modeling and numerical solution of elastoacoustic vibrations with interface damping, Int. J. Numer. Methods Eng., 46, pp. 1763–1779, 1999.
[7] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991.
[8] H. C. Chen and R. L. Taylor, Vibration analysis of fluid-solid systems using a finite element displacement formulation, Int. J. Numer. Methods Eng., 29, pp. 683–698, 1990.
[9] S. H. Chou and S. Tang, Conservative P 1 conforming and nonconforming Galerkin FEMs: effective flux evaluation via a non–mixed method approach, SIAM J. Numer. Anal., 38, pp. 660–680, 2000.
[10] H.Y. Fan, W.W. Lin, and P. Van Dooren, Normwise scaling of second order polynomial matrices. SIAM J. Matrix Anal. Appl., 26:252–256, 2004.
[11] L. Gastaldi, Mixed finite element methods in fluid structure systems, Numer. Math., 74, pp. 153–176, 1996.
[12] M. Hamdi, Y. Ouset, and G. Verchery, A displacement method for the analysis of vibrations of coupled fluid-structure systems, Int. J. Numer. Methods Eng., 13, pp. 139–150, 1978.
[13] V. Hernandez, JE Roman, A. Tomas, and V. Vidal. Krylov-Schur methods in SLEPc, Techni-cal report, Tech. Rep. CSE-2008-13, University of California, Davis, USA, 2008. Available at http://www. grycap. upv. es/slepc.
[14] T.-M. Hwang, W.-W. Lin, J.-L. Liu, and W. Wang. Jacobi–Davidson methods for cubic eigen-value problems. Numer. Linear Algebr. Appl., 12:605–624, 2005.
[15] V. Kehr-Kandille and R. Ohayon, Elastoacoustic damped vibrations. Finite element and modal reduction methods, in New Advances in Computational Structural Mechanics, O. C. Zienkiewicz and P. Ladev‘eze, eds., Elsevier, Amsterdam, pp. 321–334, 1992.
[16] R.B. Lehoucq, D.C. Sorensen, and C. Yang. ARPACK USERS GUIDE: Solution of large scale eigenvalue problems with implicitely restarted Arnoldi methods. SIAM, Philadelphia, 1998. [17] T. Li, Y.-T. Li, and W.-W. Lin. A novel restarted and refined Arnoldi-like method for computing
a few smallest eigenvalues of quadratic eigenvalue problems. Preprint, 2010.
[18] R.B. Morgan. On restarting the Arnoldi method for large non-symmetric eigenvalue problems. Math. Comput., 65(215):1213–1230, 1996.
[19] P. A. Raviart and J. M. Thomas, A mixed finite element method for second order elliptic problems, in Mathematical Aspects of Finite Element Methods, Lecture Notes in Math. 606, Springer-Verlag, Berlin, Heidelberg, pp. 292–315, 1977.
[20] R. Rodr´ıguez and J. Solomin, The order of convergence of eigenfrequencies in finite element approximations of fluid-structure interaction problems, Math. Comput., 65, pp. 1463–1475, 1996. [21] A. Ruhe. Algorithms for the nonlinear eigenvalue problem. SIAM J. Numer. Anal., 10:674–689,
1973.
[22] G.L.G. Sleijpen, A.G.L. Booten, D.R. Fokkema, and H.A. van der Vorst. Jacobi-Davidson type methods for generalized eigenproblems and polynomial eigenproblems. BIT, 36:595–633, 1996. [23] D. C. Sorensen. Implicit application of polynomial filters in a k-step Arnoldi method. SIAM J.
Matrix Anal. Appl., 13:357–385, 1992.
[24] G.W. Stewart. A Krylov–Schur algorithm for large eigenproblems. SIAM J. Matrix Anal. Appl., 23:601–614, 2001.
[25] G.W. Stewart. Addendum to “A Krylov–Schur algorithm for large eigenproblems”. SIAM J. Matrix Anal. Appl., 24:599–601, 2002.
[26] Y. Su and Z. Bai. Solving rational eigenvalue problems via linearization. Tech-nical report, 2008-13, University of California, Davis, USA, 2008. Available at http://www.cs.ucdavis.edu/research/tech-reports/2008/CSE-2008-13.pdf.
[27] F. Tisseur and K. Meerbergen. The quadratic eigenvalue problem. SIAM Review, 43:235–286, 2001.
[29] H. Voss. Iterative projection methods for computing relevant energy states of a quantum dot. J. Comp. Physics, 217:824–833, 2006.
[30] X. Wang and K. J. Bathe, Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems, Int. J. Numer. Methods Eng., 40, pp. 2001–2017, 1997.
