SOLVING A STRUCTURED QUADRATIC EIGENVALUE PROBLEM
BY A STRUCTURE-PRESERVING DOUBLING ALGORITHM∗
CHUN-HUA GUO† AND WEN-WEI LIN‡
Abstract. In studying the vibration of fast trains, we encounter a palindromic quadratic eigen-value problem (QEP) (λ2AT +λQ + A)z = 0, where A, Q ∈ Cn×n and QT =Q. Moreover, the
matrixQ is block tridiagonal and block Toeplitz, and the matrix A has only one nonzero block in the upper-right corner. So most of the eigenvalues of the QEP are zero or infinity. In a linearization approach, one typically starts with deflating these known eigenvalues for the sake of efficiency. How-ever, this initial deflation process involves the inverses of two potentially ill-conditioned matrices. As a result, large error might be introduced into the data for the reduced problem. In this paper we propose using the solvent approach directly on the original QEP, without any deflation process. We apply a structure-preserving doubling algorithm to compute the stabilizing solution of the ma-trix equationX + ATX−1A = Q, whose existence is guaranteed by a result on the Wiener–Hopf factorization of rational matrix functions associated with semi-infinite block Toeplitz matrices and a generalization of Bendixson’s theorem to bounded linear operators on Hilbert spaces. The doubling algorithm is shown to be well defined and quadratically convergent. The complexity of the doubling algorithm is drastically reduced by using the Sherman–Morrison–Woodbury formula and the special structures of the problem. Once the stabilizing solution is obtained, all nonzero finite eigenvalues of the QEP can be found efficiently and with the automatic reciprocal relationship, while the known eigenvalues at zero or infinity remain intact.
Key words. palindromic quadratic eigenvalue problem, nonlinear matrix equation, stabilizing solution, structure-preserving, doubling algorithm
AMS subject classifications. 15A24, 65F15, 65F30 DOI. 10.1137/090763196
1. Introduction. In this paper we consider a quadratic eigenvalue problem
(QEP) occurring in the vibration analysis of rail tracks under excitation arising from high speed trains [14, 15, 17]. This problem has provided much of the motivation for the study of palindromic polynomial eigenvalue problems in [22] and subsequent papers [5, 16, 18, 19, 20, 23]. Yet the problem itself has not been solved satisfactorily. The standard approach for solving a QEP is to use a proper linearization and solve a generalized eigenvalue problem of twice the size. Another approach for solving a QEP is through finding a solution (called a solvent) of a related matrix equation. This solvent approach has been explored in [6] and more recently in [13] and [26]. The difficulty associated with the solvent approach is obvious. It is possible that the related matrix equation does not have a solvent. Even if a solvent exists, the computation of a solvent can still be a difficult task. As a result, the solvent approach can outperform the linearization approach only for special types of QEPs [9, 11].
So far every method for the special QEP here starts with the linearization ap-proach. For the sake of efficiency, a deflation process is used in the beginning. This
∗Received by the editors June 26, 2009; accepted for publication (in revised form) August 31,
2010; published electronically October 14, 2010.
http://www.siam.org/journals/simax/31-5/76319.html
†Department of Mathematics and Statistics, University of Regina, Regina, SK S4S 0A2, Canada
([email protected]). The work of this author was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada and by the National Center for Theoretical Sciences in Taiwan.
‡Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan
([email protected]). The work of this author was partially supported by the National Science Council and the National Center for Theoretical Sciences in Taiwan.
2784
initial deflation process involves the sucessive application of the inverses of two po-tentially ill-conditioned matrices (see [5], for example). As a result, large error might be introduced into the data for the reduced problem. Several methods have been proposed recently to solve the reduced QEP. In particular, the solvent approach is used in [5]. However, some major issues associated with the solvent approach remain unsolved in [5]. Another efficient method is proposed and compared to two other methods in [16]. These methods continue to use the linearization approach for the reduced QEP. The computational work for all these different methods are dominated by that of the same deflation process. The accuracy of the computed solution is thus the main issue here.
In this paper we will see that the QEP arising in the study of high speed trains is very amenable for the solvent approach if used directly on the original QEP, without any deflation process. At first sight, the solvent approach applied to the original QEP would also be very expensive. In this paper we will show that the solvent approach can be implemented to have a complexity roughly the same as that for other efficient methods using the linearization approach and the initial deflation process. Numerical experiments show that our solvent approach, applied to the original QEP, produces better accuracy in the computed results.
2. The quadratic eigenvalue problem. The vibration analysis of rail tracks
can be performed through a finite element model, in which we generate [5] two real symmetric matrices M and K of the form
M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M0 MT 1 0 · · · 0 M1 M1 M0 MT 1 0 0 0 . .. . .. . .. . .. ... .. . 0 . .. . .. . .. 0 0 . .. M1 M0 M1T MT 1 0 · · · 0 M1 M0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ m×m , (2.1) K = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ K0 KT 1 0 · · · 0 K1 K1 K0 KT 1 0 0 0 . .. . .. . .. . .. ... .. . 0 . .. . .. . .. 0 0 . .. K1 K0 KT 1 KT 1 0 · · · 0 K1 K0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ m×m , (2.2)
where each block in M and K is q× q. So M, K ∈ Rn×nwith n = mq. The matrices M and K are thus block Toeplitz (actually block circulant). This special structure is not used in [5], and the notation used there for M and K is more general. A matrix D (the damping matrix) is then taken to be a positive linear combination of M and
K. That is, D = c1M + c2K with c1, c2 > 0. So D has the same structure as M
and K. We write M = Mt+ Mc+ McT, where Mtis the block tridiagonal part of M
and Mc is the matrix with M1 in the upper-right corner and zero blocks elsewhere. Similarly, we have K = Kt+ Kc+ KcT and D = Dt+ Dc+ DTc. We also denote by
ω > 0 the frequency of the external excitation force.
For the vibration analysis, one needs to solve the palindromic QEP [5]
(2.3) P (λ)z = 0, z = 0,
with
(2.4) P (λ) = λ2AT + λQ + A,
where
(2.5) Q = Kt+ iωDt− ω2Mt
with i =√−1 (so QT = Q) and
(2.6) A = Kc+ iωDc− ω2Mc.
The set of eigenvalues of the quadratic P (λ) demonstrates a “symplectic” behavior (i.e., a symmetry with respect to the unit circle, which is denoted by T throughout this paper). More precisely, a number ξ is an eigenvalue of the quadratic P (λ) if and only if ξ−1 is so, and they have the same partial multiplicities (see [22, Theorem 2.2]). Let 2 be the usual Hilbert space of all square summable sequence of complex numbers, and let q2be the Cartesian product of q copies of 2. The infinite matrices
(2.7) TM = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ M0 MT 1 M1 M0 MT 1 M1 M0 . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, TK = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ K0 KT 1 K1 K0 KT 1 K1 K0 . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ are then seen to be inB(q2), the set of all bounded linear operators on q2. They are also self-adjoint operators inB(q2). It is well known that the spectrum of a self-adjoint operator is real [27].
By the way the matrices M and K are generated in the finite element model, we know that TK is positive semidefinite, written TK ≥ 0, in the sense that TKf, f ≥ 0
for all f ∈ q2. We also know that TM ≥ I (i.e., TM−I ≥ 0) for the identity operator I
and some > 0. These properties on TKand TM can also be verified independently (in case significant errors are introduced in setting up the matrices M0, M1, K0, K1). As noted in [7], TK≥ 0 if and only if ψK(λ) = K0+ λK1+ λ−1K1T is positive semidefinite on T. So TM ≥ I for some > 0 if and only if ψM(λ) = M0+ λM1+ λ−1M1T is positive definite onT. The latter holds if and only if the matrix equation
(2.8) X + M1TX−1M1= M0
has a positive definite solution X with ρ(X−1M1) < 1 (see [7]), where ρ(·) denotes the spectral radius. The equation (2.8), where M0is symmetric positive definite, has been well studied (see [4, 7, 9, 10, 21, 25, 30]). In particular, instead of checking that ψM(λ) is positive definite onT, one can attempt to find the maximal positive definite solution of (2.8) by the cyclic reduction method in [25] or the doubling algorithm in [21]. These methods are very efficient, and the computational work involved is only a small fraction of that for solving the QEP, which involves mq× mq matrices while the matrices in (2.8) are q× q. Recall that D = c1M + c2K for c1, c2> 0. So TD≥ c1I, where (2.9) TD= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ D0 DT 1 D1 D0 DT 1 D1 D0 . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, D0= c1M0+ c2K0, D1= c1M1+ c2K1.
From now on we assume that TD ≥ ηI for some η > 0, which can be verified as we
have described for TM.
3. Theoretical results for the solvent approach. We first show that the
QEP does not have any eigenvalues on T. The following result, due essentially to Bendixson [2], can be found in [28].
Lemma 3.1 (Bendixson’s theorem). Let X and Y be any k× k Hermitian
ma-trices. Suppose that the eigenvalues of X are λ1≤ λ2≤ · · · ≤ λk and the eigenvalues
of Y are μ1≤ μ2 ≤ · · · ≤ μk. Then the eigenvalues of X + iY are contained in the
rectangle [λ1, λk]× [μ1, μk] in the complex plane.
Theorem 3.2. The quadratic P (λ) in (2.4) has no eigenvalues on T.
Proof. The quadratic P (λ) has eigenvalues on T if and only if det P (λ) = 0
for some λ ∈ T or equivalently det(λAT + λ−1A + Q) = 0 for some λ ∈ T. This is impossible since we can show that for each fixed λ ∈ T all eigenvalues of the matrix λAT + λ−1A + Q have positive imaginary parts. In fact, λAT + λ−1A + Q = X(λ) + iY (λ) with Hermitian matrices
X(λ) = Kt− ω2Mt+ λ−1(Kc− ω2Mc) + λ(Kc− ω2Mc)T,
Y (λ) = ω(Dt+ λ−1Dc+ λDTc).
By Bendixson’s theorem, we need only to show that Y (λ) > 0 on T or equivalently
TD≥ I for some > 0, where
(3.1) TD= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Dt DcT Dc Dt DTc Dc Dt . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
The latter is true since TD is precisely (a partition of) TDin (2.9).
We now consider the matrix equation
(3.2) X + ATX−1A = Q,
where Q and A are given in (2.5) and (2.6). Suppose X is a solution of (3.2). Then we have the factorization for the quadratic P (λ) in (2.4):
λ2AT + λQ + A = (λAT+ X)X−1(λX + A).
So the eigenvalues of the quadratic P (λ) are a collection of the eigenvalues of the pencils λAT + X and λX + A. We have already shown that P (λ) has no eigenvalues on T. Suppose a solution X of (3.2) can be found such that the eigenvalues of the pencil λX + A (or equivalently the eigenvalues of the matrix−X−1A) are inside T. Then the remaining eigenvalues of P (λ) are obtained by taking the reciprocals of these eigenvalues. Such a solution X is called a stabilizing solution of (3.2). In this process, the known eigenvalues of P (λ) at zero or infinity remain intact, regardless of the accuracy of the computed X.
There are two advantages of the solvent approach over the linearization approach. First, in the linearization approach, a deflation procedure is used for the sake of effi-ciency, which involves the inverses of two potentially ill-conditioned matrices. When the QEP is reduced to a smaller QEP (even if in a structure-preserving manner), the input data obtained in the smaller QEP could be significantly different from the true data. In the solvent approach, the ill-conditioning of those matrices may also affect the accuracy of the solution X computed by some efficient method, but we can always use Newton’s method as a correction method afterward, as in [10]. Second, in the linearization approach, the eigenvalues of the smaller QEP range in modulus from to −1, where is close to 0, while in the solvent approach the eigenvalues of λX + A range in modulus from to 1. The situation in the solvent approach is easier to handle, and the symplectic structure of the eigenvalues of P (λ) is preserved automatically.
The success of the solvent approach hinges on the existence of a stabilizing so-lution of (3.2) and an efficient method for its computation. In this section we prove the existence of a stabilizing solution. In the next section we show that a doubling algorithm can be used to compute it efficiently.
We start with a generalization of Bendixson’s theorem to bounded linear operators in Hilbert spaces. It seems that such a generalization has not been given before, although special cases of this are being proved in recent literature.
Lemma 3.3 (generalization of Bendixson’s theorem). Let B and C be self-adjoint bounded linear operators on a Hilbert space. Suppose that the spectrum of B is
con-tained in [u1, u2] and the spectrum of C is contained in [v1, v2]. Then the spectrum of
B + iC is contained in the rectangle [u1, u2]× [v1, v2] in the complex plane.
Proof. Some special cases have been proved in the literature. For example, it is proved in [1] (see Corollary 4 there) that B + iC is invertible if B is invertible and positive definite (or equivalently [27, section 7.4, Corollary 2] if min σ(B) ≥ for some > 0). Also, it is proved in [12] (see the proof of Lemma 3.1 there) that B +iC is invertible if C ≥ I for some > 0 (or equivalently [27, section 7.4, Corollary 2] if min σ(C) ≥ ). Note that the result in [12] follows from the result in [1] by a multiplication with i. The general statement in Lemma 3.3 can also be proved quickly using the special case proved in [1]. We need only to prove that each point a + bi in the spectrum of B + iC satisfies a≥ u1(the rest can be proved by multiplying B + iC
with−1 or i). We may assume u1> 0 by shifting B to B + ηI for some η > 0. Since
σ(B + iC) is a compact set [27, Theorem 5.14], the distance between the imaginary axis and σ(B + iC) is attained for a point a∗+ b∗i in σ(B + iC). We need to show a∗ ≥ u1. Suppose a∗< u1. Then B− a∗I ≥ (u1− a∗)I with u1− a∗> 0, and thus
(B + iC)− (a∗+ b∗i)I = (B − a∗I) + i(C − b∗I) is invertible by [1, Corollary 4]. This is a contradiction since a∗+ b∗i is in σ(B + iC).
To prove the existence of a stabilizing solution of (3.2), we consider the semi-infinite block Toeplitz matrix
(3.3) T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Q AT A Q AT A Q . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦.
Associated with T is the rational matrix function φ(λ) = λA + Q + λ−1AT. It is clear that T is inB(n2), and we will show that T is invertible.
By (2.5) and (2.6) we have T = B + iC with B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Bt BcT Bc Bt BTc Bc Bt . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
and C = ω TD, where Bt = Kt− ω2Mt, Bc = Kc− ω2Mc, and TD is given in (3.1). Note that B and C are self-adjoint operators inB(n2). Since TD is a partition of TD
and TD≥ ηI for some η > 0, we have C ≥ ωηI, and thus min σ(C) ≥ ωη. By Lemma
3.3, 0 /∈ σ(B + iC). So T = B + iC is invertible. We can now prove the following result.
Theorem 3.4. The equation (3.2) has a unique stabilizing solution, and the solution is complex symmetric. Moreover, the dual equation of (3.2)
(3.4) X + A X−1AT = Q
also has a unique stabilizing solution, and the solution is complex symmetric.
Proof. Since T is invertible, we know from a result on linear operators (see
[8, Chapter XXIV, Theorem 4.1] and [24]) that φ(λ) has the so-called Wiener–Hopf factorization
(3.5) φ(λ) = (I − λ−1L)D(I − λU)
with D invertible, ρ(L) < 1, and ρ(U ) < 1. From (3.5) we see that
A = −DU, AT =−LD, Q = D + LDU.
Thus
(3.6) D + ATD−1A = Q
with ρ(D−1A) < 1 and ρ(ATD−1) < 1. So D is a stabilizing solution of (3.2). We will see in the next section that the pencil N0− λL0 defined by (4.1) has exactly
n eigenvalues inside T and that for any stabilizing solution Xs of (3.2) the column
space of [I XsT]T is, by (4.3), the (necessarily unique) deflating subspace of the pencil N0− λL0corresponding to its n eigenvalues insideT. It follows that (3.2) has exactly one stabilizing solution. Now taking transpose in (3.6) gives DT+ AT(DT)−1A = Q. Note that ρ((DT)−1A) = ρ(ATD−1) < 1. So DT is also a stabilizing solution of (3.2). The uniqueness of stabilizing solutions implies that DT = D.
The statements about the dual equation can be proved in a similar way. The only difference is that we now need to show that the self-adjoint operator inB(n2),
(3.7) T D= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Dt Dc DT c Dt Dc DT c Dt . .. . .. . .. ⎤ ⎥ ⎥ ⎥ ⎥ ⎦,
is such that TD≥ I for some > 0. This is true since TDis related to TD in (3.1) by
TD= W TDW,
where W = ⎡ ⎢ ⎣ V V . .. ⎤ ⎥ ⎦ , V = ⎡ ⎣ · · Iq · Iq ⎤ ⎦ m×m ,
and Iq is the q×q identity matrix. Thus for any f ∈ n2, TDf, f = TD(W f ), W f ≥
η W f 2= η f 2. So TD≥ ηI.
4. Computation of the stabilizing solution. A doubling algorithm has been
studied in [21] for (3.2) with a real A and a real symmetric positive definite Q. In our case, A is complex and Q is complex symmetric. However, the more general presentation in [4] can be used directly.
Let (4.1) N0= A 0 Q −I , L0= 0 I AT 0 .
Then the pencil N0− λL0 is a linearization of the T -palindromic polynomial λ2AT−
λQ + A. It is easy to verify that the pencil N0− λL0 is T -symplectic, i.e.,
N0JNT 0 = L0JLT0 for J = 0 I −I 0 . We can define the sequences{Nk} and {Lk}, where
(4.2) Nk= A k 0 Qk −I , Lk= −Pk I AT k 0 , by the following doubling algorithm [4] if no breakdown occurs.
Algorithm 4.1. Let A0= A, Q0= Q, P0= 0. For k = 0, 1, . . . , compute
Ak+1= Ak(Qk− Pk)−1Ak,
Qk+1= Qk− ATk(Qk− Pk)−1Ak,
Pk+1= Pk+ Ak(Qk− Pk)−1ATk.
We now show that this algorithm will not break down, and Qkconverges
quadrat-ically to the stabilizing solution of (3.2).
Theorem 4.1. Let A and Q be given by (2.6) and (2.5). Let Xsbe the stabilizing
solution of (3.2) and Xs be the stabilizing solution of the dual equation (3.4). Then
(a) the sequences{Ak}, {Qk}, {Pk} in Algorithm 4.1 are well defined, and Qk and
Pk are complex symmetric.
(b) Qk converges to Xs quadratically, Ak converges to 0 quadratically, Q− Pk
converges to Xs quadratically, with
lim sup
k→∞
2k Q
k− Xs ≤ (ρ(Xs−1A))2, lim sup
k→∞ 2k A k ≤ ρ(Xs−1A), lim sup k→∞ 2k Q − Pk− Xs ≤ (ρ(Xs−1A))2,
where · is any matrix norm.
Proof. Let Tk be the leading principal block k× k submatrix of T in (3.3) and
write Tk = Bk+ iCk, where Bk and Ck are Hermitian. So
Ck= ω ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Dt DTc Dc Dt . .. . .. . .. DT c Dc Dt ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ k×k .
Since TD ≥ ηI, we have TDf, f ≥ η f 2 for all f ∈ q2. Taking f = g0 with
g ∈ Ckmq, we know that for all g ∈ Ckmq,
Ckg, g = ωTDf, f ≥ ωη f 2= ωη g 2.
Thus Ck is positive definite for each k≥ 1. It then follows from Bendixson’s theorem
that Tk is invertible for each k≥ 1.
Let Wk= Qk− Pk in Algorithm 4.1. Then the sequence{Wk} satisfies
Wk+1= Wk− ATkWk−1Ak− AkWk−1ATk, W0= Q.
It follows from [3, Theorem 13; see also equation (9)] that Wk is nonsingular for each
k ≥ 0. The sequences {Ak}, {Qk}, {Pk} are then well defined. It is easy to see by
induction that Qkand Pkare complex symmetric since Q is complex symmetric. This proves (a).
To prove (b), we start with the easily verified relation
(4.3) N0 I Xs = L0 I Xs X−1 s A.
From the discussions in [4] we have for each k≥ 0
(4.4) Nk I Xs = Lk I Xs (Xs−1A)2k. Substituting (4.2) into (4.4) yields
(4.5) Ak = (Xs− Pk)(Xs−1A)2 k , Qk− Xs= ATk(Xs−1A)2 k . Similarly we have N0 I Xs = L0 I Xs X−1 s AT, where N0= AT 0 Q −I , L0= 0 I A 0 .
The pencil N0−λL0is a linearization of λ2A−λQ+AT, which has the same eigenvalues
as λ2AT− λQ + A. It follows that Xs−1AT and Xs−1A have the same eigenvalues, and
thus ρ( Xs−1AT) = ρ(Xs−1A). For each k ≥ 0 we now have
(4.6) Nk I Xs = Lk I Xs ( Xs−1AT)2k,
where Nk= Ak 0 Qk −I , Lk= − Pk I AT k 0
are defined by Algorithm 4.1, starting with A0 = AT, Q0= Q, P0 = 0. It is easy to prove by induction that for all k≥ 0
(4.7) Ak= ATk, Pk= Q− Qk, Qk= Q− Pk.
Indeed, assuming (4.7) for k, we have
Qk+1= Qk− ATk( Qk− Pk)−1Ak= Q− Pk− Ak(Qk− Pk)−1ATk = Q− Pk+1,
and similarly we have Ak+1= ATk+1and Pk+1= Q− Qk+1. By (4.6) and (4.7) we now have
(4.8) ATk = ( Xs− Pk)( Xs−1AT)2k, Qk− Xs= Ak( Xs−1AT)2k. By (4.5), (4.8), and (4.7) we have Qk− Xs= ATk(Xs−1A)2k = ( Xs− Pk)( Xs−1AT)2 k (Xs−1A)2k = (Qk− Xs+ (Xs+ Xs− Q))( Xs−1AT)2 k (Xs−1A)2k, from which we obtain
(4.9) (Qk− Xs)(I− ( Xs−1AT)2k(Xs−1A)2k) = (Xs+ Xs− Q)( Xs−1AT)2 k (Xs−1A)2k. It follows that lim sup k→∞ 2k Q k− Xs ≤ ρ( Xs−1AT)ρ(Xs−1A) = (ρ(Xs−1A))2< 1.
So Qk converges to Xs quadratically. Then we know{ Pk} is bounded and have by
the first equation in (4.8) that lim sup
k→∞
2k A
k ≤ ρ(Xs−1A) < 1.
So Ak converges to 0 quadratically. By the second equations in (4.7) and (4.8) we get
lim sup
k→∞
2k
(Q − Pk)− Xs ≤ (ρ(Xs−1A))2< 1.
So Q− Pk converges to Xs quadratically. This completes the proof of (b).
Algorithm 4.1 is said to be structure-preserving since for each k≥ 0, Nk and Lk
have the structures given in (4.2), and the pencil Nk− λLk is T -symplectic.
The complexity of Algorithm 4.1 can be reduced drastically by using the special structure of the matrix A given by (2.6). Write Qk= Q− Rk. Then it is easy to see
by induction that the matrices Ak, Rk, and Pk have the special forms
Ak = ⎡ ⎢ ⎢ ⎣ Ek 0 · · · 0 ⎤ ⎥ ⎥ ⎦ , Rk = ⎡ ⎢ ⎢ ⎢ ⎣ 0 . .. 0 Fk ⎤ ⎥ ⎥ ⎥ ⎦, Pk = ⎡ ⎢ ⎢ ⎢ ⎣ Gk 0 . .. 0 ⎤ ⎥ ⎥ ⎥ ⎦,
where the q×q matrices Ek, Fk, and Gk can be determined by the following simplified algorithm, in which (4.10) Q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ H0 HT 1 H1 H0 . .. . .. . .. HT 1 H1 H0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ m×m is given by (2.5), with H0= K0+ iωD0− ω2M0, H1= K1+ iωD1− ω2M1. Algorithm 4.2. Let E0= H1, F0= 0, G0= 0. For k = 0, 1, . . . , compute ⎡ ⎢ ⎢ ⎢ ⎣ Sk,1 Sk,2 .. . Sk,m ⎤ ⎥ ⎥ ⎥ ⎦= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Q − ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Gk 0 . .. 0 Fk ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ −1⎡ ⎢ ⎢ ⎢ ⎣ Ek 0 .. . 0 ⎤ ⎥ ⎥ ⎥ ⎦, (4.11) ⎡ ⎢ ⎢ ⎢ ⎣ Tk,1 Tk,2 .. . Tk,m ⎤ ⎥ ⎥ ⎥ ⎦= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝Q − ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Gk 0 . .. 0 Fk ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ −1⎡ ⎢ ⎢ ⎢ ⎣ 0 .. . 0 ET k ⎤ ⎥ ⎥ ⎥ ⎦, (4.12)
where all matrix blocks on the left side of (4.11) and (4.12) are q×q, and then compute
(4.13) Ek+1= EkSk,m, Fk+1= Fk+ EkTSk,1, Gk+1= Gk+ EkTk,m.
The main task of Algorithm 4.2 is to solve the large sparse linear systems in (4.11) and (4.12). We rewrite the common matrix in (4.11) and (4.12) as
(4.14) Q − ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Gk 0 . .. 0 Fk ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = Q− BkCkT with (4.15) Bk = GT k 0 · · · 0 0 · · · 0 FkT T , CT k = Iq 0 · · · 0 0 · · · 0 Iq , and solve the linear systems by the Sherman–Morrison–Woodbury formula (4.16) (Q− BkCkT)−1= Q−1+ Q−1Bk(I2q− CkTQ−1Bk)−1CkTQ−1.
Let Q = UHR be a qr-factorization of Q, where U is unitary and R is upper triangular. Since QT = Q, a linear system QX = B can be solved by
(4.17) X = R−1UB or X = UTR−TB.
Write Q = [Qij]mi,j=1 with
(4.18) ⎧ ⎨ ⎩ Qii= H0, i = 1, · · · , m, Qi+1,i= QTi,i+1= H1, i = 1, · · · , m − 1, Qij = 0q, |i − j| > 1,
where 0qis the q×q zero matrix. Let Ukdenote the set of all k×k unitary matrices and Δk the set of all k× k upper triangular matrices. The following algorithm computes the qr-factorization of Q in a sparse way.
Algorithm 4.3. sqr-factorization: [U,R]= sqr(Q). Input: Q as in (4.18).
Output: U =U(i,i+1)∈ U2q, i = 1 : m − 1, U(m,m)∈ Uq and R∈ Δn.
Set R← 0n;
For i = 1 : m− 1,
compute the qr-factorization
Qii Qi+1,i → (U(i,i+1))H Qii 0 ,
where U(i,i+1)∈ U2q and the new Qii is in Δq;
compute Qij Qi+1,j ← U(i,i+1) Qij Qi+1,j , j = i + 1 : min{i + 2, m};
Compute the qr-factorization Qmm→ (U(m,m))HQm,m,
where U(m,m)∈ Uq and the new Qmm is in Δq;
Rij ← Qij, i = 1 : m, j = i : min{i + 2, m}.
The above algorithm gives the qr-factorization Q = UHR, where the unitary matrix U = [Ui,j]mi,j=1, with Ui,j ∈ Cq×q, is given in a sparse factored form. More precisely, U = U(m,m) 1i=m−1U (i,i+1) with U(i,i+1) and U(m,m) being the extensions of U(i,i+1) and U(m,m), respectively, by adding appropriate 1’s on the diagonal. We now use the sqr-factorization of Q to solve the linear system QX = B with B =
[Iq, 0q, . . . , 0q]T. Note that Iq appears in the top position in B. In this process,
Ui,1 (i = 1 : m) are obtained explicitly for later use. Algorithm 4.4. [X1, Xm, U1:m,1] = Solt(U, R). Input: The output from Algorithm 4.3.
Output: The first and last submatrices of the solution X = [Xi]mi=1∈ Cn×q
with Xi∈ Cq×q for the linear system QX = [Iq, 0q, . . . , 0q]T and the first
block column of U . Set B1← Iq; For i = 1 : m − 1, compute Bi Bi+1 ← U(i,i+1)(1 : 2q, 1 : q) Bi; Compute Bm← U(m,m)Bm; Set Ui,1← Bi, i = 1 : m; For i = m :−1 : 1, compute Xi= R−1ii Bi− min{i+2,m}j=i+1 RijXj .
For the linear system QX = B with B = [0q, . . . , 0q, Iq]T, it is possible to compute
X1 and Xmdirectly without computing Xk(k = 2, . . . , m− 1).
Algorithm 4.5. [X1, Xm] = Solb(U, R, Um,1).
Input: The output from Algorithm 4.3 and Um,1 from Algorithm 4.4.
Output: The first and the last submatrices of the solution X = [Xi]mi=1∈ Cn×q
with Xi∈ Cq×q for the linear system QX = [0q, . . . , 0q, Iq]T.
Set Bm−1 Bm ← U(m−1,m)(1 : 2q, q + 1 : 2q);
Compute Bm← U(m,m)Bm; Xm← R−1m,mBm (by the first equation in (4.17));
Compute X1← Um,1T R−Tm,m (by the second equation in (4.17)).
The following algorithm gives a more detailed implementation of Algorithm 4.2 and computes the stabilizing solutions of (3.2) and (3.4) by Theorem 4.1.
Algorithm 4.6. Computation of Xs and Xs.
Input: H0,H1 ∈ Cq×q, tolerance τ .
Output: The solutions Xs∈ Cn×n for (3.2) and Xs∈ Cn×n for (3.4).
Take Q in (4.18), E0= H1, F0= 0, G0= 0; Call [U, R] = sqr(Q); [Y1, Ym, U1:m,1] = Solt(U, R); [Z1, Zm] = Solb(U, R, Um,1); For k = 0, 1, . . . Xk,1= [Y1Gk, Z1Fk], Xk,m= [YmGk, ZmFk]; [Xk,1f , Xk,mf ] = [Y1Ek, YmEk]; [Xk,1g , Xk,mg ] = [Z1ETk, ZmEkT];
Sk,i= Xk,if + Xk,i
I2q− !X k,1 Xk,m " −1#Xf k,1 Xk,mf $ , i = 1, m, Tk,m= Xk,mg + Xk,m I2q− !X k,1 Xk,m " −1#Xg k,1 Xk,mg $ ; Ek+1= EkSk,m, Fk+1= Fk+ EkTSk,1, Gk+1= Gk+ EkTk,m; If EkTSk,1 ≤ τ Fk and EkTk,m ≤ τ Gk , then Xs← Q, Xs(n: n, n : n)← H0− Fk+1, Xs← Q, Xs(1 : q, 1 : q)← H0− Gk+1,
where n= (m− 1)q + 1, and stop.
5. Numerical results. The sqr-factorization in Algorithm 4.3 requires about
86
3mq3flops. The linear system solvers in Algorithms 4.4 and 4.5 require 9mq3and 4q3 flops, respectively. Each iteration of the doubling algorithm in Algorithm 4.6 requires about 1543 q3 flops. Algorithm 4.6 is efficient since no more than 10 iterations are typically needed for convergence. For large q and m the total computational work for Algorithm 4.6 is thus roughly 1133 mq3, assuming that the number of iterations for the doubling algorithm is bounded independent of q and m. We note that Algorithms 4.3, 4.4, and 4.5 presented in this paper can also be used in the initial deflation procedure [5] for the linearization approach. So the deflation procedure can be completed in about 1133 mq3flops as well.
In this section we present numerical results to illustrate the efficiency and accu-racy of the solvent approach for computing the eigenpairs of the QEP (2.3), through computing the solvent Xs by Algorithm 4.6.
We first explain how the eigenpairs can be computed after the solvent Xs is
obtained. By Algorithm 4.6 we see that
(5.1) Q − Xs= [0, . . . , 0, Iq]T[0, . . . , 0, F∞],
where F∞ = limk→∞Fk. Write A =
0n−q H1t 0 0q
, where H1t = [H1T, 0, . . . , 0]T. Ap-plying U (given implicitly in a sparse factored form) in Algorithm 4.3 to A and Xs,
respectively, we have (5.2) UA = 0n−q H 1t 0 Φ1 , UXs= X1 X2 0 Φ2 ,
where X1= R(1 : n−q, 1 : n−q) and X2(1 : n−3q, 1 : q) = 0. From the factorization
P (λ) = (λAT+ X
s)Xs−1(λXs+ A), the nonzero stable eigenpairs (λs, zs) of P (λ) are
those of λXs+ A and can be computed by
Φ1zs,2=−λsΦ2zs,2, (5.3) zs,1=−X1−1(X2zs,2+ λ−1s H 1tzs,2), zs= zs,1 zs,2 (5.4)
for s = 1, . . . , q. Recall that the first block column of U is known from Algorithm 4.4. So Φ1= Um,1H1and H1t= U1:m−1,1H1.
If we are only interested in the eigenvalues, then we can find all nonzero stable eigenvalues from (5.3) and get all finite unstable eigenvalues by taking the reciprocals of the nonzero stable ones. The cost is O(q3) flops. Eigenvectors corresponding to nonzero stable eigenvalues can be found from (5.4) with a cost of 7mq3 flops, noting that X1 is block 3-banded upper triangular and that H1tzs,2= U1:m−1,1(H1zs,2).
Some further work is required if the eigenvectors corresponding to finite unstable eigenvalues are also needed. We first compute all left eigenvectors of λΦ2+ Φ1by
(5.5) yTsΦ1=−λsysTΦ2,
for s = 1, . . . , q, at a cost of O(q3) flops. The finite unstable eigenpairs (λu, zu) of
P (λ) satisfy (5.6) P (λu)zu≡ P ! 1 λs " zu= λ12 s % AT+ λ sXs&Xs−1(Xs+ λsA) zu= 0.
From (5.2) and (5.5) it follows that (5.7) %AT + λsXs&UT 0 ys = ! 0 0 H1 ΦT1 + λsX1T 0 λsX2T λsΦT2 " 0 ys = 0. From (5.6) and (5.2) the eigenvector zu corresponding to λu = λ−1s can be found by
solving the linear system
(5.8) (Xs+ λsA) zu= Xs ! UT 0 ys " = 0 ΦT2ys .
Premultiplying (5.8) with U and using (5.2) again, we see that the finite unstable eigenpairs (λu, zu) of P (λ) can be computed by
ζu,1 ζu,2 = U 0 ΦT2ys , zu,2= (Φ2+ λsΦ1)−1ζu,2, (5.9) zu,1= X1−1 ' ζu,1− X2+ λsH 1t zu,2 ( , zu= zu,1 zu,2 (5.10)
for u = 1, . . . , q. Note that Φ2+ λsΦ1is nonsingular since Φ2+ λΦ1has only unstable eigenvalues by (5.3). The vectors zu,2 in (5.9) can be found in O(q3) flops via a Hessenberg-triangular form of the pair (Φ2, Φ1) obtained by the qz-algorithm. So (5.9) requires O(q3) flops, while (5.10) requires 7mq3flops.
In the linearization approach, the computation of stable and unstable eigenvectors involves the successive application of the inverses of the two potentially ill-conditioned matrices used in the initial deflation process [5]. In our solvent approach, the matrix Xs used in the computation of stable eigenpairs is usually well conditioned. So we expect to have better accuracy in the computed results, at least for stable eigenpairs, when using the solvent approach proposed in this paper.
We now present numerical results on three sets of test data generated by a fi-nite element package, with (q, m) = (159, 11), (303, 19), (705, 51), respectively. The matrices M and K are given by (2.1) and (2.2), and we take D = 0.8M + 0.2K. All numerical experiments are carried out in MATLAB 2008b with machine precision
eps ≈ 2.22 × 10−16.
Our solvent approach is efficient since we have fully exploited the sparse struc-ture in the QEP. The only uncertainty is the number of iterations needed for the convergence of{Fk} and {Gk} in Algorithm 4.6. In Tables 5.1–5.3 we give Fk+1−
Fk 2/ Fk 2 for the three pairs of (q, m) and for ω = 100, 1000, 3000, 5000, respec-tively. The values ρ = ρ(Xs−1A) are also given for the ω values. From the tables we can see that the sequence{Fk} converges within 10 iterations for each ω. The
conver-gence behavior of{Gk} is roughly the same, as indicated by Theorem 4.1. There is no
significant difference in the performance of Algorithm 4.6 for different values of (q, m). Table 5.1
Fk+1− Fk2/Fk2 for differentω values, (q, m) = (159, 11).
ω = 100 ω = 1000 ω = 3000 ω = 5000 k ρ = 0.9622 ρ = 0.8831 ρ = 0.8080 ρ = 0.7569 1 2.0e − 02 1.4e − 02 9.7e − 03 8.1e − 03 2 3.4e − 03 1.8e − 03 1.7e − 03 1.5e − 03 3 7.4e − 04 6.0e − 04 2.9e − 04 1.5e − 04 4 3.2e − 04 8.2e − 05 9.3e − 06 1.8e − 06 5 1.0e − 04 1.6e − 06 9.3e − 09 2.4e − 10 6 8.5e − 06 5.4e − 10 9.5e − 15 2.3e − 18 7 6.1e − 08 6.2e − 17 0
8 3.2e − 12 0
9 1.1e − 22
10 0
Table 5.2
Fk+1− Fk2/Fk2 for differentω values, (q, m) = (303, 19).
ω = 100 ω = 1000 ω = 3000 ω = 5000 k ρ = 0.9307 ρ = 0.7933 ρ = 0.6692 ρ = 0.5953 1 2.2e − 02 1.3e − 02 1.1e − 02 8.4e − 03 2 3.9e − 03 1.9e − 03 9.3e − 04 5.9e − 04 3 1.2e − 03 2.4e − 04 3.8e − 05 9.5e − 06 4 2.3e − 04 6.0e − 06 6.2e − 08 2.4e − 09 5 2.3e − 05 3.6e − 09 1.6e − 13 1.4e − 16
6 2.3e − 07 1.3e − 15 0 0
7 2.4e − 11 0
8 1.3e − 20
9 0
Table 5.3
Fk+1− Fk2/Fk2 for differentω values, (q, m) = (705, 51).
ω = 100 ω = 1000 ω = 3000 ω = 5000 k ρ = 0.9593 ρ = 0.8745 ρ = 0.7925 ρ = 0.7406 1 1.1e − 01 1.0e − 01 7.0e − 02 5.7e − 02 2 2.8e − 02 1.2e − 02 1.0e − 02 8.8e − 03 3 4.7e − 03 3.6e − 03 1.5e − 03 7.8e − 04 4 2.1e − 03 4.2e − 04 3.8e − 05 6.4e − 06 5 5.7e − 04 5.8e − 06 2.2e − 08 4.3e − 10 6 4.0e − 05 1.1e − 09 7.7e − 15 2.9e − 19 7 1.9e − 07 3.5e − 17 0 8 4.6e − 12 0 9 0 10−20 10−10 100 1010 1020 10−35 10−30 10−25 10−20 10−15 | λ |
Relative residuals of eigenpairs
(q,m) = (159,11), ω = 1000
SDA_CHLW SA_HLQ SDA_GL
Fig. 5.1.
To show numerically that our method has better accuracy than existing methods, we compare our method (SDA GL) to the method in [5] (SDA CHLW) and the method SA-I in [16] (SA HLQ). The latter method has been shown in [16] to have better accuracy than two other methods compared there.
To measure the accuracy of an approximate eigenpair (λ, z) for P (λ) we use the relative residual
(5.11) RRes = λ
2ATz + λQz + Az
2
(|λ|2 A F+|λ| Q F+ A F) z 2.
In Figures 5.1–5.3 we plot for ω = 1000 the relative residuals of approximate eigen-pairs for (q, m) = (159, 11), (303, 19), (705, 51), respectively. Indeed, our new method (SDA GL) has significantly better accuracy for stable eigenpairs.
10−30 10−20 10−10 100 1010 1020 1030 10−40 10−35 10−30 10−25 10−20 10−15 | λ |
Relative residuals of eigenpairs
(q,m) = (303,19), ω = 1000 SDA_CHLW SA_HLQ SDA_GL Fig. 5.2. 10−30 10−20 10−10 100 1010 1020 1030 10−40 10−35 10−30 10−25 10−20 10−15 | λ |
Relative residuals of eigenpairs
(q,m) = (705,51), ω = 1000
SDA_CHLW SA_HLQ SDA_GL
Fig. 5.3.
6. Conclusion. We have solved a structured quadratic eigenvalue problem
effi-ciently and accurately by using a structure-preserving doubling algorithm in the sol-vent approach. The doubling algorithm has fast convergence and exploits the sparsity of the QEP. Theoretical issues involved in this solvent approach are settled satisfac-torily. In particular, we present a generalization of the classical Bendixson’s theorem to bounded linear operators in infinite-dimensional Hilbert spaces, which could also be useful elsewhere. We also mention that the solvent approach studied in this paper
can also be applied to QEPs with more general sparsity structures, such as the QEPs arising in SAW-filter simulations [29].
Acknowledgments. We thank Dr. Chin-Tien Wu from National Chiao Tung
University for discussions about the finite element model that leads to the quadratic eigenvalue problem and for generating the test data used in our numerical experiments. We are also grateful to two referees for their very helpful comments.
REFERENCES
[1] A. G. Baskakov, Dichotomy of the spectrum of non-self-adjoint operators, Sibirsk. Mat. Zh., 32 (1991), pp. 24–30 (in Russian); Siberian Math. J., 32 (1992), pp. 370–375 (in English). [2] I. Bendixson, Sur les racines d’une ´equation fondamentale, Acta Math., 25 (1902), pp. 359–365
(in French).
[3] D. A. Bini, L. Gemignani, and B. Meini, Computations with infinite Toeplitz matrices and polynomials, Linear Algebra Appl., 343–344 (2002), pp. 21–61.
[4] C.-Y. Chiang, E. K.-W. Chu, C.-H. Guo, T.-M. Huang, W.-W. Lin, and S.-F. Xu, Con-vergence analysis of the doubling algorithm for several nonlinear matrix equations in the critical case, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 227–247.
[5] E. K.-W. Chu, T.-M. Hwang, W.-W. Lin, and C.-T. Wu, Vibration of fast trains, palin-dromic eigenvalue problems and structure-preserving doubling algorithms, J. Comput. Appl. Math., 219 (2008), pp. 237–252.
[6] G. J. Davis, Numerical solution of a quadratic matrix equation, SIAM J. Sci. Stat. Comput., 2 (1981), pp. 164–175.
[7] J. C. Engwerda, A. C. M. Ran, and A. L. Rijkeboer, Necessary and sufficient conditions for the existence of a positive definite solution of the matrix equationX + A∗X−1A = Q, Linear Algebra Appl., 186 (1993), pp. 255–275.
[8] I. Gohberg, S. Goldberg, and M. A. Kaashoek, Classes of Linear Operators, Vol. II, Oper. Theory Adv. Appl. 63, Birkh¨auser, Basel, 1993.
[9] C.-H. Guo, Numerical solution of a quadratic eigenvalue problem, Linear Algebra Appl., 385 (2004), pp. 391–406.
[10] C.-H. Guo and P. Lancaster, Iterative solution of two matrix equations, Math. Comp., 68 (1999), pp. 1589–1603.
[11] C.-H. Guo and P. Lancaster, Algorithms for hyperbolic quadratic eigenvalue problems, Math. Comp., 74 (2005), pp. 1777–1791.
[12] U. Haagerup and S. Thorbjørnsen, A new application of random matrices: Ext(Cred∗ (F2))
is not a group, Ann. of Math. (2), 162 (2005), pp. 711–775.
[13] N. J. Higham and H.-M. Kim, Numerical analysis of a quadratic matrix equation, IMA J. Numer. Anal., 20 (2000), pp. 499–519.
[14] A. Hilliges, Numerische L¨osung von quadratischen eigenwertproblemen mit Anwendung in der Schienendynamik, Master’s thesis, Technical University Berlin, Berlin, Germany, 2004. [15] A. Hilliges, C. Mehl, and V. Mehrmann, On the solution of palindromic eigenvalue prob-lems, in Proceedings of the 4th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS), Jyv¨askyl¨a, Finland, 2004.
[16] T.-M. Huang, W.-W. Lin, and J. Qian, Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration of fast trains, SIAM J. Matrix Anal. Appl., 30 (2009), pp. 1566–1592.
[17] I. C. F. Ipsen, Accurate eigenvalues for fast trains, SIAM News, 37 (2004), pp. 1–2.
[18] D. Kressner, C. Schr¨oder, and D. S. Watkins, Implicit QR algorithms for palindromic and even eigenvalue problems, Numer. Algorithms, 51 (2009), pp. 209–238.
[19] P. Lancaster, U. Prells, and L. Rodman, Canonical structures for palindromic matrix polynomials, Oper. Matrices, 1 (2007), pp. 469–489.
[20] R.-C. Li, W.-W. Lin, and C.-S. Wang, Structured backward error for palindromic polynomial eigenvalue problems, Numer. Math., 116 (2010), pp. 95–122.
[21] W.-W. Lin and S.-F. Xu, Convergence analysis of structure-preserving doubling algorithms for Riccati-type matrix equations, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 26–39. [22] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann, Structured polynomial eigenvalue
problems: Good vibrations from good linearizations, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 1029–1051.
[23] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann, Numerical methods for palindromic eigenvalue problems: Computing the anti-triangular Schur form, Numer. Linear Algebra Appl., 16 (2009), pp. 63–86.
[24] C. van der Mee, G. Rodriguez, and S. Seatzu,LDU factorization results for bi-infinite and semi-infinite scalar and block Toeplitz matrices, Calcolo, 33 (1996), pp. 307–335. [25] B. Meini, Efficient computation of the extreme solutions of X + A∗X−1A = Q and X −
A∗X−1A = Q, Math. Comp., 71 (2002), pp. 1189–1204.
[26] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev., 43 (2001), pp. 235–286.
[27] J. Weidmann, Linear Operators in Hilbert Spaces, Grad. Texts in Math., 68, Springer-Verlag, New York, 1980.
[28] H. Wielandt, On eigenvalues of sums of normal matrices, Pacific J. Math., 5 (1955), pp. 633– 638.
[29] S. Zaglmayr, Eigenvalue Problems in SAW-filter Simulations, Master’s thesis, Johannes Ke-pler University, Linz, Austria, 2002.
[30] X. Zhan, Computing the extremal positive definite solutions of a matrix equation, SIAM J. Sci. Comput., 17 (1996), pp. 1167–1174.
