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 Ω⊂ R^d(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.

The dynamic variables of our model problem are the fluid pressure P and the

2.2 The Model Problem

displacement field U, which satisfy ([8, 35])

ρ∂²U

∂t² + ∇P = 0 in Ω, (2.1)

P = −ρc²divU 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 : Ω → Cⁿ and p : Ω→ C, (u, p) 6= (0, 0) such that

ρλ²u + ∇p = 0 in Ω, (2.5)

p = −ρc² 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 ∈ L²(∂Ω) and v· n = 0 on Γ^R}.

Here we employ standard Sobolev spaces notation. For example, H(div, Ω) stands for the space of all L² vector functions v on Ω with L² 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

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 ∈ H¹(Ω) :={p ∈ L²(Ω) : ∇p ∈ L²(Ω)} such that 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 R^d, d = 2, 3 (see, for instance, [10, 50]). For simplicity we will consider only the two dimensional case. Let {T^h} be a regular family of triangulations of Ω indexed by h, the maximum diameter of the elements. Let

V^h :={v^h ∈ H(div, Ω) : v^h|^T ∈ P0^d⊕ P⁰x ∀ T ∈ T^h and vh· n = 0 on Γ^R} ⊂ V,

where d = 2 andP^kdenotes the set of polynomials of degree at most k. Thus locally v_h takes the form (a + sx, b + sy)^⊤. The discrete problem associated with (2.9) is : Find λ∈ C and u^h ∈ V^h, uh 6= 0, such that

λ² Z

Ω

ρuh· ¯v^h+ λ Z

Γ_A

βuh· n ¯v^h· n + Z

Γ_A

αuh· n ¯v^h· n + Z

Ω

ρc²divuh div¯vh= 0, ∀ v^h∈ V^h. (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

divuh = 0 on Ω and uh · n = 0 on ∂Ω.

2.2 The Model Problem

Since uh = (a + sx, b + sy)^⊤ on T ∈ T^h, 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 = (n₁, n₂)^⊤ be a unit normal to an edge e of T , so t = (−n2, n₁)^⊤ is a unit tangent vector to e. We see that

u_h· 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 N^j, 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

∇^⊥Sh :={∇^⊥Ψh : Ψh is continuous piecewise linear and vanishes on the boundary}.

Thus we have just shown the zero eigenspaceE⁰ is contained∇^⊥Sh and the opposite inclusion is also easily checked. Hence

E⁰ =∇^⊥Sh.

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,

2.2 The Model Problem

suppose PN j=1

cj∇^⊥Ψj = 0. Then this implies PN j=1

cjΨj = c for some constant c. Hence cj = c by the condition Ψj(Nk) = δkj. Consequently, c(P

jΨj − 1) = 0. But we know

PN j=1

Ψj 6= 1 due to the vanishing boundary condition. Thus c^j = 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 :={p^h ∈ H¹(Ω) : ph|^T ∈ P1 ∀ T ∈ T^h}.

This is the subspace of H¹(Ω) consisted of continuous piecewise linears. The alter-native discrete problem in terms of the approximate pressure field is: Find λ ∈ C and ph ∈ H^h such that

λ² c²

Z

Ω

phq¯h+ λ² α + λβ

Z

ΓA

ρphq¯h+ Z

Ω∇p^h· ∇¯q^h = 0 ∀ q^h ∈ H^h. (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, ∇ · u^h = 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}ⁿj=1¹

denote the cardinal basis of V^h, 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 · nⁱdς = δij. For uh ∈ V^h, we write uh =

n1

P

j=1

ujφj and denote u = [u^⊤₁,· · · , 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≡ (λ²Mu+ (α + λβ)Au+ Ku)u = 0, (2.13)

where Mu ≡ [Mij^u] and Ku ≡ [Kij^u] are mass and stiffness matrices, respectively, and Au ≡ [A^uij] is used to describe the effect of the absorbing wall. Here

M_ij^u = Z

Ω

ρφi· ¯φj, K_ij^u = Z

Ω

ρc² divφi div ¯φj, A^u_ij = 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

µ²Mfu+ µ eDu+ eKu

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

A^u =

and linearize (2.15) into the GEP

A^uϕ= 1

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⁻¹u A^u)ϕ = µ⁻¹ϕ and (A^uBu⁻¹)ψ = µ⁻¹ψ,

2.3 Linearization of Nonlinear Eigenvalue Problems

respectively, where ψ =B^uϕ. Therefore, from (2.17) and (2.18) we have

(Q-SEP1) B⁻¹u A^u

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}ⁿj=1² be a nodal basis of Hh. For ph ∈ H^h, we write ph = Ap ≡ [A^pij] describes the effect of the absorbing wall. Here,

M_ij^p =

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.

By applying the long division, the rational term in (2.22) can be simplified into the following

2.3 Linearization of Nonlinear Eigenvalue Problems

the REP in (2.20) can be reformulated as

µ²Mfp+ µ eDp+ eKp



p− µ²Lpq= 0. (2.29)

Using (2.28) and (2.29), we get the GEP

A^pϕ ≡ 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 (Bp⁻¹A^p)ϕ = µ⁻¹ϕand (A^pB⁻¹p )ψ =

2.3 Linearization of Nonlinear Eigenvalue Problems

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 (µ⁻¹₁ ,h

v1

u1

i) be an

approximate eigenpair of (2.19.1) and h

f11

f12

i be the associated residual vector. That

is,

2.3 Linearization of Nonlinear Eigenvalue Problems

i) be an approximate eigenpair of (2.19.2) andh

f21

f22

i

be the associated residual vector. That is,



Now, we derive residual bounds of approximate eigenpairs for REP (2.20) by us-ing linearizations R-SEP1 and R-SEP2, respectively. Let (µ⁻¹₁ , [s^⊤₁, p^⊤₁, q^⊤₁]^⊤) be

2.3 Linearization of Nonlinear Eigenvalue Problems

Substituting (2.36) into (2.34), s1 can be represented by

s₁ =−µ1Mfpp₁+ µ²₁ µ1̺⁻¹ϑ− 1
−1

2.3 Linearization of Nonlinear Eigenvalue Problems

Substituting (2.40) and (2.41) into (2.39), we have

µ²₂MfpKe_p⁻¹t₂+ µ2DepKe_p⁻¹t₂+ t2− µ²2 ϑ− ̺µ⁻¹2

2.3 Linearization of Nonlinear Eigenvalue Problems

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/m³, c = 340 m/s, α = 5× 10⁴ N/m³, and β = 200 Ns/m³. In addition, we choose σ =−25 + 600πi as the shift value. Then

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

h²

respectively. Hence, by (2.16) the infinity norm of eKu can be estimated by k eKuk^∞ ≈ kK^uk^∞ =O(ρc²) = O(10⁵). 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

h²

2.4 Arnoldi Method with Schur-restarting

respectively. Hence, by (2.26) we have that k eKpk∞ ≈ kK^pk∞ =O(1). If the eigenvalue λ is one of the desired eigenvalues in Figure 2.2, then with µ = λ−σ we have

4× 10⁷ <

  µ²

 β

σβ + α + 1 µ

−1< 3.1× 10¹⁰

and

0.001 <

 

µ² α²β (σβ + α)⁴

 β

σβ + α + 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 (µ₁+σ, p₁) of REP (2.20) by solving R-SEP1 is larger than that of (µ₂+ σ, p₂) 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 V_m+1 is generated to satisfy

AVm = VmHm+ hm+1,mv_m+1e^⊤_m, (2.43)

where Hm ∈ C^m×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_k^H U_ℓ^H



 (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ℓ



 + h^m+1,mv_m+1(e^⊤_m[Uk Uℓ]),

which implies that

A ˜Vk= ˜VkTk+ ˜v_k+1t^H_k, (2.45)

where ˜Vk ≡ V^mUk, ˜v_k+1 = vm+1 and t^H_k ≡ h^m+1,me^⊤_mUk.

Let Q1 be a Householder matrix with t^H_kQ1 = τ e^⊤_k. Then (2.45) can be rewritten as

A( ˜VkQ₁) = ( ˜VkQ₁)(Q^H₁ TkQ₁) + τ ˜v_k+1e^⊤_k. (2.46)

The matrix Q^H₁ TkQ1 can be reduced to a new Hessenberg matrix H_k⁺ by using Householder matrices Qi for i = 2, . . . , k− 1 with

Q^H_k−1· · · Q^H2 (Q^H₁ TkQ1)Q2· · · Qk−1 = H_k⁺ e^⊤_kQ2· · · Q^k−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

AV_k⁺ = V_k⁺H_k⁺+ τ v⁺_k+1e^⊤_k

is obtained where V_k⁺ := ˜VkQ₁· · · Qk−1, v_k+1⁺ = ˜v_k+1 = v_m+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 (µ⁻¹_i , zi) with Hkzi = µ⁻¹_i zi, i = 1, . . . , m and sorting Ritz values so that {(µ1, z₁), . . . , (µk, zk)} are wanted.

5: for i = 1, . . . , k do

6: Check convergence by |hm+1,m||e^⊤mzi| < tol^A.

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 t^H_k := hm+1,me^⊤_mUk.

11: Compute Householder transformation Q1 such that t^H_kQ1 = τ e^⊤_k.

12: Reduce Q^H₁ 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,kv_k+1e^⊤_k. (2.47)

14: end if

15: until ( desired k eigenpairs are convergent or r ≥ r^max )

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

iare partitioned with compatible sizes. Using

the first equation of (2.48), we can eliminate Vm1z in the second equation and get

kQ(λ)u¹k

ku1k = k(µ²Mfu+ µ eDu+ eKu)u1k

ku1k = |µ| |h^m+1,m|e^⊤_mz ζ₁

ku1k ≡ q¹(µ), (2.49)

where u1 = Vm2z, λ = µ + σ and ζ1 = kµvm+1,1+ eKuv_m+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

 (2.49) and q₂(µ) 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 (µ⁻¹, z) be a Ritz pair and satisfy Hmz = µ⁻¹z. From (2.43), and SEP1,

R-2.4 Arnoldi Method with Schur-restarting

SEP2 in (2.31) we have

 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 residuals while Algorithm 2.1 is applied to solve REPs (2.20).

Applying Algorithm 2.1 to solve QEPs (2.13) and REPs (2.20) are summarized

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) = (|σ + µ⁻¹i |²kM^uk + |α + (σ + µ⁻¹i )β|kA^uk + kK^uk).

7: Check convergence of QEP by qℓ(µi)

ϕ(µi) < tolQ

ϕ(µi) < tolQ

在文檔中 二次暨有理特徵值問題中高效能Arnoldi型態演算法 (頁 22-109)