應
用
數
學
系
博
士
論
文
二次暨有理特徵值問題中高效能 Arnoldi 型態演算法
Efficient Arnoldi-Type Algorithms for
Quadratic and Rational Eigenvalue Problems
研 究 生：黃韋強
二次暨有理特徵值問題中高效能 Arnoldi 型態演算法
Efficient Arnoldi-Type Algorithms for
Quadratic and Rational Eigenvalue Problems
研 究 生：黃韋強
Student：Wei-Qiang Huang
指導教授：林文偉
Advisor：Wen-Wei Lin
國 立 交 通 大 學 應 用 數 學 系
博 士 論 文
A Thesis
Submitted to Department of Applied Mathematics College of Science
National Chiao Tung University in partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
in
Applied Mathematics August 2012
Hsinchu, Taiwan, Republic of China
二次暨有理特徵值問題中高效能 Arnoldi 型態演算法
學生：黃韋強 指導教授：林文偉 教授
國立交通大學應用數學系博士班
摘
要
本論文探討求解二次特徵值問題及有理特徵值之高效能 Arnoldi 型態演算法。其研究
主題可分為兩部分：（一）流固系統中非線性特徵值問題之 Arnoldi 型態演算法之比較；
（二）求解二次特徵值問題中的半正交廣義 Arnoldi 法。
我們探討並分析一個具有耗散聲能吸音牆密閉空間中聲場的阻尼振動模態。利用有限
元素法，我們可由位移場的棱邊離散化將問題轉變為一個求解二次特徵值問題。另一方
面，若考慮壓力節點的離散則會獲得一個有理特徵值問題。透過線性化的技巧，我們可
將這兩個非線性特徵值問題分別改寫成型態為
x
x
的廣義特徵值問題。該問題可以
用 Arnoldi 演算法處理兩種不同型態係數矩陣,
1及
1, 的標準特徵值問題。數值
結果顯示利用 Arnoldi 法求解
1具有較高的精準度。
對於求解二次特徵值問題中絕對值較靠近零之特徵值所對應的特徵對，我們發展了一
個正交投影法－半正交廣義 Arnoldi 法。此外，我們更進一步提出可精化、可重啟動的
半正交廣義 Arnoldi 法。相較於將二次特徵值問題線性化後再利用傳統隱式重啟動
Arnoldi 法求解，數值實驗顯示隱式重啟動半正交廣義 Arnoldi 法
（不論是否有精化過程）
Efficient Arnoldi-Type Algorithms for
Quadratic and Rational Eigenvalue Problems
student：Wei-Qiang Huang Advisors：Dr. Wen-Wei Lin
Department of Applied Mathematics
National Chiao Tung University
ABSTRACT
In this dissertation, we consider two themes related to Arnoldi-type algorithms for
solving nonlinear eigenvalue problems.
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 edge-based finite elements for the displacement field results in a
quadratic eigenvalue problem. On the other hand, the discretization in terms of
pressure nodal finite elements results in a rational eigenvalue problem. We use the
linearization technique to transform these nonlinear eigenvalue problems, respectively,
into generalized eigenvalue problems
x
x
and apply Arnoldi algorithm to two
different types of single matrices
and
. Numerical accuracy shows that the
application of Arnoldi on
is better than that on
.
For computing a few eigenpairs with smallest eigenvalues in absolute value of
quadratic eigenvalue problems, we develop the semiorthogonal generalized Arnoldi
method, an orthogonal projection technique. Furthermore, we propose refinable and
restartable variations of this method to improve the accuracy and efficiency. Numerical
examples demonstrate that the implicitly restarted semiorthogonal generalized Arnoldi
method with or without refinement has superior convergence behaviors than the
implicitly restarted Anoldi method applied to the linearized quadratic eigenvalue
problem.
誌 謝
本篇論文得以完成，首先得感謝我的指導教授，林文偉教授。謝謝老師您近 四年來對我的栽培、鼓勵、鞭策、體諒與照顧。您對學術研究的熱誠、堅持、洞 察力與創造力始終是我敬佩及效仿的對象。儘管我並無太多數值計算的背景知識， 您仍然在旁一步步引領著我進入這豐富的研究領域。當面臨學術或人生所遭遇的 瓶頸與問題，您所分享的研究心得與人生歷練提供我更多不同的思考角度。您對 我生活上的幫助更是讓我無顧之憂的學習。能在您的教導下完成學位論文，是一 件榮耀的事。接著，我要感謝口試委員林松山教授、朱景華教授、王振男教授和 王偉仲教授，於學位論文口試時的提問及建議。特別感謝林松山教授協助我申請 研發替代役讓我得以逐步邁向下個階段中的學習規劃。 謝謝泓勳、育豪、麻將、青松、偉碩、德軒、柏任、函恩、芷瑄、昌翰、劭 芃、安怡、美亨、淑如、明誠學長、Apostol 學長、恭儉學長、文貴學長、明杰 學長、光暉學長、建綸學長、瑜堯、郁傑、建智學長、其棟學長、宏橡（清章） 學長等在我五年博士求學階段中，一同奮鬥過的學長姐、同學夥伴以及學弟妹。 遇見你們並一起在交大奮鬥是一件令人開心與值得回味的事。 感謝王辰樹教授、黃聰明教授、李勇達教授、李宗錂教授、張書銘教授、郭 岳承教授、蔣俊岳教授、黃建智博士、林敏雄教授、吳金典教授、王偉仲教授、 黃印良教授等 Lin Group 的前輩在我求學這段期間帶著我了解、剖析問題，並與 我分享心情和各種大小趣事。此外我得感謝於交大求學這段時間內，每個教導過 我或一同參與研究過程的各位老師，包括莊重教授、賴明治教授、王夏聲教授、 李明佳教授、吳慶堂教授、盧鴻興教授、周所向教授、張企教授、陳泰賓教授等。 同時我也得感謝各位系助理：陳盈吟小姐、張麗君小姐、小張姐、崔妮臻小姐、 宋雅鈴小姐、Billing、李珍韶學姐等在這段時間內，於公於私對我的幫助。 感謝父親黃丁福先生，母親許錦雀女士以及可愛的妹妹春萍、于娗。謝謝父 母養我、育我，並讓我能任性的完成學業。請允許我目前以一張文憑回報你們對 我最無私的關心、鼓勵與支持。 接著，我要謝謝我的女友蕙甄。感謝妳一直以來的陪伴，並體諒我將大多數 的時間投入於研究與學習的無底洞中。謝謝妳陪著我經歷每個無助、徬徨、失望、 開心與平凡的日子。家人與妳是支持我一路來最大的力量也是我最後的避風港。 誠如林文偉教授常告誡我們的一句話：「拿到博士學位沒什麼了不起，學術 研究的一切才剛要開始」。也許迎接我的未來充滿太多的未知數，但我將滿懷感 恩，不怕苦、不怕難，心甘情願地往前走下去。由衷地感謝每個關心我的大家， 謝謝你們。Contents
1 Introduction and Preliminaries 1
1.1 Notations . . . 2
1.2 The Arnoldi Method for Standard Eigenvalue Problems . . . 3
1.3 The Generalized Arnoldi Method for Generalized Eigenvalue Problems 6 1.4 Quadratic Eigenvalue Problems and Linearizations . . . 7
1.5 Rational Eigenvalue Problems and the Trimmed Linearization . . . . 9
2 Arnoldi-Type Algorithms for Nonlinear Eigenvalue Problems Aris-ing in Fluid-Solid Systems 11 2.1 Introduction . . . 12
2.2 The Model Problem . . . 14
2.2.1 Spectral approximation . . . 17
2.3 Linearization of Nonlinear Eigenvalue Problems . . . 20
2.3.1 Linearization of quadratic eigenvalue problems . . . 20
2.3.2 Trimmed linearization for rational eigenvalue problems . . . . 22
2.3.3 Error analysis . . . 25
2.4 Arnoldi Method with Schur-restarting . . . 30
2.4.1 Stopping criteria . . . 33
2.4.2 Computational costs . . . 35
2.5 Numerical Results . . . 38
2.6 Summary . . . 46
3 The Semiorthogonal Generalized Arnoldi (SGA) Method for Quadratic Eigenvalue Problems 49 3.1 Introduction . . . 50
3.2 The SGA Decomposition . . . 54
3.2.1 Existence and uniqueness . . . 55
3.2.2 The SGA algorithm . . . 60
Contents
3.3.1 The SGA method . . . 66
3.3.2 The projection subspace . . . 67
3.4 Refined SGA Method . . . 72
3.5 Implicit Restarting of the SGA Method . . . 74
3.5.1 The IRSGA method and the IRRSGA method . . . 74
3.5.2 The selection of shifts . . . 78
3.6 Numerical Results . . . 79
3.7 Summary . . . 86
4 Conclusions and Future Work 90
List of Figures
2.1 Fluid in a cavity with one absorbing wall and initial mesh . . . 40 2.2 The distribution of the ten desired eigenvalues λ1, . . . , λ10. . . 42 2.3 The #Its of Q1 and Q2 with different shift values. “o” denotes desired
eigenvalues λ1, . . . , λ10. “(i, j)” denotes the #Its for Q1 and Q2, respectively. . . 43 2.4 The relative residuals of computed eigenpairs, obtained by Q1, Q2
for QEP (2.13) and R1, R2 for REP (2.20) with (nℓ, nw) = (768, 576). 43 2.5 The distribution of the ten desired eigenvalues λ1, . . . , λ10 for
Exam-ple 2.2. . . 47 2.6 Relative residuals of computed eigenpairs obtained by R1 and R2 for
REP in Example 2.2 with (nℓ, nw) = (768, 576). . . 47 3.1 Convergence histories for methods ℓ-IRA, r-IRA, IRSGA and
IR-RSGA in Example 3.1. . . 81 3.2 Convergence histories for methods ℓ-IRA, r-IRA, IRSGA and
IR-RSGA in Example 3.2. . . 83 3.3 Convergence histories for methods ℓ-IRA, r-IRA, IRSGA and
List of Tables
2.1 Computational costs of the j-th Arnoldi step of Algorithm 2.1 for
Q-SEP2 and R-SEP2. . . 38
2.2 Dimension information and convergence rates of λ1. . . 41
2.3 #Its for λ1, . . . , λ10 of Q-SEPs and R-SEPs with/without scaling. . . 42
2.4 Iteration numbers and CPU time for Q2 and R2. . . 44
3.1 Iteration numbers and CPU time in Example 3.1. . . 81
3.2 Iteration numbers and CPU time in Example 3.2. . . 82
List of Algorithms
2.1 Arnoldi method with Schur-restarting for solving Ax = λx . . . 32
2.2 Arnoldi method with Schur-restarting for solving QEP in (2.13) . . . 35
2.3 Arnoldi method with Schur-restarting for solving REP in (2.20) . . . 36
3.1 The SGA Algorithm . . . 64
3.2 The SGA method . . . 68
3.3 The RSGA method . . . 74
1
Introduction and Preliminaries
Contents
1.1 Notations . . . 2
1.2 The Arnoldi Method for Standard Eigenvalue Problems 3
1.3 The Generalized Arnoldi Method for Generalized
Eigen-value Problems . . . 6
1.4 Quadratic Eigenvalue Problems and Linearizations . . . 7
1.5 Rational Eigenvalue Problems and the Trimmed
1.1 Notations
The theme explored in this thesis is to develop and exploit efficient Arnoldi-type methods to solve the quadratic and rational eigenvalue problems. This chapter will briefly introduce some basic notions, mathematical notations and conventional methods of the so-called “eigenvalue problems”. We then, in Chapter 2, develop and analyze efficient methods for quadratic and rational eigenvalues arising from computing damped vibration modes of an acoustic fluid confined in a cavity with absorbing walls capable of dissipating acoustic energy. In Chapter 3, we will propose an orthogonal projection method for solving quadratic eigenvalue problems. Finally, conclusions and the future work of this thesis will be discussed in Chapter 4.
1.1
Notations
The following notations are frequently used in this thesis. Other notations will be clearly defined whenever they are used.
• i =√−1.
• We use the symbol ∀ to mean ‘for all’ throughout the thesis.
• R denotes the set of real numbers and C denotes the set of complex numbers. • Re(λ) and Im(λ), respectively, denote the real part and the complex part of
the scalar λ∈ C.
• 0 denotes zero vectors and matrices with appropriate size. • In denotes the n× n identity matrix.
• ej denotes the jth column of the identity matrix In with specified n.
• We use ·⊤ and ·H to denote the transpose and conjugate transpose for vectors or matrices.
1.2 The Arnoldi Method for Standard Eigenvalue Problems
• ⊗ denotes the Kronecker product.
• k · k2, k · kF and k · k∞ respectively denote the 2-norm, Frobenius norm and infinity norm for vectors or matrices.
• We adopt the following MATLAB notations:
v(i : j) denotes the subvector of the vector v that consists of the ith to the jth entries of v;
A(i : j, k : ℓ) denotes the submatrix of the matrix A that consists of the intersection of the rows i to j and the columns k to ℓ;
A(i : j, :) denotes the rows of A from i to j and A(:, k : ℓ) denotes the columns of A from k to ℓ.
1.2
The Arnoldi Method for Standard Eigenvalue
Problems
Given a large sparse matrix A∈ Cn×n, the Arnoldi method [1] is a well known and very prevalent algorithm for solving the so-called standard eigenvalue problem (SEP)
Ax = λx. (1.1)
That is, to find a scalar λ (real or complex) and a nonzero n-vector x satisfying the equations (1.1). In this case, we say that λ is an eigenvalue of A and x is called an eigenvector of A with respect to λ. Moreover, the pair (λ, x) is said to be an eigenpair of A.
Starting with a unit vector v1, the Arnoldi method successively constructs a sequence of unitary vectors v2, v3, . . . , vm which forms a unitary basis of the Krylov
1.2 The Arnoldi Method for Standard Eigenvalue Problems
subspace Km(A, v1)≡ span{v1, Av1, . . . , Am−1v1} with m ≪ n such that hj+1,jvj+1 = Avj − j P i=1 hijvi, j = 1, 2, . . . , m, vH s vt= 0, ∀s 6= t and vsHvs= 1, ∀s, or equivalently, AVm = VmHm+ hm,m+1vm+1e⊤m, h VH m vH m+1 i [Vm vm+1] = Im 0 0 1 , (1.2)
where Vm is an n× m matrix with column vectors v1, v2, . . . , vm, Hm is an m× m upper Hessenberg matrix. After building the factorization (1.2), called the Arnoldi decomposition, we then reduce A into the upper Hessenberg Hm through the unitary transformation VH
mAVm = Hm. The eigenvalues and corresponding eigenvectors of the reduced SEP Hmz = µz can be solved by the classical eigenvalue techniques, such as the QR algorithm (also named the Francis algorithm [18, 19]). Moreover, we see that if (θ, y) is an eigenpair of Hm then (θ, Vmy) is called a Ritz pair of A – an approximate eigenpair of A with the residual norm
k(A − θIn)Vmyk = |hm+1,m||e⊤my|.
For more details on the practical realization and theoretical analysis of the Arnoldi method, we refer to [2, 14, 22, 49, 54, 67, 73].
There are some variations of the Arnoldi method. In practice, a small number of eigenvalues that are nearest to a target σ or located in a prescribed region of the complex plane and the corresponding eigenvectors are often of interest. Under the assumption that σ is not an eigenvalue of the SEP (1.1) but not to far away from the wanted eigenvalues, the shift-and-invert Arnoldi method [48, 54] tends to solve
1.2 The Arnoldi Method for Standard Eigenvalue Problems
the transformed eigenproblem
(A− σIn)−1x= νx, (1.3)
where the scalar value σ is called a shift. It is easy to verify that (1.3) and (1.1) are mathematically equivalent since (ν, x) is an eigenpair of (1.3) if and only if (σ +1ν, x) is an eigenpair of (1.1).
The restarted Arnoldi method aims to overcome the increasing storage as well as the computational cost of the Arnoldi decomposition (1.2) as m is increasing. In [53], Saad coped with these difficulties by developing the explicitly restarted Arnoldi iteration. The idea of this strategy is to compute another mth order Arnoldi decomposition with a “better” initial vector which is a linear combination of some wanted Ritz vectors. The implicitly restarted Arnoldi method [59] and Krylov-Schur algorithm [28, 61, 63], on the other hand, are two remarkable implicitly restarting schemes. These schemes are called implicit due to the fact that the initial vector is sequentially constructed by using the implicitly shifted QR algorithm [18, 19] on the Hessenberg matrix Hm in (1.2). We will review the implementation of the Krylov-Schur restarting in Section 2.4.
Another possible problem is that even though some desirable eigenvalues com-puted by the Arnoldi method already attempt to converge, the corresponding ap-proximate eigenvectors may converge very slowly and even fail to converge. The refined Arnoldi method [32] proposed by Jia gave an alternative approach to rem-edy this problem by computing refined approximate eigenvectors. See also [33]. We will mimic this idea and design a refinement strategy for our Arnoldi-type method in Section 3.4. Other variations of the Arnoldi method include the block-Arnoldi method [55], the inexact Arnoldi method [56], the residual Arnoldi method [37, 38], and so on.
1.3 The Generalized Arnoldi Method for Generalized Eigenvalue Problems
1.3
The Generalized Arnoldi Method for
General-ized Eigenvalue Problems
The generalized eigenvalue problem (GEP) for the matrix pencil A− λB of two square matrices A and B with size n is to determine scalars λ and n-vectors x6= 0 such that
Ax = λBx. (1.4)
If B is nonsingular, the GEP (1.4) can be transformed into SEPs
(B−1A)x = λx (1.5)
or
(AB−1)y = λy, y= Bx (1.6)
and subsequently solved by the standard Arnoldi method. Alternatively, the QZ algorithm [45], an analog of the QR algorithm for the GEP, is the method of choice for dealing with the GEP (1.4) with small dense coefficient matrices.
The truncated QZ method proposed by Sorensen [60] is one of the approaches for solving large-scale GEPs. For m ≪ n, this method constructs a generalization of the standard Arnoldi decomposition (1.2),
AZm = YmHm+ hm+1,mym+1e⊤m, BZm = YmRm, ZH mZm = Im, YmHYm = Im, YmHym+1 = 0, (1.7)
which is called the generalized Arnoldi reduction in [60], and deals with the small-sized GEP Hmv= µRmv of the m× m upper Hessenberg-triangular pair (Hm, Rm) to approximate eigenpairs of the original large-scale GEP (1.4).
1.4 Quadratic Eigenvalue Problems and Linearizations
1.4
Quadratic Eigenvalue Problems and
Lineariza-tions
In this section, we consider the quadratic eigenvalue problem (QEP) of the form
Q(λ)x≡ (λ2M + λD + K)x = 0, (1.8)
where M, D and K are n× n large and sparse matrices. The QEP is a special case of the polynomial eigenvalue problem (PEP)
P (λ)x ≡ d X i=0 λiPi ! x= 0, (1.9)
where Pi are constant matrices of size n for 1 ≤ i ≤ d. P (λ) ≡ d P i=0
λiP
i is called a matrix polynomial (in λ) of degree d. Obviously, for d = 0, 1 and 2, the PEP (1.9) is, respectively, indeed the case of SEP (1.1), GEP (1.4) and QEP (1.8).
The “linearization” is a typical and most widely used technique to solve the QEP in which the problem is reformulated into a linear one which doubles the order of the system. By selecting suitable matrices A, B ∈ C2n×2n and the vector ϕ ∈ C2n, we can convert (1.8) into the GEP
(A − λB)ϕ = 0 (1.10)
satisfying the relation
E(λ)(A − λB)F(λ) = Q(λ) 0 0 In ,
1.4 Quadratic Eigenvalue Problems and Linearizations
determinants. In this case,
det(A − λB) = det(λ2M + λD + K)
indicates the eigenvalues of the original QEP (1.8) coincide with the eigenvalues of the enlarged GEP (1.10). As a result, the linearization technique of QEPs makes classical methods for GEPs as well as SEPs can be used.
There are many choices of (A, B)’s, but probably the most popular ones in prac-tice are the so-called companion forms [21]: the first companion form
A = −D −K In 0 and B = M 0 0 In
as well as the second companion form
A = −D In −K 0 and B = M 0 0 In . (1.11)
There are some drawbacks, however, of the linearization technique to solve QEPs. For instance, the doubling size of the problem increases the computational cost and the original structures of the coefficient matrices (M, D, K) such as symmetry and positive definiteness may be lost. To circumvent these drawbacks, one may expect to solve the QEP (1.8) directly. The QEP is projected onto a properly chosen low-dimensional subspace in order to lower the matrix sizes of the coefficient matrices in (1.8). The reduced QEP can then be solved by a standard approach for dense matrices. Methods of this type include the residual iteration method [27, 43, 47], the Jacobi-Davidson method [57, 58], a Krylov-type subspace method [40], the nonlinear Arnoldi method [68], the second-order Arnoldi method [3, 41, 72] and an iterated
1.5 Rational Eigenvalue Problems and the Trimmed Linearization
shift-and-invert Arnoldi method [75]. While these methods use a similar projection process, the main difference between them is the selection of projection subspaces. Convergence analysis of projection methods to approximate eigenpairs of the QEP (1.8) has recently been invented in [29].
In Chapter 2, we consider a QEP arising from a finite element model and convert it into SEPs in (1.5) and (1.6) through an equivalent second companion form (1.11). We will report theoretical and numerical comparisons of these two SEPs (1.5) there. In Chapter 3, we combine the generalized Arnoldi reduction (1.7) and the second companion form linearization (1.11) of the QEP (1.8) to develop a projection method to solve the QEP (1.8) directly.
1.5
Rational Eigenvalue Problems and the Trimmed
Linearization
The rational eigenvalue problem (REP) concerns the problem of finding (λ, x) with x6= 0 satsifying the equation
R(λ)x≡ P (λ)− r X j=1 sj(λ) tj(λ) Cj ! x= 0, (1.12)
where P (λ) is an n×n matrix polynomial in λ, sj(λ) and tj(λ) are scalar polynomials in λ, and Cj are n× n constant matrices. The simulation of the three-dimensional pyramid quantum dot heterostructure [31] produces a REP. For more examples related to the REPs, see [9, 44, 64].
To solve the REP (1.12), one may immediately multiply (1.12) by the scalar polynomial
r Q j=1
tj(λ) to convert it into a PEP. Subsequently, the PEP can be lin-earized to a GEP. The nonlinear eigensolver, on the other hand, is another way to
1.5 Rational Eigenvalue Problems and the Trimmed Linearization
solve this problem. The nonlinear Jacobi-Davidson method [70] and the nonlinear Arnoldi method [68] fall into this category. Yet, these approaches also have the prob-lem that restricts advantages of the underlying matrix structures and properties of REPs.
Trimmed linearization [64] is a recent linearization-based approach to solve the REP (1.12), especially when matrices Cj in (1.12) have the low-rank property. This method utilizes and preserves the structure and property of the REP (1.12) as much as possible to transform it into a GEP (and hence a SEP), and it only slightly increases the size of the GEP, compared to the size of the original REP (1.12).
The REP discussed in this thesis is a quadratic matrix polynomial (P (λ) in (1.9) with d = 2) together with low-rank rational terms (see Eq. (2.20)). We will use the trimmed linearizing skill and the shift-and-invert Arnoldi method with Krylov-Schur restarting to detect desired eigenpairs.
2
Arnoldi-Type Algorithms for
Nonlinear Eigenvalue Problems
Arising in Fluid-Solid Systems
Contents
2.1 Introduction . . . 12
2.2 The Model Problem . . . 14
2.3 Linearization of Nonlinear Eigenvalue Problems . . . 20
2.4 Arnoldi Method with Schur-restarting . . . 30
2.5 Numerical Results . . . 38
2.1 Introduction
2.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 methods to compute the modes. This phase involves correct discretization of the mathemati-cal formulation and computation of large smathemati-cale 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 mathemat-ical 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 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 [4, 12, 20, 24, 71], among others.
One of the discretizations we will be using in this chapter is the edge-based or Raviart-Thomas finite elements for the displacement field, following [5, 7]. The main concerns in [5, 51] are pure mathematical issues of proving that their numerical approximation is free of spurious modes and has second order convergence rate.
2.1 Introduction
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 chapter 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 [5, 7], but we further develop efficient methods for solving the associated QEP. However, we show in Section 2.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. While there is an extensive literature on QEPs problems [66], REPs are much less studied [64, 68, 69]. Although on the surface the REP (Eq. (2.20)) 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 2.3. The organization of this chapter is as follows. We describe the underlying model fluid-solid problem of this chapter in Section 2.2, where the edge-based displace-ment 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 Sec-tion 2.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
2.2 The Model Problem
is as usual turned into a generalized eigenvalue problem (GEP), from which two types of standard eigenvalue problems (SEP) (2.19.1) and (2.19.2) are derived. The REP is trimmed-linearized into two types of three by three block SEPs (2.31.1) and (2.31.2). The important issue of residual error bound analysis is addressed here. We then apply Arnoldi method with Schur-restarting described in Section 2.4 to 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 Section 2.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. Summaries are included in Section 2.6.
2.2
The Model Problem
Let us consider a simple model of a rigid container filled with an inviscid com-pressible 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 Ω⊂ Rd(d = 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 ΓR is the remaining part of Γ. An example of the setup is in Figure 2.1(i) on Section 2.5, where the top boundary is absorbing and the remaining boundary is rigid.
2.2 The Model Problem
displacement field U, which satisfy ([8, 35])
ρ∂ 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, the acoustic speed, and n, the unit outer normal vector along ∂Ω. 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λtp(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.
2.2 The Model Problem
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 L2 vector functions v on Ω with L2 integrable divergence.
Testing (2.5) by ¯v ∈ V and integrating by parts, we obtain a variational for-mulation of problem (2.5)–(2.8) involving only the displacement variable: Find λ∈ C and u ∈ V, u 6= 0, such that
λ2 Z Ω ρu· ¯v + λ Z ΓA βu· n¯v · n + Z ΓA αu· n¯v · n + Z Ω ρc2 divu div¯v= 0 ∀ v ∈ V. (2.9)
This is a quadratic eigenvalue problem. Note that λ = 0 is an eigenvalue and the dimension of its eigenspace
N := {u ∈ V : div u = 0 in Ω and u · n = 0 on ∂Ω}
is infinity. All nonzero eigenvalues have finite multiplicity (the dimension of the eigenspace is finite) [6]. It is shown in [6] that all the other solutions of (2.9), the decay rate is strictly negative. That is, if an eigenpair 06= λ ∈ 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(Ω) :={p ∈ L2(Ω) : ∇p ∈ L2(Ω)} such that
λ2 c2 Z Ω p¯q + λ 2 α + λβ Z ΓA ρp¯q + Z Ω∇p · ∇¯q = 0 ∀ q ∈ H 1(Ω). (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
2.2 The Model Problem
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.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. [6] successfully demonstrated that the spurious modes can be avoided by using the lowest order Raviart-Thomas elements in Rd, d = 2, 3 (see, for instance, [10, 50]). 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 ∈ P0d⊕ P0x ∀ T ∈ Th and vh· n = 0 on ΓR} ⊂ V,
where d = 2 andPkdenotes the set of polynomials of degree at most k. Thus locally vh takes the form (a + sx, b + sy)⊤. The discrete problem associated with (2.9) is : Find λ∈ C and uh ∈ Vh, uh 6= 0, such that
λ2 Z Ω ρuh· ¯vh+ λ Z ΓA βuh· n ¯vh· n + Z ΓA αuh· n ¯vh· n + Z Ω ρc2divuh div¯vh= 0, ∀ vh∈ Vh. (2.11)
Theorem 2.1. The dimension of the zero eigenspaceE0associated with (2.11) equals the number of interior nodes in the triangulation.
Proof. Setting vh = uh and λ = 0 in (2.11), we see that
2.2 The Model Problem
Since uh = (a + sx, b + sy)⊤ on T ∈ Th, the divergence free condition implies that uh is a constant vector (a, b)⊤ 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)⊤ 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 = ∇ψT · t = ∂ψT
∂t .
So if an edge e is common to T1 and 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) ⊤ and define ∇⊥S
h :={∇⊥Ψh : Ψh is continuous piecewise linear and vanishes on the boundary}.
Thus we have just shown the zero eigenspaceE0 is contained∇⊥Sh and the opposite inclusion is also easily checked. Hence
E0 =∇⊥Sh.
We now find the dimension of ∇⊥S
h. 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,
2.2 The Model Problem
suppose N P j=1
cj∇⊥Ψj = 0. Then this implies N P j=1
cjΨj = c for some constant c. Hence cj = c by the condition Ψj(Nk) = δkj. Consequently, c(PjΨj − 1) = 0. But we know
N P 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 dimE0 = dim∇⊥Sh equals the number of interior nodes in the mesh.
Define the conforming P1 finite element space
Hh :={ph ∈ H1(Ω) : ph|T ∈ P1 ∀ T ∈ Th}.
This is the subspace of H1(Ω) consisted of continuous piecewise linears. The alter-native discrete problem in terms of the approximate pressure field is: Find λ ∈ C and ph ∈ Hh such that
λ2 c2 Z Ω phq¯h+ λ2 α + λβ Z ΓA ρphq¯h+ Z Ω∇p h· ∇¯qh = 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 (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 addition we show in the remaining sections that its associated eigenvalue problem can be efficiently solved. A minor remark is in order here.
Remark 2.2. 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 linear equations again so as to maintain the advantage of the pressure formulation. It should be given by a simple
2.3 Linearization of Nonlinear Eigenvalue Problems
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. Con-sequently, ∇ · 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 [13] 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.
2.3
Linearization of Nonlinear Eigenvalue Problems
In this section we start to address the computational issues related to the dis-placement approximation (2.11) and the pressure approximation (2.12).
2.3.1
Linearization of quadratic eigenvalue problems
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, Reiφj · nidς = δij. For uh ∈ Vh, we write uh = n1 P j=1 ujφj and denote u = [u⊤1,· · · , u⊤n1] ⊤. 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:
Q(λ)u≡ (λ2Mu+ (α + λβ)Au+ Ku)u = 0, (2.13)
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, (2.14)
2.3 Linearization of Nonlinear Eigenvalue Problems
for i, j = 1, . . . , n1. For this problem, we are only interested in eigenvalues that are located in the interior of the spectrum. Suppose that the eigenvalues near σ are of interest. Accordingly, the QEP (2.13) is shifted into
µ2Mfu+ µ eDu+ eKu u= 0 (2.15) with µ = λ− σ and f Mu = Mu, e Du = 2σMu+ βAu, e Ku = σ2Mu+ (α + σβ)Au+ Ku. (2.16)
On the one hand, one can numerically solve (2.15) without transforming it further. Among such direct methods we mention the second-order Arnoldi (SOAR) method [3] and the Jacobi-Davidson algorithm applied to polynomial eigenvalue problems [57]. On the other hand, it is more common to transform or linearize (2.15) into a SEP [66]. In this chapter, we let
Au = 0 −fMu In1 − eDu , Bu = In1 0 0 Keu (2.17)
and linearize (2.15) into the GEP
Auϕ= 1 µBuϕ with ϕ≡ −µfMuu u ≡ v u . (2.18)
The matrix eKu in (2.17) is nonsingular owing to the fact that the shift value σ is not an eigenvalue of (2.13). Furthermore, the GEP (2.18) can then be transformed into two types of SEPs of the forms (B−1
2.3 Linearization of Nonlinear Eigenvalue Problems
respectively, where ψ =Buϕ. Therefore, from (2.17) and (2.18) we have
(Q-SEP1) B−1 u Au v u = 0 − fMu e K−1 u − eKu−1Deu v u = 1µ v u (2.19.1) and (Q-SEP2) AuBu−1 v w = 0 − fMuKeu−1 In1 − eDuKeu−1 v w =µ1 v w , w = eKuu.(2.19.2)
Note that the SEPs of (2.19.1) and (2.19.2) derived by the QEP in (2.15), 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 Section 2.4.
2.3.2
Trimmed linearization for rational eigenvalue problems
Let {ψj}nj=12 be a nodal basis of Hh. For ph ∈ Hh, we write ph =
n2
P j=1
pjψj and denote p = [p1,· · · , pn2]
⊤. Then, the discrete problem (2.12) can be written as the following REP: R(λ)p≡ λ2 c2Mp+ Kp+ λ2 λβ + αAp p= 0, (2.20)
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 (2.21) for i, j = 1, . . . , n2.
To solve REP (2.20), one approach is to multiply equation (2.20) by the scalar λβ+α and expand it into a cubic polynomial eigenvalue problem, and then solve it by
2.3 Linearization of Nonlinear Eigenvalue Problems
Jacobi-Davidson method [30]. An alternative approach is to treat (2.20) as nonlinear eigenvalue problem and solve it by a nonlinear eigensolver, such as Newton’s method, nonlinear Arnoldi method, or nonlinear Jacobi-Davidson method [52, 68, 69]. Re-cently, a trimmed linearization is proposed in [64] which linearizes (2.20) into a GEP so that the standard Arnoldi method can be applied. We introduce the trimmed linearization below.
Given a shift value σ. With µ = λ− σ, the rational λ-matrix R(λ) in (2.20) 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σ c2 Mp + σ2 c2Mp+ Kp +µ 2+ 2µσ + σ2 µβ + σβ + α Ap. (2.22)
By applying the long division, the rational term in (2.22) can be simplified into the following µ2+ 2µσ + σ2 µβ + σβ + α = µ 2 α2 (σβ + α)3 + µ σ2β + 2σα (σβ + α)2 + σ 2 σβ + α − µ 2 (σβ + α)3 α2 + (σβ + α)4 α2βµ −1 .
This implies that
R(λ) = µ2 1 c2Mp + α2 (σβ + α)3Ap + µ 2σ c2 Mp+ σ2β + 2σα (σβ + α)2 Ap + σ2 c2Mp + Kp+ σ2 σβ + αAp − µ2 (σβ + α)3 α2 + (σβ + α)4 α2βµ −1 Ap = µ2Mfp+ µ eDp+ eKp− µ2 ϑ− ̺µ−1 −1 LpR⊤p, (2.23)
2.3 Linearization of Nonlinear Eigenvalue Problems where f Mp = 1 c2Mp+ α2 (σβ + α)3Ap, (2.24) e Dp = 2σ c2Mp+ σ2β + 2σα (σβ + α)2Ap, (2.25) e Kp = σ2 c2Mp+ Kp+ σ2 σβ + αAp, (2.26) ϑ = (σβ + α) 3 α2 , ̺ = − (σβ + α)4 α2β , (2.27)
and LpR⊤p = Ap is the full-rank decomposition of Ap with Lp, Rp ∈ Rn2×ℓ, ℓ ≪ n2. Introducing an auxiliary vector
q= µ
ϑµ− ̺R ⊤
pp, (2.28)
the REP in (2.20) can be reformulated as
µ2Mfp+ µ eDp+ eKp
p− µ2Lpq= 0. (2.29)
Using (2.28) and (2.29), we get the GEP
Apϕ ≡ 0 −fMp Lp In2 − eDp 0 0 −Rp⊤ ϑIℓ ϕ= 1 µ In2 0 0 0 Kep 0 0 0 ̺Iℓ ϕ≡ 1 µBpϕ, (2.30) where ϕ = [((µ−1Ke
p + eDp)p)⊤, p⊤, q⊤]⊤. As before, the matrix eKp in (2.26) is nonsingular due to the fact that the shift value σ is not an eigenvalue of (2.20). As in (2.19.1) and (2.19.2), the GEP (2.30) can then be, respectively, transformed into the following two types of the SEPs of the forms (B−1
2.3 Linearization of Nonlinear Eigenvalue Problems µ−1ψ where ψ =B pϕ. Consequently, we have (R-SEP1) B−1 p Apϕ= 0 − fMp Lp e Kp−1 − eKp−1Dep 0 0 −̺−1R⊤p ̺−1ϑIℓ ϕ= 1 µϕ, (2.31.1) and (R-SEP2) ApBp−1ψ= 0 − fMpKep−1 ̺−1Lp In2 − eDpKe −1 p 0 0 −R⊤ pKep−1 ̺−1ϑIℓ ψ = 1 µψ, ψ=Bpϕ. (2.31.2)
Note that the SEPs of (2.31.1) and (2.31.2) derived by the REP in (2.29) are called R-SEP1 and R-SEP2, respectively.
2.3.3
Error analysis
In this subsection, we will discuss residuals of QEP (2.13) and REP (2.20) by using linearizations (2.19) and (2.31), respectively.
We first derive residual bounds of approximate eigenpairs for QEP (2.13) by by using linearizations Q-SEP1 and Q-SEP2, respectively. Let (µ−11 ,hv1
u1
i
) be an approximate eigenpair of (2.19.1) and hf11
f12
i
be the associated residual vector. That is, f11 f12 = 0 −fMu e K−1 u − eKu−1Deu v1 u1 −µ1 1 v1 u1 = 1 µ1 −v1 − µ1 f Muu1 e K−1 u (µ1v1 − µ1Deuu1− eKuu1) .
2.3 Linearization of Nonlinear Eigenvalue Problems
It follows that
µ12Mfuu1+ µ1Deuu1+ eKuu1 = µ1(−v1− µ1f11) + µ1v1− µ1Keuf12 = −µ2
1f11− µ1Keuf12.
Let λ1 = µ1+ σ. From (2.13) we have kQ(λ1)u1k ku1k = kµ 2 1Mfuu1+ µ1Deuu1+ eKuu1k ku1k ≤ |µ1|2kf11k + |µ1|k eKukkf12k ku1k . (2.32)
On the other hand, let (µ−12 , h
v2
w2
i
) be an approximate eigenpair of (2.19.2) andhf21
f22
i be the associated residual vector. That is,
f21 f22 = 0 −fMuKeu−1 I − eDuKeu−1 v2 w2 − µ1 2 v2 w2 = −fMuKe −1 u w2− µ12v2 v2 − eDuKeu−1w2− µ12w2 . It follows that µ22MfuKeu−1w2+ µ2DeuKeu−1w2+ w2 = µ2(−v2 − µ2f21) + µ2v2− µ2f22 = −µ22f21− µ2f22.
Letting u2 = eKu−1w2 and λ2 = µ2+ σ. From (2.13) we have, kQ(λ2)u2k ku2k = kµ 2 2Mfuu2+ µ2Deuu2+ eKuu2k ku2k ≤ |µ2|2kf21k + |µ2|kf22k ku2k . (2.33)
Now, we derive residual bounds of approximate eigenpairs for REP (2.20) by us-ing linearizations R-SEP1 and R-SEP2, respectively. Let (µ−11 , [s⊤
2.3 Linearization of Nonlinear Eigenvalue Problems
an approximate eigenpair of (2.31.1) and [g⊤
11, g⊤12, g⊤13]⊤ be the associated residual vector. That is,
g11 g12 g13 = 0 −fMp Lp e K−1 p − eKp−1Dep 0 0 −̺−1R⊤ p ̺−1ϑIℓ s1 p1 q1 − 1 µ1 s1 p1 q1 .
This implies that
s1 = −µ1Mfpp1+ µ1Lpq1− µ1g11, (2.34) g12 = Kep−1s1− eKp−1Depp1− 1 µ1 p1, (2.35) q1 = µ1̺−1ϑ− 1 −1 µ1 g13+ ̺−1Rp⊤p1 . (2.36)
Substituting (2.36) into (2.34), s1 can be represented by
s1 =−µ1Mfpp1+ µ21 µ1̺−1ϑ− 1 −1
Lpg13+ ̺−1LpR⊤pp1
− µ1g11. (2.37)
Substituting (2.37) into (2.35) and taking λ1 = µ1+ σ. From (2.23) and (2.27),
R(λ1)p1 = µ21Mfpp1+ µ1Depp1+ eKpp1− µ21 (σβ + α)3 α2 + (σβ + α)4 α2βµ 1 −1 LpR⊤pp1 = −µ2 1g11− µ1Kepg12− µ21 β σβ + α + 1 µ1 −1 Lpg13
which implies that kR(λ1)p1k kp1k ≤ 1 kp1k ( |µ1|2kg11k + |µ1|k eKpk kg12k + µ 2 1 β σβ + α + 1 µ1 −1 kLpk kg13k ) . (2.38)
2.3 Linearization of Nonlinear Eigenvalue Problems
On the other hand, let (µ−12 , [s⊤
2, t⊤2, q⊤2]⊤) be an approximate eigenpair of (2.31.2) and [g⊤
21, g⊤22, g23⊤]⊤ be the associated residual vector. That is, g21 g22 g23 = 0 −fMpKep−1 ̺−1Lp In2 − eDpKe −1 p 0 0 −R⊤ pKep−1 ̺−1ϑIℓ s2 t2 q2 − 1 µ2 s2 t2 q2 .
This implies that
g21 = −fMpKep−1t2+ ̺−1Lpq2− 1 µ2 s2, (2.39) s2 = DepKep−1t2+ 1 µ2 t2+ g22, (2.40) q2 = ̺−1ϑ− 1 µ2 −1 R⊤pKep−1t2+ g23 . (2.41)
Substituting (2.40) and (2.41) into (2.39), we have
µ22MfpKep−1t2+ µ2DepKep−1t2+ t2− µ22 ϑ− ̺µ−12 −1 LpRp⊤Kep−1t2 = −µ22g21− µ2g22+ µ22 ϑ− ̺µ−12 −1 Lpg23.
Letting p2 = eKp−1t2 and setting λ2 = µ2+ σ. From (2.23) we get
R(λ2)p2 = µ22Mfpp2+ µ2Depp2+ eKpp2 −µ22 (σβ + α)3 α2 + (σβ + α)4 α2βµ 2 −1 LpR⊤pp2 = −µ22g21− µ2g22+ µ22 α2β (σβ + α)4 β σβ + α + 1 µ2 −1 Lpg23.
2.3 Linearization of Nonlinear Eigenvalue Problems Hence, kR(λ2)p2k kp2k ≤ 1 kp2k ( |µ2|2kg21k + |µ2|kg22k + µ 2 2 α2β (σβ + α)4 β σβ + α + 1 µ2 −1 kLpkkg23k ) . (2.42)
Remark 2.3. In order to check the tightness of upper bounds in (2.32) and (2.33), as well as, (2.38) and (2.42) for residuals, respectively, we refer to the coefficient matrices generated in Example 2.1 of Section 2.5. For (2.9) we adopt the data as in [6] by setting ρ = 1kg/m3, c = 340 m/s, α = 5× 104 N/m3, and β = 200 Ns/m3. In addition, we choose σ =−25 + 600πi as the shift value. Then
(i) from (2.14), the element mass and stiffness matrices are
h2 6 ρ 2 −1 0 −1 2 0 0 0 2 and ρc 2 2 2 2√2 2 2 2√2 2√2 2√2 4 ,
respectively. Hence, by (2.16) the infinity norm of eKu can be estimated by k eKuk∞ ≈ kKuk∞ =O(ρc2) = O(105). From (2.32) and (2.33), we conclude that the upper bound for the residual of the approximate eigenpair (µ1+ σ, u1) of (2.13) by solving Q-SEP1 is larger than that of the approximate eigenpair (µ2+ σ, u2) of (2.13) by solving Q-SEP2.
(ii) From (2.21), the element mass and stiffness matrices are
h2 24 2 1 1 1 2 1 1 1 2 and 1 −1/2 −1/2 −1/2 1/2 0 −1/2 0 1/2 ,
2.4 Arnoldi Method with Schur-restarting
respectively. Hence, by (2.26) we have that k eKpk∞ ≈ kKpk∞ =O(1). If the eigenvalue λ is one of the desired eigenvalues in Figure 2.2, then with µ = λ−σ we have 4× 107 < µ 2 β σβ + α + 1 µ −1 < 3.1× 10 10 and 0.001 < µ 2 α2β (σβ + α)4 β σβ + α + 1 µ −1 < 0.8.
Clearly, from (2.38) and (2.42) we conclude that the upper bound for the resid-ual of the approximate eigenpair (µ1+σ, p1) of REP (2.20) by solving R-SEP1 is larger than that of (µ2+ σ, p2) of (2.20) by solving R-SEP2.
2.4
Arnoldi Method with Schur-restarting
The Arnoldi method is the most popular method for solving large sparse SEPs: Ax = λx. In Arnoldi process, an orthonormal matrix Vm+1 is generated to satisfy
AVm = VmHm+ hm+1,mvm+1e⊤m, (2.43)
where Hm ∈ Cm×m is an upper Hessenberg matrix. If the dimension of the Krylov subspace span{Vm} is larger than a certain value, then the process of Arnoldi de-composition will be restarted.
For the restarting process, we can use an implicit restart scheme [46, 59]. The package ARPACK [39] includes a very successful implementation of the implicitly restarted Arnoldi algorithm. It has been used by numerous 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 [28, 61, 63] that relaxes the need to
2.4 Arnoldi Method with Schur-restarting
preserve the structure of the Arnoldi decomposition and therefore ease the compli-cations of the purging and deflating.
We state the Schur-restarting scheme as follows. Let
Hm = [Uk Uℓ] Tk Tf 0 Tℓ U H k UH ℓ (2.44)
be a Schur decomposition of Hm where Tk and Tℓ are upper triangular, and the eigenvalues of Tk are of interest. Substituting (2.44) into (2.43), we see that
A(Vm[Uk Uℓ]) = (Vm[Uk Uℓ]) Tk Tf 0 Tℓ + hm+1,mvm+1(e⊤m[Uk Uℓ]),
which implies that
A ˜Vk= ˜VkTk+ ˜vk+1tHk, (2.45)
where ˜Vk ≡ VmUk, ˜vk+1 = vm+1 and tHk ≡ hm+1,me⊤mUk.
Let Q1 be a Householder matrix with tHkQ1 = τ e⊤k. Then (2.45) can be rewritten as
A( ˜VkQ1) = ( ˜VkQ1)(Q1HTkQ1) + τ ˜vk+1e⊤k. (2.46)
The matrix QH
1 TkQ1 can be reduced to a new Hessenberg matrix Hk+ by using Householder matrices Qi for i = 2, . . . , k− 1 with
QHk−1· · · QH
2 (Q
H
1 TkQ1)Q2· · · Qk−1 = Hk+ e⊤kQ2· · · Qk−1 = e⊤k.
2.4 Arnoldi Method with Schur-restarting
Multiplying (2.46) by Qi, i = 2, . . . , k− 1, a new Arnoldi decomposition of order k
AVk+ = Vk+Hk++ τ v+k+1e⊤k
is obtained where V+
k := ˜VkQ1· · · Qk−1, vk+1+ = ˜vk+1 = vm+1 and the Arnoldi process can be applied to generate it to order m in (2.43). One repeats the above process until the desired eigenvalues are convergent. The process is summarized in Algorithm 2.1.
Algorithm 2.1 Arnoldi method with Schur-restarting for solving Ax = λx
Input: A: coefficient matrix, tolA: tolerance for convergence, rmax: maximum num-ber of Schur-restartings.
Output: The desired k eigenpairs.
1: Build an initial Arnoldi decomposition of order k as in (2.43) and set r = 0.
2: restart
3: Extend Arnoldi decomposition of order k to order m = k + ℓ and set r = r + 1. 4: Compute all Ritz pairs (µ−1i , zi) with Hkzi = µ−1i zi, i = 1, . . . , m and sorting
Ritz values so that {(µ1, z1), . . . , (µk, zk)} are wanted.
5: for i = 1, . . . , k do
6: Check convergence by |hm+1,m||e⊤mzi| < tolA.
7: end for
8: if ( Not all m desired eigenvalues are convergent and r < rmax ) then
9: Compute the Schur decomposition of Hm as in (2.44), where the eigenvalues of Tk are of interest.
10: Set Vk:= VmUk, vk+1 := vm+1 and tHk := hm+1,me⊤mUk.
11: Compute Householder transformation Q1 such that tHkQ1 = τ e⊤k.
12: Reduce QH
1 TkQ1 to a new Hessenberg matrix Hk by using Householder transformations Qi for i = 2, . . . , k− 1.
13: Set Vk := VkQ1· · · Qk−1 and hk+1,k = τ to get the new Arnoldi decomposi-tion with order k:
AVk = VkHk+ hk+1,kvk+1e⊤k. (2.47)
14: end if
15: until ( desired k eigenpairs are convergent or r ≥ rmax )
Now, we will apply the Algorithm 2.1 to solve QEP (2.13) and REP (2.20), respectively, by setting A to be the coefficient matrices in (2.19) and (2.31), respec-tively.
2.4 Arnoldi Method with Schur-restarting
2.4.1
Stopping criteria
Let (µ−1, z) be a Ritz pair and satisfy H
mz= µ−1z. From (2.43) and Q-SEP1 in (2.19.1) we have 0 −fMu e K−1 u − eKu−1Deu Vm1 Vm2 z = µ1 Vm1 Vm2 z + hm+1,m vm+1,1 vm+1,2 e⊤mz, (2.48) where Vm = h Vm1 Vm2 i and vm+1 = h vm+1,1 vm+1,2 i
are partitioned with compatible sizes. Using the first equation of (2.48), we can eliminate Vm1z in the second equation and get
kQ(λ)u1k ku1k = k(µ 2Mf u+ µ eDu+ eKu)u1k ku1k = |µ| |hm+1,m| e⊤ mz ζ1 ku1k ≡ q 1(µ), (2.49)
where u1 = Vm2z, λ = µ + σ and ζ1 = kµvm+1,1+ eKuvm+1,2k. Without ambiguity by using the same notations as above in Algorithm 2.1, from (2.43) and Q-SEP2 in (2.19.2) we also have 0 −fMuKeu−1 In1 − eDuKe −1 u Vm1 Vm2 z = µ1 Vm1 Vm2 z + hm+1,m vm+1,1 vm+1,2 e⊤mz and kQ(λ)u2k ku2k = k(µ 2Mf u+ µ eDu+ eKu)u2k ku2k = |µ| |hm+1,m| e⊤ mz ζ2 ku2k ≡ q 2(µ), (2.50)
where u2 = eKu−1Vm2z, λ = µ + σ and ζ2 =kµvm+1,1+ vm+1,2k. Therefore, q1(µ) in (2.49) and q2(µ) in (2.50), respectively, can be used as stopping criteria for residuals while Algorithm 2.1 is applied to solved QEPs (2.13).
Similarly, we can apply Algorithm 2.1 to solve REPs (2.20). As above, we let (µ−1, z) be a Ritz pair and satisfy H
R-2.4 Arnoldi Method with Schur-restarting SEP2 in (2.31) we have 0 − fMp Lp e K−1 p − eKp−1Dep 0 0 −̺−1R⊤p ̺−1ϑIℓ Vm1 Vm2 Vm3 z= 1 µ Vm1 Vm2 Vm3 z+ hm+1,m vm+1,1 vm+1,2 vm+1,3 e⊤mz (2.51) and 0 − fMpKep−1 ̺−1Lp In2 − eDpKep−1 0 0 −R⊤ pKep−1 ̺−1ϑIℓ Vm1 Vm2 Vm3 z= 1 µ Vm1 Vm2 Vm3 z+ hm+1,m vm+1,1 vm+1,2 vm+1,3 e⊤mz, (2.52)
where Vm = [Vm1⊤, Vm2⊤, Vm3⊤]⊤ and vm+1 = [v⊤m+1,1, vm+1,2⊤ , v⊤m+1,3]⊤ are partitioned with compatible sizes. Using the first and the third equations of (2.51) and (2.52), we can eliminate V1z and V3z in the second equation of (2.51) and (2.52), respectively, and get kR(λ)p1k kp1k = k[µ 2Mf p + µ eDp + eKp− µ2(ϑ− ̺µ−1)−1Ap]p1k kp1k = |µ| |hm+1,m| e⊤ mz ξ1 kp1k ≡ r1 (µ), (2.53) where p1 = Vm2z, λ = µ + σ and ξ1=kµvm+1,1+ eKpvm+1,2− ̺µ 2 ϑµ−̺Lpvm+1,3k, and kR(λ)p2k kp2k = k h µ2Mf p + µ eDp + eKp− µ2(ϑ− ̺µ−1)−1Ap i p2k kp2k = |µ| |hm+1,m| e⊤ mz ξ2 kp2k ≡ r 2(µ), (2.54) where p2 = eKp−1Vm2z, λ = µ + σ and ξ2 = kµvm+1,1 + vm+1,2 − µ2 ϑµ−̺Lpvm+1,3k.
Therefore, r1(µ) in (2.53) and r2(µ) in (2.54) can be used as stopping criteria for residuals while Algorithm 2.1 is applied to solve REPs (2.20).
2.4 Arnoldi Method with Schur-restarting
in Algorithm 2.2 and Algorithm 2.3, respectively.
Algorithm 2.2 Arnoldi method with Schur-restarting for solving QEP in (2.13) Input: Coefficient matrices Mu, Du and Ku, parameters c, α and β, σ: shift value,
tolQ: tolerance for convergence, rmax: maximum number of Schur-restartings. Output: The desired eigenpairs (λi, ui) for i = 1, . . . , k.
1: Construct matrices fMu, eDu and eKu defined in (2.16) and set r = 0.
2: Compute initial Arnoldi decomposition in Line 1 of Algorithm 2.1 with A in
Q-SEP1 or Q-SEP2.
3: restart
4: Do the steps in Lines 3 and 4 of Algorithm 2.1. 5: for i = 1, . . . , k do
6: Compute
ϕ(µi) = (|σ + µ−1i |2kMuk + |α + (σ + µ−1i )β|kAuk + kKuk).
7: Check convergence of QEP by
qℓ(µi) ϕ(µi)
< tolQ with qℓ(µi) in (2.49) or (2.50), ℓ = 1, 2.
8: end for
9: if ( Not all k desired eigenvalues are convergent and r < rmax ) then
10: Do the Schur-restarting in Lines 9–13 of Algorithm 2.1.
11: end if
12: until ( desire m eigenpairs are convergent or r≥ rmax )
13: Set λi = σ + µ−1i and ui = Vm2zi for i = 1, . . . , k.
14: if Q-SEP2 is solved then 15: ui ← eKu−1ui, i = 1, . . . , k.
16: end if
2.4.2
Computational costs
In this subsection, we compare the computational costs of the j-th Arnoldi step of Algorithm 2.1 for solving Q-SEPs (2.19) and R-SEPs (2.31), respectively. This is of general interest, because a comparison of the CPU time is sensible only if the number of outer iterations of Algorithm 2.2 or 2.3 is the same for each algorithm.
2.4 Arnoldi Method with Schur-restarting
Algorithm 2.3 Arnoldi method with Schur-restarting for solving REP in (2.20) Input: Coefficient matrices Mp, Kp and Ap, parameters c, α and β, σ: shift value,
tolR: tolerance for convergence, rmax: maximum number of Schur-restartings. Output: The desired eigenpairs (λi, pi) for i = 1, . . . , k.
1: Construct matrices fMp, eDp and eKp defined in (2.24), (2.25) and (2.26), respec-tively, 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 2.1 with A in
R-SEP1 or R-SEP2.
4: restart
5: Do the steps in Lines 3 and 4 of Algorithm 2.1. 6: for i = 1, . . . , m do 7: Compute ψ(µi) =| (σ + µ−1i )2 c2 |kMpk + kKpk + | (σ + µ−1i )2 α + (σ + µ−1i )β |kApk. 8: Check convergence by rℓ(µi) ψ(µi) < tolR with rℓ(µi) in (2.53) or (2.54), ℓ = 1, 2. 9: end for
10: if ( Not all k desired eigenvalues are convergent and r < rmax ) then
11: Do the Schur-restarting in Lines 9–13 of Algorithm 2.1.
12: end if
13: until ( desire m eigenpairs are convergent or r≥ rmax )
14: Set λi = σ + µ−1i and pi = Vm2zi for i = 1, . . . , k.
15: if R-SEP2 is solved then
16: pi ← eKp−1pi, i = 1, . . . , k.
2.4 Arnoldi Method with Schur-restarting
From (2.47), the unit vector vj+1 is generated by
Avj = j X
i=1
hjivi+ hj+1,jvj+1,
where hji = v∗iAvj for i = 1, . . . , j and hj+1,j = kAvj − j P i=1 hjivik2. For conve-nience, we let vj = h vj1 vj2 i
with vj1, vj2 ∈ Cn. The matrix-vector product Avj in Algorithm 2.2 for solving QEP (2.13) by Q-SEP1 (2.19.1) and Q-SEP2 (2.19.2) can be, respectively, represented by
B−1 u Auvj = −fMuvj2 e K−1 u (vj1− eDuvj2) and AuBu−1vj = −fMugu vj1− eDugu
with gu = eKu−1vj2. This implies that Algorithm 2.2 for Q-SEP1 and Q-SEP2 needs the same computational costs for generating the unit vector vj+1 for each j.
On the other hand, by letting vj = [v⊤j1, v⊤j2, vj3⊤]⊤ with vj1, vj2 ∈ Cn and vj3 ∈ Cℓ, the matrix-vector product Av
j in Algorithm 2.3 for solving REPs by R-SEP1 (2.31.1) and R-SEP2 (2.31.2) can be, respectively, represented by
B−1 p Apvj = Lpvj3− fMpvj2 e K−1 p (vj1− eDpvj2) ̺−1ϑv j3− ̺−1R⊤pvj2 and ApB −1 p vj = ̺−1L pvj3− fMpgp vj1− eDpgp ̺−1ϑv j3− R⊤pgp
with gp = eKp−1vj2. Consequently, the computational cost of ApBp−1vj needs an extra cost for the computation of ̺−1L
pvj3 compared to that Bp−1Apvj. The cost for generating the unit vector vj+1 by SEP1 is slightly cheaper than that by R-SEP2. We summarize the computational costs of generating vj+1 for by Q-SEP2 and R-SEP2 in Table 2.1.
2.5 Numerical Results
Q-SEP2 (2.19.2) R-SEP2 (2.31.2)
Solving linear system Keuxu = bu Kepxp = bp
Matrix-vector products Mfubu, eDubu Mfpbp, eDpbp, Lpcp, R⊤pc⊤p
Inner products j + 1 j + 1
Saxpy operators j + 1 j + 2
Scale-vector product 1 1
Table 2.1: Computational costs of the j-th Arnoldi step of Algorithm 2.1 for Q-SEP2 and R-Q-SEP2.
gp = eKp−1vj2 for j = 1, . . . , k can be saved in Gu ≡ [ eKu−1v12 · · · eKu−1vm2] and Gp ≡ [ eKp−1v12 · · · eKp−1vm2], respectively, so that the vectors u2, p2 in (2.50) and (2.54) can be computed by u2 = Guz and p2 = Gpz directly. Hence, it requires the same computational costs for computing u1, u2 in (2.49) and (2.50), as well as, p1, p2 in (2.53) and (2.54), respectively. Consequently, the computational costs of Q-SEP1 for the convergence test in Algorithm 2.2 need one extra matrix-vector product eKuvm+1,2 than those of Q-SEP2 in computing ζ1 and ζ2. Similarly, the computational costs of R-SEP1 for the convergence test in Algorithm 2.3 need one extra matrix-vector product eKpvm+1,2 than those of R-SEP2 in computing ξ1 and ξ2. Therefore, we conclude that Algorithm 2.2 for Q-SEP1 and Q-SEP2, as well as, Algorithm 2.3 for R-SEP1 and R-SEP2, respectively, almost have the same computational costs provided that they have the same outer iterations.
2.5
Numerical Results
We conduct numerical experiments to evaluate performance and accuracy of the eigenvalue solvers described in Section 2.4. To distinguish between various eigenvalue problems, we use notations Q1, Q2, R1 and R2 defined as follows:
2.5 Numerical Results
• Q2: Applying Algorithm 2.2 to solve the QEP (2.13) with Q-SEP2 in (2.19.2). • R1: Applying Algorithm 2.3 to solve the REP (2.20) with R-SEP1 in (2.31.1). • R2: Applying Algorithm 2.3 to solve the REP (2.20) with R-SEP2 in (2.31.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. We apply Algorithms 2.2 and 2.3 to solve the following examples arising in fluid-solid systems. The order m of Arnoldi decomposition in Line 3 of Algorithm 2.1 is set m = 40, the maximum number rmax of Schur-restartings is set rmax= 15 and the number of desired eigenpairs is k = 10. The relative residuals of approximate eigenpairs (λi, ui) and (λi, pi) computed by Q1 and Q2, as well as, R1 and R2 are, respectively, defined by
kQ(λi)uik ϕ(λi)kuik
and kR(λi)pik ψ(λi)kpik ,
where ϕ(λi) and ψ(λi) are given in Algorithm 2.2 and 2.3, respectively. Tolerances for relative residuals of QEPs and REPs are chosen by tolQ = tolR= 5× 10−15. The linear systems in Algorithms 2.2 and 2.3 are solved by LU-factorization with the shift value σ =−25 + 600πi. Fronbenius norm for matrices and 2-norm for vectors are used.
Example 2.1. [6] We take the geometrical data: the domain Ω = [0m, 1m] × [−0.75m, 0m], ΓA= [0m, 1m]× {0m} given in Figure 2.1(i) and the following physical data: ρ = 1kg/m3, c = 340 m/s, α = 5× 104 N/m3, and β = 200 Ns/m3.
The rectangular domain Ω is uniformly partitioned into nℓ by nw rectangles and each rectangle is further refined into two triangles, see Figure 2.1(ii). The dimen-sions of coefficient matrices in QEP (2.13) and REP (2.20) are (3nℓ− 1) × nw and
2.5 Numerical Results 0.75 m b= 1.00 m a= absorbing wall ΓA R Γ n rigid walls
(i) Fluid in a cavity with one absorbing wall. (ii) Initial mesh.
Figure 2.1: Fluid in a cavity with one absorbing wall and initial mesh
(nℓ+ 1)×(nw+ 1), respectively. Figure 2.2 plots the analytic solutions of the desired eigenvalues λ1, . . . , λ10 of (2.5)–(2.8) (see [6]) with the lowest positive vibration fre-quencies satisfying 0 < Im(λi)
2π < 600Hz.
Convergence test: We first demonstrate convergence rates of Q2 and R2 while computing the desired eigenvalues in Figure 2.2. To measure the convergence rate, we run the test over the five successively refined meshes (See the first column of Table 2.2) 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, 2, 3,
where λ[i,j]for j = 1, . . . , 5 denote the approximate eigenvalues computed by Q2 and R2 corresponding to λi obtained from the meshes described in Table 2.2. The 5-th and the 6-th columns of Table 2.2 illustrate the quadratic convergence of rate[1,j] j = 1, 2, 3 for λ1 of QEP (2.13) and REP (2.20), respectively. In our numerical experiment, the convergence rate are always close to 2 for all desired eigenvalues, λi, i = 1, . . . , 10, computed by Q2 and R2 as well as Q1 and R1.
