STRUCTURE-PRESERVING ALGORITHMS FOR PALINDROMIC QUADRATIC EIGENVALUE PROBLEMS ARISING FROM VIBRATION OF FAST TRAINS

全文

(1)c 2009 Society for Industrial and Applied Mathematics . Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. SIAM J. MATRIX ANAL. APPL. Vol. 30, No. 4, pp. 1566–1592. STRUCTURE-PRESERVING ALGORITHMS FOR PALINDROMIC QUADRATIC EIGENVALUE PROBLEMS ARISING FROM VIBRATION OF FAST TRAINS∗ TSUNG-MING HUANG† , WEN-WEI LIN‡ , AND JIANG QIAN§ Abstract. In this paper, based on Patel’s algorithm (1993), we propose a structure-preserving algorithm for solving palindromic quadratic eigenvalue problems (QEPs). We also show the relationship between the structure-preserving algorithm and the URV-based structure-preserving algorithm by Schr¨ oder (2007). For large sparse palindromic QEPs, we develop a generalized -skew-Hamiltonian implicitly restarted shift-and-invert Arnoldi algorithm for solving the resulting -skew-Hamiltonian pencils. Numerical experiments show that our proposed structure-preserving algorithms perform well on the palindromic QEP arising from a finite element model of high-speed trains and rails. Key words. palindromic quadratic eigenvalue problem, -symplectic pencil, -skew-Hamiltonian pencil AMS subject classifications. 65F15, 15A18, 15A57 DOI. 10.1137/080713550. 1. Introduction. In this paper, we consider the palindromic quadratic eigenvalue problem (QEP) of the form (1.1). P(λ)x ≡ (λ2 A 1 + λA0 + A1 )x = 0,. where λ ∈ C, x ∈ Cn \{0} and A1 , A0 ∈ Cn×n with A 0 = A0 . Note that the superscript “” denotes the complex transpose. The scalar λ and the nonzero vector x in (1.1) are the eigenvalue and the associated eigenvector of P(λ), respectively. The underlying matrix polynomial P(λ) has the property that reversing the order of the coefficients, followed by taking the transpose, leads back to the original matrix polynomial, which explains the word “palindromic.” Consequently, taking the transpose of (1.1), we easily see that the eigenvalues of P(λ) satisfy the “symplectic” property; that is, they are paired with respect to the unit circle, containing both an eigenvalue λ and its reciprocal 1/λ (with 0 and ∞ considered to be reciprocal). The palindromic QEP (1.1) was first raised in the study of the vibration in the structural analysis for fast trains in Germany [3, 4], associated with the company SFE GmbH in Berlin. Existing fast train systems, like the Japanese Shinkansen, the French TGV, and the German ICE, are being modernized and expanded. Vibration is produced from the interaction between the wheels of a train and the rails underneath. Due to the ever increasing speed (currently up to 300 km/hr) of modern trains, the study of its vibration becomes an important task. Research does not only contribute ∗ Received. by the editors January 17, 2008; accepted for publication (in revised form) by F. Tisseur September 3, 2008; published electronically January 16, 2009. http://www.siam.org/journals/simax/30-4/71355.html † Department of Mathematics, National Taiwan Normal University, Taipei, 116, Taiwan (min@ math.ntnu.edu.tw). This author’s work was partially supported by the National Science Council and the National Center for Theoretical Sciences in Taiwan. ‡ Department of Applied Mathematics, National Chiao-Tung University, Hsinchu, 300, Taiwan ([email protected]). § School of Sciences, Beijing University of Posts and Telecommunications, Beijing, 100876, China; Department of Mathematics, National Tsinghua University, Hsinchu, 300, Taiwan (jqian104@gmail. com). This author’s research was partly supported by Project NSFC 10571007. 1566. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(2) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1567. towards the increased comfort of passengers, in terms of lower noise and vibration levels. More importantly, the safety in the operation of the trains will be improved, and the operational and construction costs will be optimized [4, 5, 12, 13]. In addition, innovative designs of railway bridges, embedded rail structures, and train suspension systems require accurate resolution of the vibration. A standard approach for solving the palindromic QEP (1.1) is to transform it into a 2n × 2n linear eigenvalue problem (1.2). 0 A1. I A0. . x λx. . =λ. I 0. 0 −A 1. . x λx. . and compute its generalized Schur form (see [23]). However, the symplectic property of the eigenvalues of (1.1) is not preserved by computation, generally, producing large numerical errors ([5]). Recently, some pioneering work [4, 12, 13] proposed a good linearization which linearizes the palindromic QEP (1.1) into the form λZ +Z, which preserves symplecticity to some extent, and suggested some structure-preserving solution methods. This leads to a vast improvement over previous approaches. Later, a QR-like algorithm [19] and a Jacobi-type method [4] combined with the Laub trick, a preprocessing step of the generalized Schur form [11], have been developed for solving the palindromic linear pencil λZ + Z. However, the latter method works well, only if there are no eigenvalues near ±1. The Jacobi method typically needs about O(n3 log(n)) flops and the QR-like algorithm is of O(n4 ) flops. Recently, a URVdecomposition-based structured method of cubic complexity was developed in [20] to solve the palindromic linear pencil λZ + Z, producing eigenvalues which are paired to working precision. In section 3, we will show that the URV-based method [20] is mathematically equivalent to applying the structure-preserving algorithm in section 2 to the enlarged 2n × 2n palindromic quadratic pencil ζ 2 Z + ζ0 + Z (with ζ 2 = λ). On the other hand, a structure-preserving doubling algorithm was developed in [1] via the computation of a solvent of a nonlinear matrix equation associated with (1.1). The numerical results show much promise but the convergence theory holds only when the algorithm does not break down. As mentioned before, the linearization (1.2) generally cannot preserve the symplectic structure. Fortunately, the special linearization for (1.1) (see [1] or [10]) (1.3). (M − λL)z ≡. obtained by setting y = satisfies (1.4). 1 λ A1 x. A1 −A0. 0 0 −λ −I A1. I 0. x =0 y. and multiplying the second equation of (1.3) by λ MJ M = LJ L ,. where J ≡ J2n is the 2n × 2n matrix [ −I0 n I0n ]. In other words, the pencil M − λL or the matrix pair (M, L) in (1.3) preserves the symplectic structure of (1.4) and is said to be -symplectic. For a real matrix pair (M, L) satisfying (1.4), a structure-preserving (S + S −1 )transform for the computation of all its eigenvalues is proposed by [9] and a numerically stable algorithm for reducing the transformed pair to a block triangular condensed form by using only orthogonal transformations was developed by Patel [16]. It is perfectly suitable for the -symplectic pair, but not applicable to the complex conjugate. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(3) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1568. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. symplectic pair (i.e., MJ MH = LJ LH ). In this paper, we adapt Patel’s approach to solve the -symplectic pencil in (1.3) resulting from the palindromic QEP (1.1). Only unitary transformations are used and the symplectic structure is fully preserved, which make the method attractive. It is worth mentioning that the (S + S −1 )-transform is, in general, a nonlinear transform as in solving the discrete-time optimal control problem [9, 16]. However, the special form in (1.3) leads to a linear (S + S −1 )-transform without involving any matrix multiplication. In some applications, the matrices A1 and A0 in (1.1) (and hence M and L in (1.3)) can be large and sparse and only the eigenvalues in a specified region are required. To accomplish this, the shift-and-invert (implicitly restarted) Arnoldi algorithm [7, 17, 21] is one of the most widely used standard techniques for computing selected eigenvalues of the large sparse matrix pencil M − λL. In this approach, the corresponding shifted and inverted matrix is reduced to a Hessenberg form which no longer has the desirable symplectic structure. Mehrmann and Watkins [15] developed a structure-preserving skew-Hamiltonian, isotropic, implicitly restarted shift-and-invert Arnoldi algorithm (SHIRA) for the computation of eigenpairs of a large sparse real skew-Hamiltonian/Hamiltonian pencil by transforming the pencils to a skew-Hamiltonian operator. In fact, SHIRA can be straightforwardly extended to solve a skew-Hamiltonian/Hamiltonian pencil in the complex transpose case (not in the complex conjugate case), referred to as SHIRA. We first transform the -symplectic pencil to a -skew-Hamiltonian eigenvalue problem by using the (S + S −1 )-transform, then SHIRA is applied to the resulting -skew-Hamiltonian matrix. On the other hand, to avoid explicitly forming the -skew-Hamiltonian matrix in the above transformation, we also develop a generalized -skew-Hamiltonian implicitly restarted shift-and-invert Arnoldi algorithm (GSHIRA) for solving the -skew-Hamiltonian pencil resulting from the (S + S −1 )transform of the symplectic pencil M − λL. We introduce some definitions that will be used frequently in this paper. Definition 1.1. (i) A matrix A ∈ Cn×n is called -symmetric or -skew-symmetric if it satisfies A = A or A = −A, respectively. (ii) A matrix U ∈ C2n×2n is called -symplectic if U J U = J ; a pencil M − λL ∈ C2n×2n or the matrix pair (M, L) is called -symplectic if MJ M = LJ L . (iii) A matrix H ∈ C2n×2n is called -Hamiltonian or -skew-Hamiltonian if it satisfies (HJ ) = HJ or (HJ ) = −HJ , respectively. (iv) A pencil K − λN ∈ C2n×2n or the matrix pair (K, N ) is called -skewHamiltonian if K and N are -skew-Hamiltonian. (v) Let X, Y ∈ C2n×m (1 ≤ m ≤ n); X is called -isotropic if X J X = 0m ; and X and Y are called -bi-isotropic if X J Y = 0m . Throughout this paper, A and AH denote the transpose and conjugate transpose of a matrix A, respectively. We denote the m × n zero matrix by 0m,n , and the zero and identity matrices of order n by 0n and In , respectively. The ith column of In is denoted by ei . 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, :) and A(:, k : ) select the rows i to j and the columns k to , respectively, of A. The paper is organized as follows. In section 2, we briefly present the structurepreserving algorithm based on Patel’s method [16] for solving palindromic QEPs. In. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(4) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1569. section 3, we show the relationship between the structure-preserving algorithm and the URV-based structured method proposed by Schr¨ oder [20]. In section 4, based on the SHIRA developed in [15], we introduce the -skew-Hamiltonian implicitlyrestarted shift-and-invert Arnoldi algorithm (SHIRA) for solving the resulting skew-Hamiltonian matrix. In section 5, a generalized -skew-Hamiltonian implicitlyrestarted shift-and-invert Arnoldi algorithm (GSHIRA) for solving the resulting -skew-Hamiltonian pencils is developed. We present some numerical results of the proposed algorithms, using examples from a finite element model of fast trains [1], in section 6. Conclusions are given in section 7. 2. Structure-preserving algorithm I. We adapt Patel’s algorithm [16] applying to the (S + S −1 )-transform of a -symplectic matrix pair for the computation of all its eigenpairs. Let (M, L) be a -symplectic pair. The (S + S −1 )-transform (Ms , Ls ) of (M, L) is defined by (see [9]) Ms ≡ MJ L + LJ M ,. (2.1). Ls ≡ LJ L .. We first give the relationship between eigenpairs of a -symplectic pencil and its (S + S −1 )-transform. Theorem 2.1. Let (M, L) be a -symplectic pair and (Ms , Ls ) be its (S + S −1 )transform. Then (i) μ is a double eigenvalue of (Ms , Ls ) if and only if ν, ν1 are eigenvalues of (M, L), where ν, ν1 are two roots of the quadratic equation λ + λ1 = μ. (ii) Let x and y be linearly independent eigenvectors of (L , M ) corresponding to ν and ν1 , respectively, i.e., (L −νM )x = 0 and (L − ν1 M )y = 0. Then x and y are two linearly independent eigenvectors of (Ms , Ls ) corresponding to μ = ν + ν1 . (iii) Furthermore, from (ii), if zs = αx + βy (with αβ = 0) is an eigenvector of (Ms , Ls ) corresponding to μ = ν + ν1 (μ = ±2), i.e., (Ms − μLs )zs = 0, then J (L − ν1 M )zs and J (L − νM )zs are the eigenvectors of (M, L) corresponding to ν and ν1 , respectively. Proof. (i) As in [9], since MJ M = LJ L , by (2.1) it holds that 1 Ms − μLs = MJ L + LJ M − ν + LJ L ν 1 = (M − νL)J L − M ν 1 = M − L J (L − νM ). (2.2) ν Hence (i) follows. (ii) From the last two equations of (2.2), it follows that 1 (Ms − μLs ) x = M − L J (L − νM ) x = 0, ν and. 1 (Ms − μLs ) y = (M − νL)J L − M y = 0. ν. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(5) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1570. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. (iii) By applying the last two equations of (2.2) again, it remains to show only that J (L − ν1 M )zs = 0 and J (L − νM )zs = 0. From (ii) we have 1 1 1 J L − M zs = J L − M (αx + βy) = αJ L − M x = 0. ν ν ν Similarly, J (L − νM )zs = J (L − νM )(αx + βy) = βJ (L − νM )y = 0. Theorem 2.2. Let (M, L) be the -symplectic pair as in (1.3) and (Ms , Ls ) be its (S + S −1 )-transform. If zs = [z1 , z2 ] with z1 , z2 ∈ Cn is an eigenvector of (Ms , Ls ) corresponding to μ = ν + ν1 (μ = ±2), then z1 + ν1 z2 and z1 + νz2 are eigenvectors of P(λ) in (1.1) corresponding to ν and ν1 , respectively. Proof. From (iii) of Theorem 2.1 we compute (2.3) 1 1 J L − M zs (1 : n) = z1 + z2 , (J (L − νM )zs )(1 : n) = z1 + νz2 . ν ν Then, from (1.3) and (2.3), it follows that P(ν)(z1 + ν1 z2 ) = 0 and P( ν1 )(z1 + νz2 ) = 0. Note that from (1.3), we have (Ms , Ls ) = (MJ L + LJ M , LJ L ) 0 −A1 A1 − A A0 1 = , −A0 A1 − A A1 0 1 0 −A1 A0 A1 − A1 = , J 0 −A A1 − A A0 1 1 (2.4). ≡ (K, N )J .. From (2.4), if z is an eigenvector of (K, N ) corresponding to μ, then zs = J z is the eigenvector of (Ms , Ls ) corresponding to the same μ. Remark 2.1. (i) The (S + S −1 )-transform (Ms , Ls ) in (2.1) of a -symplectic pair, in general, is a nonlinear (quadratic) transformation. For instance, the (S + S −1 )A 0 I G transform of the symplectic pair of the form (M, L) ≡ ([ −H I ], [ 0 A ]) with H = H and G = G arisen from discrete-time optimal control problems produces a quadratic (S + S −1 )-transform which involves matrix multiplications and is not backward stable. However, the special form of the -symplectic pair (M, L) in (1.3) leads to a linear (S + S −1 )-transform as in (2.4) and does not involve any matrix multiplication. (ii) The eigenvectors of P(λ) corresponding to ν and 1/ν can be obtained from the eigenvectors of (K, N ) directly (see Theorem 2.2), not requiring us to solve any linear system or perform any matrix-vector multiplications. It is easily seen that K and N in (2.4) are both -skew-Hamiltonian. Patel [16] introduced two types of transformations that preserve the skew-Hamiltonian structure. The first type involves similarity transformations on K and N , respectively, using Given rotations G0 (i, c, s¯) := G(i, n + i, c, s¯). The second type involves equivalence transformations on K and N , respectively, by the left transformation ⊕ V ) and the right transformation Z0 := (V ⊕ U ), where the unitary Q 0 := (U. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(6) 1571. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. U, V ∈ Cn×n represent the application of Givens rotations. One can easily verify that the new transforming K and N are still -skew-Hamiltonian. Based on Patel’s approach [16] with these two types of transformations, we may reduce (K, N ) to a block triangular structure; that is, (2.5). . K11 K := Q KZ = 0 . K12 , K11. . N11 N := Q N Z = 0 . N12 , N11. where K11 ∈ Cn×n is upper Hessenberg, N11 ∈ Cn×n is upper triangular, and Q, Z are unitary satisfying Q = J ZJ .. (2.6). From (2.5), we see that the pair (K11 , N11 ) contains half of the eigenvalues of (K, N ). We then apply the QZ algorithm to (K11 , N11 ) for computing all eigenpairs {(μi , yi )}ni=1 . Consequently, {(μi , Z[ y0i ])}ni=1 are n eigenpairs of (K, N ). From (2.4), {(μi , zi (≡ J Z[ y0i ]))}ni=1 are eigenpairs of (Ms , Ls ). Finally, we compute all eigenvalues and the associated eigenvectors of P(λ) by Theorem 2.2. Algorithm 2.1 (structure-preserving algorithm I (SA I)). Input: A palindromic quadratic pencil P(λ) ≡ λ2 A 1 + λA0 + A1 with A0 , A1 ∈ Cn×n and A = A . 0 0 Output: All eigenvalues and eigenvectors of P(λ). Step 1. Form the pair (K, N ) as in (2.4); Step 2. Reduce (K, N ) to block upper triangular forms in (2.5) using unitary transformations. (See a pseudocode in Appendix A.1.); Step 3. Compute eigenpairs {(μi , yi )}n i=1 of (K11 , N11 ) defined in (2.5) by using the QZ algorithm; Step 4. Compute zi = J Z[ y0i ] ≡ [ zzi1 ], i = 1, 2, . . . , n; i2 Step 5. Compute eigenvalues νi and ν1i of P(λ) by solving ν 2 − μi ν + 1 = 0; Compute eigenvectors xi1 ≡ zi1 + ν1i zi2 , xi2 ≡ zi1 + νi zi2 corresponding to νi , ν1i , respectively, for i = 1, 2, . . . , n.. Remark 2.2. The SA I requires approximately 27n3 flops for the eigenvalues, and an additional 23n3 flops for the eigenvectors. While the QZ algorithm is applied to (M, L) directly, it requires approximately 120n3 flops for the eigenvalues and an 3 additional 260 3 n flops for the eigenvectors. Here and hereafter a flop is a floating point multiplication and addition for complex numbers, which involves 6 real flops. 3. Structure-preserving algorithm II vs. URV-based method. Recently in [4, 12, 13], a “good” linearization of the palindromic quadratic pencil (1.1) was proposed: (3.1). . λZ + Z ≡ λ. . A 1 A 1. A0 − A1 A 1. . +. A1 A0 − A 1. A1 A1. .. This preserves the “symplecticity” of the eigenvalues. In order to solve the palindromic linear eigenvalue problem of (3.1), we rewrite it into a new palindromic quadratic pencil (3.2). Q(ζ) ≡ ζ 2 Z + ζ02n + Z. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(7) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1572. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. with ζ 2 = λ. We then apply the SA I algorithm proposed in section 2 to solve the palindromic QEP of (3.2). As in (2.4), we form 0 Z − Z −Z 0 K= , N = (3.3) . Z − Z 0 0 −Z Va J4n such that By (2.5) there are unitary Ua , Va ∈ C4n×4n with Ua = J4n a a a a N11 K12 N12 a = K11 N V = (3.4) , U , Ua KV a a a a 0 (K11 ) 0 (N11 ) a where K11 ∈ C2n×2n is upper Hessenberg with {0, 2, . . . , 2n − 2}-diagonals being a zeros, N11 ∈ C2n×2n is upper triangular with {1, 3, . . . , 2n − 1}-diagonals being zeros, a a and K12 and N12 ∈ C2n×2n are skew symmetric with {1, −1, . . . , 2n − 1, −(2n − 1)}diagonals and with {2, −2, . . . , 2n − 2, −(2n − 2)}-diagonals, respectively, being zeros. Here the -diagonal of a matrix A ≡ [aij ]ni,j=1 consists of the entries {aij } with a a a a j − i = . Note that the extra zeros in K11 , N11 , K12 , and N12 are obtained by performing some suitable permutations on the special forms of (3.3) without any calculation. (See Appendix A.2 for details.) Denote. P2n = [e1 , en+1 , e2 , en+2 , . . . , en , e2n ].. (3.5) Let U. (3.6) Then we have. (3.7a). (3.7b). . =. P2n 0. 0 P2n. . . Ua ,. ⎡. 0 ⎢ R =⎢ 2 U KV ⎣ 0 0 ⎡ R3 ⎢ V = ⎢ 0 U N ⎣ 0 0. V = Va. R1 0 0 0. T1 0 0 R1. 0 R4 0 0. 0 T3 R3 0. P2n 0. 0 P2n. .. ⎤ 0 −T2 ⎥ ⎥, R2 ⎦ 0 ⎤ −T3 0 ⎥ ⎥, 0 ⎦ R4. where R1 ∈ Cn×n is upper Hessenberg, R2 , R3 , R4 ∈ Cn×n are upper triangular, T1 , T2 ∈ Cn×n are skew symmetric, and T3 ∈ Cn×n . From (3.7), we see that in order N ) it suffices to compute those to compute the eigenvalues and the eigenvectors of (K, of the matrix pair (3.8). (R1 R4−1 R2 , R3 ).. We apply the periodic QZ algorithm [2, 18] to the matrix pair in (3.8) without forming the product explicitly, which gives the n eigenpairs {(γi , yi )}ni=1 , where yi ∈ Cn . Let √ μi = γi (one branch of the square root of γi ), ηi := μi R1−1 R3 yi , and yi = [yi , ηi ] . N ). Write zi = [ zi1 , zi2 ] It follows that {(μi , zi (≡ V[ y0i ]))}ni=1 are n eigenpairs of (K, 1 2 2 and solve νi and νi for ν + (2 − μi )ν + 1 = 0. By Theorem 2.2 and (3.1), we compute the eigenvectors (3.9a). xi1 = x i1 (1 : n) + x i1 (n + 1 : 2n),. xi2 = x i2 (1 : n) + x i2 (n + 1 : 2n). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(8) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. of P(λ) corresponding to νi and (3.9b). 1 νi ,. 1573. respectively, where. 1 x i1 := zi2 − √ zi1 , νi. x i2 := zi2 −. √ νi zi1 .. Algorithm 3.1 (structure-preserving algorithm II (SA II)). Input: A palindromic quadratic pencil P(λ) ≡ λ2 A 1 + λA0 + A1 with A0 , A1 ∈ Cn×n and A = A . 0 0 Output: All eigenvalues and eigenvectors of P(λ). N ) as in (3.3); Step 1. Form the pair (K, Step 2. Reduce (K, N ) to block upper triangular forms as in (3.7) using unitary transformations of (3.4)–(3.6); Step 3. Compute eigenpairs {(γi , yi )}ni=1 of (R1 R4−1 R2 , R3 ) in (3.8) by the periodic QZ algorithm [18]; I Step 4. Compute zi = V[ y0i ] ≡ [ zzi1 ], where yi = [ √γi Rn−1 R3 ]yi i2 1 for i = 1, 2, . . . , n; Step 5. Compute νi and ν1i by solving ν 2 + (2 − γi )ν + 1 = 0; Compute eigenvectors xi1 and xi2 of P(λ) as in (3.9a) corresponding to νi , ν1i , respectively, for i = 1, 2, . . . , n. Remark 3.1. (i) In Step 3, since R1 , R4 , R2 , and R3 are already in Hessenberg-triangular form, the first step in the periodic QZ algorithm is not needed. (ii) The SA II requires 62n3 flops for the eigenvalues, and an additional 23n3 flops for the eigenvectors. Recently a URV-decomposition-based structured method was proposed in [20] for solving the palindromic linear pencil (3.1). From [20] there are unitary U, V ∈ C2n×2n such that. Πn Πn 0 R 0 −R 4 2 (3.10a) U ZV = 3 Πn T3 Πn , V (Z − Z )V = Πn R 2 Πn T2 Πn Πn R and. (3.10b). U (Z − Z)U =. 0 1 Πn R. 1 Πn −R Πn T1 Πn. ,. 1 ∈ Cn×n is upper Hessenberg, R 2 , R 3 , R 4 ∈ Cn×n are where Πn = [en , . . . , e1 ], R n×n n×n upper triangular, T1 , T2 ∈ C are skew symmetric, and T3 ∈ C . Define ⎡ ⎤ 0 0 Πn 0 ⎢ ⎥ U 0 0 0 Π 0 n ⎥ ⎢ U0 := ⎣ (3.11) U0 J4n . , V0 := J4n 0 V 0 0 In 0 ⎦ −In 0 0 0 0 Then it is easily seen that U0H KV “\hat” being over all submatrices. P2n (3.12) Ub = 0. V0 have the same forms as in (3.7) with and U0H N Furthermore, if we define 0 Ub J4n , U0 , Vb := J4n P2n. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(9) 1574. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. then we have (3.13). b= UbH KV. . b K11 0. b K12 b (K11 ). ,. Vb = UbH N. . b N11 0. b N12 b (N11 ). ,. b b b b where K11 , K12 , N11 , and N12 are of the same forms as in (3.4). a b a b are unreduced, and N11 and N11 are nonsingular Theorem 3.1. If K11 and K11 (see (3.4) and (3.13)), then the SA II is mathematically equivalent to the URV-based structured method. Va J4n , Proof. Denote Va := [V1a , V2a ] with Via ∈ C4n×2n (i = 1, 2). Since Ua = J4n a a it holds that Ua = [J4n V2 , −J4n V1 ]. From (3.4), it follows that. a = J4n V a K a , KV 1 2 11. (3.14). V a = J4n V a N a . N 1 2 11. This implies that a=N V a (N a )−1 K a . KV 1 1 11 11 Since the first columns of V1a and V1b (Vb ≡ V1b , V2b ) are both e1 , by applying the implicit Q-theorem to (3.15), the matrices Ua and Va are uniquely determined, and Ua = Ub and Va = Vb . (3.15). 4. -skew-Hamiltonian Arnoldi method. Based on SHIRA [15], in this section we briefly introduce the structure-preserving -skew-Hamiltonian Arnoldi algorithm to compute the desired eigenpairs of a -skew-Hamiltonian B. As in (2.4), using the (S + S −1 )-transform, we transform M − λL of (1.3) into a -skew-Hamiltonian pencil K − μN by (4.1) K − μN ≡ (LJ M + MJ L ) − μLJ L J . / σ(M, L). Next, we derive the shift-invert transformation of K − μN . Let λ0 ∈ Then, from Theorem 2.2(i), we have μ0 ≡ λ0 + λ10 ∈ / σ(K, N ). Define the shift-invert −μ for K − μN with μ transformation K N = 1 and . μ−μ0. A 1 0. 0 A1. . (4.2a). ≡ −λ0 N = −λ0 LJ L J = λ0 K. (4.2b). ≡ −λ0 (K − μ0 N ) = −λ0 (LJ M + MJ L − μ0 LJ L )J . N. Substituting μ0 = λ0 +. (4.3). 1 λ0. ,. can be factorized as into (4.2b), N. 1 N = −λ0 LJ M + MJ L − λ0 + LJ L J λ0 = (M − λ0 L) J M − λ0 L J ≡ N1 N2 ,. where (4.4). N1 = M − λ0 L,. N2 = J (M − λ0 L )J . = are nonsingular and satisfy N2 J = J N1 . The generalized eigenvalue problem Kz μ N z is then equivalent to the eigenvalue problem By = μ y, where y = N2 z and (4.5). −1 . B ≡ N1−1 KN 2. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(10) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1575. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. = JK and N J = J N1 , we find that B satisfies Using the facts that KJ 2 N − = N −1 J K N − = N −1 KN −1 J = BJ , N − = N −1 KJ J B = J N2− K 1 1 1 1 1 1 2 and hence B is again -skew-Hamiltonian. We now define the Krylov matrix with respect to u1 and j (1 ≤ j ≤ n) by (4.6). Kj ≡ Kj [B, u1 ] = [u1 , Bu1 , . . . , B j−1 u1 ]. and state two useful theorems from [15]. Note that these theorems are slightly different from the originals, but the proofs are almost identical to the ones in [15]. Theorem 4.1 (see [15]). Let B ∈ C2n×2n be -skew-Hamiltonian and Kj ≡ Kj [B, u1 ] (1 ≤ j ≤ n) be a Krylov matrix with rank(Kj ) = j. Then span(Kj ) is j is a QR-factorization, then -isotropic and if Kj = Uj R (4.7). j + u BUj = Uj H j+1 e j ,. j ∈ Cj×j is unreduced upper Hessenberg, Uj ∈ C2n×j is orthonormal and where H -isotropic, and u j+1 ∈ C2n is a suitable vector such that (4.8). UjH u j+1 = 0. Uj J u j+1 = 0.. and. Theorem 4.2 (see [15]). Let B ∈ C2n×2n be -skew-Hamiltonian. If rank (Kn [B, u1 ]) = n, then there is a unitary -symplectic matrix U with Ue1 = u1 such that. (4.9). H. U BU =. n H 0. . n N H. ,. n. n is unreduced upper Hessenberg and N n is -skew-symmetric. where H Based on Theorem 4.2, the jth step of the Arnoldi process is given by. (4.10). hj+1,j uj+1 = Buj −. j . hij ui ,. i=1. where hij = uH i Buj , i = 1, . . . , j, and hj+1,j > 0 is chosen so that uj+1 2 = 1. In order to ensure that the space span{u1 , . . . , uj+1 } is -isotropic to working precision, the jth step of the -isotropic Arnoldi process is modified by. (4.11). hj+1,j uj+1 = Buj −. j i=1. hij ui −. j . ¯i , tij J u. i=1. where hij = uH i Buj , tij = −ui J Buj , i = 1, . . . , j, and hj+1,j > 0 is chosen so that. uj+1 2 = 1. We present the SHIRA-method.. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(11) 1576. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. Algorithm 4.1 (SHIRA). Input: -skew-Hamiltonian matrix B with starting vector u1 . with BU = U H , U H U = I Output: U and upper Hessenberg matrix H and U J U = 0. Use (4.11) with starting vector u1 to generate the th step of the -isotropic Arnoldi factorization: + BU = U H h+1, u+1 e . For k = 1, 2, . . . , until wanted eigenpairs of B are convergent, Use (4.11) to extend the th step of the -isotropic Arnoldi factorization to the ( + p)th step of the -isotropic Arnoldi factorization: +p + BU+p = U+p H h+p+1,+p u+p+1 e +p . Use standard implicitly restarted step for the Arnoldi algorithm [8] to reform a new th step of the -isotropic Arnoldi factorization. End. Remark 4.1. (i) h+1, is set to zero if | h+1, | < tol(| h, | + | h+1,+1 |) for some stopping tolerance “tol.” , i.e., H vi = θi vi . Let yi = U vi be the Ritz (ii) Let (θi , vi ) be an eigenpair of H vector of B corresponding to the Ritz value θi . Then from (4.7) and (4.8), we have. Byi − θi yi 2 = BUvi − θi U vi 2 + u = (U H +1, e )vi − θi U vi 2 = U (H vi − θi vi ) + h+1, (e vi )u+1 2 . = | h+1, ||e vi |.. 5. Generalized -skew-Hamiltonian Arnoldi method. We now consider =μ z, where K and N are -skew-Hamiltonthe generalized eigenvalue problem Kz N can be reduced to ian given in (4.2). Based on the reduction method [16], K − μ N block triangular condensed forms K11 K12 N11 N12 V (K − μ N )U = (5.1) −μ , 0 K11 0 N11 where K11 , N11 ∈ Cn×n are, respectively, upper Hessenberg and upper triangular, and V and U ∈ C2n×2n are unitary satisfying (5.2). V = J UJ .. In order to solve a large sparse product or a periodic eigenvalue problem, recently, a product (or a periodic) Arnoldi process and a product Krylov process were, respectively, proposed by Kressner’s book [6, section 4.2.5] and Watkins’ book [24, section 9.10]. Using the result of Theorem 4.1, we adopt the idea of the periodic Arnoldi process [6, section 4.2.5] to develop a generalized -skew-Hamiltonian algorithm which preserves the structure of (5.1) for the computation of the desired =μ z. eigenpairs of Kz N −1 be -skew-Hamiltonian defined in (4.5). Let Theorem 5.1. Let B ≡ N1−1 KN 2 = N1 N2 and Kj ≡ Kj [B, u1 ] be the Krylov matrix with rank(Kj ) = j. If N (5.3). N2−1 Kj = Zj R2,j. and. N1 Kj = Yj R1,j. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(12) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1577. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. are QR-factorizations, where Zj , Yj ∈ C2n×j are orthonormal and R2,j , R1,j are nonsingular upper triangular, then we have j = Yj Hj + yj+1 e KZ j. (5.4) and. Zj = Yj Rj , N. (5.5). where Hj ∈ Cj×j is unreduced upper Hessenberg, Rj ∈ Cj×j is nonsingular upper triangular, and Yj and Zj are -bi-isotropic such that YjH yj+1 = 0. (5.6). and. Zj J yj+1 = 0. for a suitable yj+1 ∈ C2n . j be the QR-factorization of Kj with R j being nonsingular Proof. Let Kj = Uj R upper triangular. From Theorem 4.1, it follows that −1 Uj = Uj H j + u j+1 e N1−1 KN j . 2. (5.7). Substituting (5.3) into (5.7) we obtain j R−1 j = KN −1 Kj R−1 = KN −1 Uj R KZ 2 2 2,j 2,j j + N1 u −1 = (N1 Uj H j+1 e j )Rj R2,j H −1 H j R j R−1 ) + γj Yj YjH N1 u = Yj (R1,j R j+1 e j+1 e j + γj (I − Yj Yj )N1 u j j 2,j. (5.8). = Yj Hj + yj+1 e j ,. −1 where γj = e j Rj R2,j ej , (5.9). −1 H j R j R−1 + γj YjH N1 u j+1 e Hj = R1,j R j , j 2,j. and (5.10). j+1 . yj+1 = γj (I − Yj YjH ) N1 u. j , R1,j , and R2,j are nonsingular upper triangular, and H j is unreduced upper Since R Hessenberg, from (5.9) it follows that Hj is unreduced upper Hessenberg. Clearly, it holds that YjH yj+1 = 0 by (5.10). On the other hand, from (5.3), we also have Zj = N1 N2 Zj = N1 Kj R−1 = Yj R1,j R−1 ≡ Yj Rj , N 2,j 2,j −1 where Rj = R1,j R2,j is nonsingular and upper triangular. We now show that Yj and Zj are -bi-isotropic. By the fact that N2 J = J N1 and (5.3), it holds that. (5.11). − −1 − −1 Yj J Zj = R1,j Kj (N1 J N2−1 )Kj R2,j = R1,j Kj J Kj R2,j = 0.. From (5.8) and (5.10), we have Zj J yj+1 e j = Zj J (KZj − Yj Hj ) = Zj J KZj ,. which is -skew-symmetric. This implies that Zj J yj+1 = 0.. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(13) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1578. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. −1 be -skew-Hamiltonian defined in (4.5) Theorem 5.2. Let B = N1−1 KN 2 = N1 N2 . If rank(Kn [B, u1 ]) = n, then there are unitary matrices U and V and N satisfying (5.2) and Ve1 = N1 u1 / N1 u1 2 such that Hn S n Rn Tn V KU = (5.12) , V NU = , 0 Hn 0 Rn where Hn is unreduced upper Hessenberg, Rn is nonsingular upper triangular, and Sn and Tn are -skew-symmetric. Proof. Applying Theorem 5.1 for j = n, we have yn+1 being orthogonal to Yn and J Z¯n . This implies that yn+1 = 0. Then (5.4) and (5.5) become (5.13). n = Yn Hn KZ. Zn = Yn Rn , N. and. where Hn is unreduced upper Hessenberg and Rn is nonsingular upper triangular. Let U ≡ Zn −J Y¯n , V ≡ Yn −J Z¯n . Clearly, (5.14). ZnH Zn = In ,. YnH Yn = In ,. and Yn J Zn = 0n .. and N J are -skew symmetric, from (5.13)– Then U and V satisfy (5.2). Since KJ (5.14), (5.12) follows. Based on Theorem 5.2, we now introduce a generalized -isotropic Arnoldi process which produces -bi-isotropic matrices Zj and Yj+1 at the jth step. By the recursive definition of j, let us first assume that the -bi-isotropic matrices Zj−1 and Yj satisfy (5.4) and (5.5) with j := j − 1. That is, the (j − 1)th step of the generalized -isotropic Arnoldi process generates (5.15). Zj−1 = Yj−1 Rj−1 . N. Now, we compare the jth columns of both sides in (5.5) which give (5.16). zj = N. j−1 . rij yi + rjj yj .. i=1. With (5.15), (5.16) it can be rewritten as (5.17). −1 −1 yj − rjj zj = N. j−1 . rij zi ,. i=1. where (5.18). −1 −1 [ r1j , . . . , rj−1,j ] := −rjj Rj−1 [r1j , . . . , rj−1,j ] .. Since ZjH Zj = Ij , the coefficient rij in (5.17) can be evaluated by (5.19). −1 yj , rij = zjH N. i = 1, . . . , j − 1,. r1j , . . . , rj−1,j ] of (5.19) and rjj in (5.17) is chosen so that zj 2 = 1. Substituting [ into (5.18), we obtain the coefficient vector [r1j , . . . , rj−1,j ] .. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(14) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1579. In exact arithmetic, zj is orthogonal to J Y¯j automatically. As before, roundoff errors cause zj J yi , i = 1, . . . , j, to be tiny values. Thus, the jth step of the generalized -isotropic Arnoldi process for zj should be modified by (5.20a). −1 rjj zj. −1 yj − =N. where (5.20b). j−1 i=1. sij =. yi J . rij zi −. j . −1 yj − N. j−1 . sij J y¯i ,. i=1. rij zi. ,. i = 1, . . . , j.. i=1. From (5.4), similar to (4.11), the jth step of the generalized -isotropic Arnoldi process for yj+1 is given by (5.21a). j− hj+1,j yj+1 = Kz. j . hij yi −. i=1. j . tij J z¯i ,. i=1. where (5.21b). j, hij = yiH Kz. j, tij = zi J Kz. i = 1, . . . , j,. and hj+1,j > 0 is chosen so that yj+1 2 = 1. Combing (5.20) and (5.21), we state the jth step of the generalized -isotropic Arnoldi process. Algorithm 5.1 (the jth generalized -isotropic Arnoldi step). and N , upper triangular R(1 : j − 1, 1 : j − 1), Input: -skew-Hamiltonian K Yj = [y1 , . . . , yj ] and Zj−1 = [z1 , . . . , zj−1 ] with YjH Yj = Ij , H Zj−1 Zj−1 = Ij−1 , and Yj J Zj−1 = 0. Output: [h1,j , . . . , hj+1,j ], R(1 : j, 1 : j), yj+1 , and zj . Compute zj in (5.20) by using the modified Gram–Schmidt step: zj = yj ; Solve N For i = 1, . . . , j − 1 rij = ziH zj , zj = zj − rij zi End Set R(j, j) := zj −1 2 , zj := R(j, j)zj , and R(1 : j − 1, j) := −R(j, j)R(1 : j − 1, 1 : j − 1)[ r1j , · · · , rj−1,j ] ; ¯ Reorthogonalize zj to J Yj : For i = 1, . . . , j sij = yi J zj , zj = zj − sij J y¯i End Compute yj+1 in (5.21): j; Compute yj+1 = Kz For i = 1, . . . , j hij = yiH yj+1 , yj+1 = yj+1 − hij yi End Set hj+1,j := yj+1 2 and yj+1 := yj+1 /hj+1,j yj+1 ; For i = 1, . . . , j tij = zi J yj+1 , yj+1 = yj+1 − tij J z¯i End. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(15) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1580. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. 5.1. Implicitly restart. We now derive the implicitly restarted step for the (+ p)th step of the generalized -isotropic Arnoldi process. Suppose we have computed the ( + p)th step of the generalized -isotropic Arnoldi factorization: (5.22). +p = Y+p H+p + h+p+1,+p y+p+1 e KZ +p ,. (5.23). Z+p = Y+p R+p . N. Let {λ1 , . . . , λ , λ+1 , . . . , λ+p } be the eigenvalues of the matrix pair (H+p , R+p ), where {λ1 , . . . , λ } are the wanted eigenvalues. Let Qk and Vk for k = 1, . . . , p be unitary matrices computed by the implicit-QZ step [22, p. 147] for (H+p , R+p ) with the single shift λ+k . +p := QH · · · QH H+p V1 · · · Vp , R +p := QH · · · QH R+p V1 · · · Vp , Y+p := Let H p 1 p 1 +p := Z+p V1 · · · Vp . Then H +p and R +p are upper Hessenberg Y+p Q1 · · · Qp , and Z +p satisfy Y J Z +p = 0 because and upper triangular, respectively, and Y+p and Z +p of Y+p J Z+p = 0. Multiplying (5.22) and (5.23) by V1 · · · Vp , we get (5.24). +p + h+p+1,+p y+p+1 e V1 · · · Vp , Z +p = Y+p H K +p. (5.25). Z +p = Y+p R +p . N. Since e +p V1 = α+p e+p−1 + β+p e+p ,. by induction, the first − 1 entries of e +p V1 · · · Vp are zero. Hence a new th step of the generalized -isotropic Arnoldi factorization can be obtained by equating the first columns of (5.24) and (5.25): Z = Y H + K h+p+1,+p y+p+1 e , N Z = Y R . We summarize the above processes in Algorithm 5.2. Algorithm 5.2 (generalized implicitly restarted step). Input: given (Y+p , y+p+1 , Z+p , H+p , h+p+1,+p , R+p ); Output: (Y , y+1 , Z , H , h+1, , R ) formed a new th step of the generalized -isotropic Arnoldi factorization. The best eigenvalues are locked in (H , R ). Sort the eigenvalues of (H+p , R+p ) from best to worst according to the sorting criterion and take {λ+1 , . . . , λ+p } to be the p worst eigenvalues. Set v := h+p+1,+p e+p ; For k = 1, . . . , p, Compute unitary matrices Qk and Vk by the implicit-QZ step for (H+p , R+p ) with the single shift λ+k so that QH k H+p Vk and QH R V are upper Hessenberg and upper triangular, respectively; k +p k Update Y+p := Y+p Qk , Z+p := Z+p Vk , H+p := QH H +p Vk , k H R+p := QH R V , v := Z v; k +p k k End Set H = H+p (1 : , 1 : ), h+1, := e v, R = R+p (1 : , 1 : ), Y := Y+p (:, 1 : ), y+1 := y+p+1 , Z := Z+p (:, 1 : ).. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(16) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1581. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. We now present the GSHIRA. Algorithm 5.3 (GSHIRA). and N with starting vector y1 . Input: -skew-Hamiltonian matrices K Output: Z , Y , upper Hessenberg H , and upper triangular R with = Y H , N Z = Y R , Y H Y = I , Z H Z = I , and Y J Z = 0. KZ Use Algorithm 5.1 with starting vector y1 to generate an th step of the generalized -isotropic Arnoldi factorization: = Y H + h+1, y+1 e , KZ Z = Y R . N N ) are convergent, For k = 1, 2, . . . , until wanted eigenpairs of (K, Use Algorithm 5.1 to extend the th step of the generalized -isotropic Arnoldi factorization to the ( + p)th step of the generalized -isotropic Arnoldi factorization: +p = Y+p H+p + h+p+1,+p y+p+1 e , KZ +p Z+p = Y+p R+p . N Use Algorithm 5.2 to reform a new th step of the generalized -isotropic Arnoldi factorization. End Remark 5.1. (i) h+1, is set to zero if |h+1, | < tol(|h, | + |h+1,+1 |) for some stopping tolerance “tol.” (ii) Let (θi , vi ) be an eigenpair of (H , R ), i.e., H vi = θi R vi , and let zi = Z vi = μN z corresponding to the Ritz be a Ritz vector of the eigenproblem Kz value θi . Then from (5.4) and (5.5), we have zi 2 = KZ Z vi 2 i − θi N vi − θi N. Kz = (Y H + h+1, y+1 e )vi − θi Y R vi 2 = Y (H vi − θi R vi ) + h+1, (e vi )y+1 2 = h+1, (e vi )y+1 2 = |h+1, ||e vi |.. 6. Numerical study: Vibration of fast trains. In this section, we shall study the resonance phenomena of a railway track under high frequent excitation forces. We present numerical results of the vibration of fast trains to illustrate the performance of the proposed structure-preserving algorithms in sections 2–5. All numerical experiments are carried out using MATLAB 2006b with the machine precision eps ≈ 2.22 × 10−16 . Research in the vibration of fast trains contributes to the safety of operations of high-speed trains as well as new designs of train bridges, embedded rail structures (ERS), and train suspension systems. Recently, the dynamic response of the vehiclerail-bridge interaction system under different train speed was studied in [25] and a procedure for designing an optimal ERS was proposed in [14]. In both papers, the accurate numerical estimation to the resonance frequencies of the rail plays an important role. However, as mentioned by Ipsen in [5], the classical finite element packages fail to deliver correct resonance frequency for such problems. In this section, we would like to use our structure-preserving algorithms to solve the palindromic QEP (1.1) arising from the spectral modal analysis of rails under periodic excitation forces. In the model of vibration of fast trains, we assume that the rail sections between. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(17) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1582. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. consecutive sleeper bays are identical, the distance between consecutive wheels is the same, and the wheel loads are equal. The rail between two sleepers is modeled by a three-dimensional isotropic elastic solid with linear isoparametric tetrahedron finite elements. Figure 6.1 shows a three-dimensional rail model (see [1] for details).. Fig. 6.1. A three-dimensional rail model.. Based on the ERS design [14], the external force is assumed to be periodic and the displacements of two boundary cross sections of the modeled rail are assumed to have a ratio λ, which is dependent on the excitation frequency of the external force. From the virtual work principle and strain-stress relationship, the governing equation for the displacement vector q involving viscous damping can be formulated by Kq + Dq˙ + M q¨ = f (t), where K, D, and M from the finite element discretization on a uniform mesh satisfy the given boundary conditions. These matrices have the form ⎤ ⎡ 1 1,2:m−1 E E E11 m,m+1 λ ⎥ ⎢ 2:m−1 2:m−1,m E E ⎦ ⎣ E1,2:m−1 λEm,m+1 E2:m−1,m Em,m 2:m−1,m = [0n , . . . , 0n , Em−1,m ], and E 2:m−1 = = [E12 , 0n , . . . , 0n ], E in which E 1,2:m−1 m−1 n×n tridiag Ei−1,i , Ei,i , Ei,i+1 i=2 with Eij ∈ R , i, j = 1, . . . , m + 1. (See [1] for details.) Furthermore, from the spectral modal analysis, we consider q = x eiωt , where ω is the frequency of the external force and x is the corresponding eigenmode. Consequently, we get the palindromic QEP + λA0 + A 1 x (6.1) = 0, λ2 A 1. where Km,m+1 + iωDm,m+1 − ω 2 Mm,m+1 (if i = m, j = 1), 1 A = 0 otherwise, ij 2 Ki,j + iωDi,j − ω Mi,j (if i − 1 ≤ j ≤ i + 1), 0 A = 0 otherwise. ij By consulting the preprocessing procedure (see [4] or [1]) for the deflation of all trivial zero and infinite eigenvalues of (6.1), we arrive to the deflated palindromic QEP 2 λ A1 + λA0 + A1 x = 0. (6.2). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(18) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1583. Example 6.1. We first consider the deflated palindromic QEP (6.2) for high-speed trains and rails. The size of A0 and A1 after deflation is n = 303, and the excitation frequency ω is chosen as 1000. The absolute values of the eigenvalues vary from 10−20 to 1020 . We compute all eigenpairs of Example 6.1 by the SA I, SA II, and QZ algorithm. Note that as shown in section 3, SA II and the URV-based method [20] are mathematically equivalent. In practice, we compare the backward error (relative residual (RRes)) of (1.1) by SA II and the SKURV software [18]. Since SKURV gives only the eigenvalues, the associated eigenvectors are computed from (3.9) and (3.10) by inverse iteration. Numerical results show that the backward errors obtained by SA II and SKURV for Example 6.1 are slightly different. Therefore, in the following computation, we adapt SA II instead of the URV-method. To measure the accuracy of an approximate eigenpair (λ, x) for (6.2), we use the RRes (6.3). RRes ≡. λ2 A 1 x + λA0 x + A1 x 2 . 2 (|λ | A1 F + |λ| A0 F + A1 F ) x 2. As mentioned before, theoretically, the eigenvalues of (6.2) appear in pairs (λ, λ1 ). So, if we sort the eigenvalues in the ascending order by modulus, the product of the ith and (2n + 1 − i)th sorted eigenvalues should be one. Therefore, we define the reciprocities of computed eigenvalues by (6.4). |λi λ2n+1−i − 1|, i = 1, . . . , n.. The RRes of the computed eigenpairs by the SA I, SA II, and QZ algorithm for the eigenvalues with absolute values in [10−20 , 1020 ] and ω = 1000 are shown in Figure 6.2. For eigenvalues with small modulus, the SA I performs much better than the SA II and the QZ algorithm. For eigenvalues near the unit circle or with large modulus, all three algorithms have similar accuracy. The important reciprocity property of eigenvalues is shown in Figure 6.3. Clearly, SA I and SA II preserve the essential reciprocity property as expected, while the QZ algorithm has only less than 12 pairs of computed eigenvalues near the unit circle with reciprocity near zero (≈ 1.17 × 10−12 ). The average and maximal values of all reciprocities are 0.220 and 1.006, respectively. Next, we apply the SA I, SA II, and QZ algorithm to the palindromic QEP with various excitation frequency ω. Figure 6.4 shows the RRes of all computed eigenpairs with eigenvalues in [10−20 , 1020 ] by the three algorithms for 100 different ω’s uniformly chosen from 50 to 5000. We see that the RRes of the SA I are better than those of the SA II and the QZ algorithm for all ω’s. Example 6.2. We now consider the palindromic QEP (6.1) for high-speed trains and rails, with n, the size of A0 and A1 , being 5757. Computational cost. Before showing our numerical results computed by the SHIRA and GSHIRA, we compare the computational costs of one step of the -isotropic Arnoldi process and the implicitly restarted step in each algorithm. In one step of the SHIRA, it requires one matrix-vector product for B, and 3j −1 , by inner products and saxpy operations with vector length 2n. Since B = N1−1 KN 2 the definitions of K and (N1 , N2 ) in (4.2a) and (4.4), the matrix-vector of B requires solving 2 linear systems, 4 and 2 matrix-vector products for A1 and A0 , respectively, and 6 saxpy operations with vector length n.. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(19) 1584. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. SA_I QZ SA_II. Relative residuals of eigenpairs. -10 10. -12 10. -14 10. -16 10. -18 10 -20 10. -10 10. 0. 10. 10 | λ|. 20. 10. 10. Fig. 6.2. RRes of Example 6.1 (ω = 1000). 2. 10. 0. 10. reciprocity of eigenvalues. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. -8 10. -2 10 -4 10 -6 10 -8 10 -10 10 -12 10 -20 10. -10 10. 0. 10 |λ |. 10. 10. 20. 10. i. Fig. 6.3. Reciprocities of computed eigenvalues produced by the QZ algorithm (ω = 1000).. In one step of the GSHIRA, solving zj requires solving 2 linear systems, 2 matrixvector products of A0 and A1 , and 6 saxpy operations with vector length n; computing zj requires 2j − 1 inner products and saxpy operations with vector length 2n;. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(20) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. QZ Algorithm. 1585. SA_I. SA_II/URV. Fig. 6.4. The RRes of eigenvalues vs. ω. Table 6.1 Computational cost of one step of the -isotropic Arnoldi process in the SHIRA and GSHIRA algorithms.. Solving linear system Matrix-vector product for A1 Matrix-vector product for A0 Inner products Saxpy operations. SHIRA 2 4 2 6j 6j + 6. GSHIRA 2 4 2 8j − 2 8j + 4. computing yj+1 requires 2 matrix-vector products for A1 , and 2j inner products and saxpy operations with vector length 2n. We summarize the above computational costs in Table 6.1. The vector length of the inner products and saxpy operations in Table 6.1 is equal to n. On the other hand, the implicitly restarted steps in the SHIRA and GSHIRA require 2( + p − 1)p and 4( + p − 1)p saxpy operations with vector length 2n, respectively. Comparing one -isotropic Arnoldi step with one implicitly restarted step, the GSHIRA algorithm is slightly more expensive than the SHIRA algorithm. Accuracy of eigenpairs. We now compare the numerical results computed by the SHIRA and GSHIRA algorithms. Here, λω,1 , . . . , λω,10 denote target eigenvalues, and we set = 10, p = 20 in the implicitly restarted step for each algorithm. 1 ,x i ) for i = 1, . . . , 10 are shown in Figure 6.5, The RRes of (λω,i , xi ) and ( λω,i where xi and x i are the corresponding computed eigenvectors. In (a) and (b) of Figure 6.5, we show those RRes for frequency ω = 50 and ω = 2000, respectively. The notations “Δ” and “×” denote the results computed by the SHIRA and GSHIRA. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(21) 1586. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN (a). (b). -11 10. -9 10. -12 10. -10 10 -13 10. RRes. RRes. -11 10 -12 10 -13 10. -14 10. -15 10. -14 10 -16 10. -15 10 -16 10 -5 10. 0. -17 10 -5 10. 5. 10 | λ|. 10. 0. 5. 10 | λ|. 10. Fig. 6.5. The RRes of the eigenpairs computed by the SHIRA and GSHIRA algorithms. The notations “Δ” and “×” denote the results computed by the SHIRA and GSHIRA algorithms, respectively. In (a) and (b), the frequency ω is equal to 50 and 2000, respectively.. (a.1). (a.2). (b.1). (b.2). 20. 20. 20. 20. 18. 18. 18. 18. 16. 16. 16. 16. 14. 14. 14. 14. 1.0e-9. 10 8. 12 10 8. Stacked value. 12. Stacked value. Stacked value. 1.0e-10. Stacked value. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. -8 10. 12 10 8. 12 1.0e-11. 10 1.0e-12. 8. 6. 6. 6. 6. 1.0e-13. 4. 4. 4. 4. 1.0e-14. 2. 2. 2. 2. 0. 0. 100. ω. 300. 500. 0 500. 3000 ω. 5000. 0. 0. 100. ω. 300. 500. 0 500. 1.0e-15. 3000 ω. 5000. Fig. 6.6. The stacked bars of ω,k for k = 1, . . . , 7 with different ω. For each ω, all ω,k for k = 1, . . . , 7 are stacked to form a vertical bar with ordering ω,1 , ω,2 , . . . , ω,7 . Each bar is multicolored and the color corresponds to distinct ω,k . The color bar in the right position shows the relationship between color and interval Ik , which corresponds to ω,k . The results in (a) and (b) are computed by the SHIRA and GSHIRA algorithms, respectively.. algorithms, respectively. From these results, we see that the reciprocity property of the eigenvalues are preserved in both algorithms, but the accuracy of the eigenpairs computed by the GSHIRA algorithm is obviously better than that by the SHIRA algorithm. In order to give an overall comparison between the two algorithms, we compute 1 ,x i ) for i = 1, . . . , 10 with ω = 5, 10, 15, . . . , 500 the eigenpairs (λω,i , xi ) and ( λω,i and ω = 550, 600, 650, . . . , 5000. We analyze the distribution of the corresponding 20 RRes with respect to ω. We partition the interval (0, 10−9 ) into seven subintervals I1 = (0, 10−15 ], I2 = (10−15 , 10−14 ], . . . , I7 = (10−10 , 10−9 ). For fixed ω, let ω,k be the number of the RRes which belongs to the interval Ik for k = 1, . . . , 7. In Figure 6.6, for each ω, all ω,k , k = 1, . . . , 7, are stacked to form a vertical bar with ordering. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(22) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP (a). -8 10. -9 10. -9 10. -10 10. -10 10. -11 10. -11 10. Average of RRes. Average of RRes. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. -8 10. -12 10 -13 10. -13 10 -14 10. -15 10. -15 10. 0. 100. ω. 300. 500. (b). -12 10. -14 10. -16 10. 1587. -16 10 500. 3000 ω. 5000. Fig. 6.7. The average of RRes for the twelve eigenpairs computed by the SHIRA and GSHIRA algorithms. The notations “Δ” and “×” denote the results computed by the SHIRA and GSHIRA algorithms, respectively.. ω,1 , ω,2 , . . . , ω,7 . The bar height is 20 which is the sum of ω,1 , . . . , ω,7 . Each bar is multicolored and the color corresponds to distinct ω,k . The color bar in the right bottom position of Figure 6.6 shows the relationship between colors and intervals Ik corresponding to ω,k . All stacked bars of ω,k (k = 1, . . . , 7) with ω = 5, 10, 15, . . . , 500 are shown in (a.1) and (b.1) of Figure 6.6 and those with ω = 550, 600, 650, . . . , 5000 are shown in (a.2) and (b.2) of Figure 6.6. The results in (a) and (b) of Figure 6.6 are computed by the SHIRA and the GSHIRA algorithms, respectively. In the above paragraph, we show the distribution of the RRes for different ω for the comparison of the accuracy of the target eigenpairs. From another point of view, we show the average of the RRes for the target eigenpairs with each ω in Figure 6.7. The notations “Δ” and “×” in Figure 6.7 denote the results computed by the SHIRA and GSHIRA algorithms, respectively. From Figures 6.6 and 6.7, we can summarize that the accuracy of the eigenpairs computed by the GSHIRA algorithm are obviously better than that of the SHIRA algorithm for all ω in (0, 5000]. We now try to explain the different accuracies of the two algorithms. One important reason is that the SHIRA algorithm needs to solve a linear system in the extraction method of eigenvectors, while the GSHIRA algorithm needs only vector additions. The accuracy of the extracted eigenvector will be reduced if the condition number of the linear system is large. On the other hand, Theorem A.1 in Appendix A.3 may help explain this phenomenon from the viewpoint of the minimal residual. The accuracy of the eigenpair computed by the GSHIRA algorithm is better than that by the SHIRA algorithm, since the GSHIRA algorithm is a generalized Arnoldi =μ z, while the SHIRA algorithm is an Arnoldi algorithm for algorithm for Kz N −1 −1 y. N1 KN2 y = μ 7. Conclusions. In this paper, we first transform a palindromic QEP to a skew-Hamiltonian pencil by the (S + S −1 )-transform. Then, we extend Patel’s ap-. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(23) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1588. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. proach to solve the -skew-Hamiltonian pencil efficiently. We have also developed a structure-preserving generalized -skew-Hamiltonian implicitly restarted Arnoldi method (GSHIRA) for solving the large sparse -skew-Hamiltonian pencil. Numerical results show that the accuracy of desired eigenpairs computed by the GSHIRA is better than that computed by the classical SHIRA. The standard algorithms proposed in this paper are numerically stable for solving palindromic QEPs. In the future, we are motivated to develop structure-preserving algorithms for solving the antipalindromic QEP λ2 A 1 + λA0 − A1 with A0 = −A0 , efficiently. Appendix. A.1. In this section we list pseudocodes of Step 2 in Algorithm 2.1. In the following, givensl(α,β,i) returns a Givens rotation G such that G[ αβ ] = γei with γ ∈ C; givensr(α,β,i) returns a Givens rotation G such that [α β]G = γe i with γ ∈ C. The functions qr(A) and ql(A) perform the standard QR and QL factorizations. Step 2 in Algorithm 2.1. function [K, N , Q, Z] = rbutf(K, N ) Input: Matrices K, N in the form (2.4). Output: Unitary Q, Z and K, N of the form (2.5), where K and N are overwritten by QKZ and QN Z, respectively. 01: [Q1 , R] ← qr(N (1 : n, 1 : n)) 02: Q ← diag(QH 1 , In ) ¯ 1) 03: Z ← diag(In , Q 04: K ← QKZ 05: N ← QN Z 06: for j = 1 : n − 2 07: for k = j + 1 : n − 1 08: % annihilate K(n + k, j) by Givens rotation in (n + k, n + k + 1) plane 09: G ←givensl(K(n + k, j), K(n + k + 1, j), 2) 10: Q(n + k : n + k + 1, :) ← GQ(n + k : n + k + 1, :) 11: K(n + k : n + k + 1, :) ← GK(n + k : n + k + 1, :) 12: N (n + k : n + k + 1, :) ← GN (n + k : n + k + 1, :) 13: Z(:, k : k + 1) ← Z(:, k : k + 1)G 14: K(:, k : k + 1) ← K(:, k : k + 1)G 15: N (:, k : k + 1) ← N (:, k : k + 1)G 16: % annihilate N (k + 1, k) by Givens rotation in (k, k + 1) plane 17: G ←givensl(N (k, k), N (k + 1, k), 1) 18: Q(k : k + 1, :) ← GQ(k : k + 1, :) 19: K(k : k + 1, :) ← GK(k : k + 1, :) 20: N (k : k + 1, :) ← GN (k : k + 1, :) 21: Z(:, n + k : n + k + 1) ← Z(:, n + k : n + k + 1)G 22: K(:, n + k : n + k + 1) ← K(:, n + k : n + k + 1)G 23: N (:, n + k : n + k + 1) ← N (:, n + k : n + k + 1)G 24: end 25: % annihilate N (2n, j) by Givens rotation in (n, 2n) plane 26: G ←givensl(N (n, j), N (2n, j), 1) 27: Q([n 2n], :) ← GQ([n 2n], :) 28: K([n 2n], :) ← GK([n 2n], :) 29: N ([n 2n], :) ← GN ([n 2n], :) 30: Z(:, [n 2n]) ← Z(:, [n 2n])GH 31: K(:, [n 2n]) ← K(:, [n 2n])GH. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(24) 1589. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 32: N (:, [n 2n]) ← N (:, [n 2n])GH 33: for k = n : −1 : j + 2 34: % annihilate K(k, j) by Givens rotation in (k − 1, k) plane 35: G ←givensl(K(k − 1, j), K(k, j), 1) 36: Q(k − 1 : k, :) ← GQ(k − 1 : k, :) 37: K(k − 1 : k, :) ← GK(k − 1 : k, :) 38: N (k − 1 : k, :) ← GN (k − 1 : k, :) 39: Z(:, n + k − 1 : n + k) ← Z(:, n + k − 1 : n + k)G 40: K(:, n + k − 1 : n + k) ← K(:, n + k − 1 : n + k)G 41: N (:, n + k − 1 : n + k) ← N (:, n + k − 1 : n + k)G 42: % annihilate N (k, k − 1) by Givens rotation in (k − 1, k) plane 43: G ←givensr(N (k, k − 1), N (k, k), 2) 44: Q(n + k − 1 : n + k, :) ← G Q(n + k − 1 : n + k, :) 45: K(n + k − 1 : n + k, :) ← G K(n + k − 1 : n + k, :) 46: N (n + k − 1 : n + k, :) ← G N (n + k − 1 : n + k, :) 47: Z(:, k − 1 : k) ← Z(:, k − 1 : k)G 48: K(:, k − 1 : k) ← K(:, k − 1 : k)G 49: N (:, k − 1 : k) ← N (:, k − 1 : k)G 50: end 51: end A.2. To show the extra zeros of the subdiagonals of the submatrices in (3.4), let Hk and Tk be the sets of k ×k upper Hessenberg and triangular matrices, respectively, and let S2k be the set of 2k × 2k -skew symmetric matrices. Denote 0k 2k×2k (A2.1) A2k = A ∈ C P2k with ∈ Hk and ∈ Tk , A ≡ P2k 0k where P2k = [e1 , ek+1 , e2 , ek+2 , . . . , ek , e2k ], 0k 2k×2k R2k = R ∈ C (A2.2) R ≡ P2k. 0k. P2k with. ∈ Tk ,. (A2.3) B2m,2k = {B ∈ C2m×2k |Be1 = Be3 = · · · = Be2k−1 = 0}, 2m,2k = {B ∈ C2m×2k |Be 2 = Be 4 = · · · = Be 2k = 0}, (A2.4) B C2m×2k = {C ∈ C2m×2k |cij = 0, i = 1, . . . , 2m, j = 1, . . . , 2k (A2.5) D2k = {D ∈ C (A2.6). 2k×2k. and (i, j) = (1, 2k)}, |D ∈ S2k with {1, −1, 3, −3, . . . , 2k − 1, −(2k − 1)}. – diagonals being zeros}, 2k×2k D2k = {D ∈ C |D ∈ S2k with {2, −2, 4, −4, . . . , 2k − 2, −(2k − 2)}. (A2.7). – diagonals being zeros}.. After performing the section A.1, for j = 1 and ⎡ 0 × ⎢ × 0 ⎢ ⎢ 0 × ⎢ (2) (A2.8a) K11 := ⎢ 0 0 ⎢ ⎢ .. .. ⎣ . . 0. first and second steps of the N ), it produces 2) on (K, ⎤ ⎡ 0 ··· 0 0 ⎥ × ··· × ⎥ ⎢ 0 ⎢ 0 ··· 0 ⎥ ⎢ × ⎥ (2) := , K ⎢ ⎥ 12 0 ··· 0 ⎥ ⎢ .. ⎣ . .. .. ⎥ . . ⎦ × 0 0 ··· 0. SA I (i.e., Steps 07–50 in. 0 × 0 0 0 .. .. ··· ··· G2n−2. ⎤ × 0 ⎥ ⎥ ⎥ ⎥, ⎥ ⎦. 0. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(25) 1590. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. ⎡ (2) K21. (A2.8b). and. ⎡. (A2.9a). (2) N11. ⎢ ⎢ ⎢ := ⎢ ⎢ ⎣. ⎢ ⎢ ⎢ := ⎢ ⎢ ⎣. × 0 0 .. .. 0 × 0 .. .. 0. 0. 0 0 0 .. .. 0 0 0 0 0 .. .. 0. 0. × 0. H2n−2. ⎡ ⎤ 0 × ⎢ × 0 ⎥ ⎢ ⎥ ⎢ ⎥ (2) ⎥ , N12 := ⎢ 0 ⎢ .. ⎥ ⎣ . ⎦ 0. ··· ··· T2n−2. (2). (2). × 0 × .. .. ⎤ 0 × ⎥ ⎥ ⎥ ⎥, ⎥ ⎦. ··· ···. 0 ×. 02n−2. ×. (2). N22 := (N11 ) ,. N21 := 02n ,. (A2.9b). ⎤ 0 0 ⎥ ⎥ ⎥ (2) (2) ⎥ , K22 := (K11 ) , ⎥ ⎦. ··· ···. where G2n−2 and H2n−2 ∈ S2n−2 and T2n−2 ∈ T2n−2 . Let m = n − k. Suppose after 2k steps (for j = 1, 2, . . . , 2k) the SA I gives. B D2k A2k B2m,2k (2k) (2k) 2m,2k (A2.10a) K11 := , , K12 := 2m,2k C2m,2k 02m B G2m 02k 0 (2k) (2k) (2k) 2m,2k (A2.10b) K21 := , K22 := (K11 ) , 02m,2k H2m and (A2.11a). (2k). N11. :=. R2k 02m,2k (2k). N21. (A2.11b). E 2m,2k T2m. . := 02n ,. (2k). , N12. (2k). N22. . 2k D E2m,2k. :=. E2m,2k 02m. ,. (2k). := (N11 ) ,. 2m,2k , where A2k ∈ A2k , R2k ∈ R2k , C2m,2k ∈ C2m,2k , B2m,2k , E2m,2k ∈ B2m,2k , B E2m,2k ∈ B2m,2k , D2k ∈ D2k , D2k ∈ D2k , G2m , H2m ∈ S2m , and T2m ∈ T2m . By letting k = k + 1 and m = m − 1, we perform the SA I for j = 2k + 1, 2k + 1 and obtain. B D2k A2k B2m ,2k (2k ) (2k ) 2m ,2k (A2.12a) K11 := , , K12 := 2m ,2k C2m ,2k 02m B G2m 02k 0 (2k ) (2k ) (2k ) 2m ,2k , K22 := (K11 ) , (A2.12b) K21 := 02m ,2k H2m and (A2.13a). (A2.13b). (2k ) N11. :=. R2k 02m ,2k (2k ). N21. E 2m ,2k T2m := 02n ,. ,. (2k ) N12. (2k ). N22. :=. 2k D. E2m ,2k. E2m ,2k 02m. ,. (2k ) . := (N11. ) ,. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(26) STRUCTURE-PRESERVING ALGORITHMS FOR PQEP. 1591. Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. where the subblocks in (A2.12)–(A2.13) have the same forms as in (A2.10)–(A2.11) by replacing k and m by k and m , respectively, and satisfy (A2.14a). A2k = Φ 2k A2k Φ2k ,. (A2.14b). R2k = Φ 2k R2k Φ2k ,. D2k = Φ 2k D2k Φ2k , D2k = Φ 2k D2k Φ2k , . where Φ2k = [e1 , . . . , e2k ] with ei ∈ C2k , i = 1, . . . , 2k. By the inductive process above, (3.4) holds with k = n in (A2.12)–(A2.13) and with the superscript “a” in (3.4) being (2n). A.3. Theorem A.1. Let V ∈ Cn×r be a unitary matrix and A, B ∈ Cn×n . Then. AV − BV C 2 ≥ AV − BV P 2. for all C ∈ Cr×r ,. where P = (V H B H BV )−1 (V H B H AV ), or equivalently, P = (U H BV )−1 (U H AV ), where BV = U S is the QR factorization of BV . Proof. Since −C H V H B H BV P = −C H V H B H BV (V H B H BV )−1 (V H B H AV ) = −C H V H B H AV, it follows that (V H AH − C H V H B H )(AV − BV C) = V H AH AV − C H V H B H AV − V H AH BV C + C H V H B H BV C = V H AH AV + (P H − C H )V H B H BV (P − C) − P H V H B H BV P = (V H AH − P H V H B H )(AV − BV P ) + (P H − C H )V H B H BV (P − C). Obviously, (P H − C H )V H B H BV (P − C) is semidefinite. Then by Weyl’s theorem, we have λj ((AV − BV C)H (AV − BV C)) ≥ λj ((AV − BV P )H (AV − BV P )), j = 1, . . . , n. Hence. AV − BV C 2 ≥ AV − BV P 2 , since G 22 = λmax(GH G). REFERENCES [1] E. K.-W. Chu, T.-M. Hwang, W.-W. Lin, and C.-T. Wu, Vibration of fast trains, palindromic eigenvalue problems and structure-preserving doubling algorithms, J. Comput. Appl. Math., 219 (2008), pp. 237–252. [2] J. J. Hench and A. J. Laub, Numerical solution of the discrete-time periodic Riccati equation, IEEE Trans. Automat. Control, 39 (1994), pp. 1197–1210. [3] A. Hilliges, Numerische L¨ osung von quadratischen eigenwertproblemen mit Anwendungen in der Schiendynamik, Master’s thesis, Technical University Berlin, Berlin, Germany, 2004. [4] A. Hilliges, C. Mehl, and V. Mehrmann, On the solution of palindromic eigenvalue problems, in Proceedings of the 4th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS), Jyv¨ askyl¨ a, Finland, 2004. [5] I. C. F. Ipsen, Accurate eigenvalues for fast trains, SIAM News, 37, SIAM, Philadelphia, 2004.. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(27) Downloaded 04/30/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php. 1592. TSUNG-MING HUANG, WEN-WEI LIN, AND JIANG QIAN. [6] D. Kressner, Numerical methods for general and structured eigenvalue problems, Lect. Notes Comput. Sci. Eng. 46, Springer, Berlin, 2005. [7] R. B. Lehoucq and D. C. Sorensen, Deflation techniques for an implicitly restarted Arnoldi iteration, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 789–821. [8] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide: Solution of LargeScale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998. [9] W.-W. Lin, A new method for computing the closed-loop eigenvalues of a discrete-time algebraic Riccati equation, Linear Algebra Appl., 96 (1987), pp. 157–180. [10] 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. [11] 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., to appear. [12] 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. [13] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann, Vector spaces of linearizations for matrix polynomials, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 971–1004. [14] V. L. Markine, A. P. de Man, S. Jovanovic, and C. Esveld, Optimal design of embedded rail structure for high-speed railway lines, in Railway Engineering 2000, 3rd International Conference, London, 2000. [15] V. Mehrmann and D. Watkins, Structure-preserving methods for computing eigenpairs of large sparse skew-Hamiltonian/Hamiltonian pencils, SIAM J. Sci. Comput., 22 (2001), pp. 1905–1925. [16] R. V. Patel, On computing the eigenvalues of a symplectic pencil, Linear Algebra Appl., 188 (1993), pp. 591–611. [17] Y. Saad, Numerical Methods for Large Eigenvalue Problems, Manchester University Press, Manchester, UK, 1992. ¨ der, SKURV: A Matlab Toolbox for the Skew URV Decomposition of a Matrix Triple, [18] C. Schro http://www.math.tu-berlin.de/∼schroed/Software/skurv. ¨ der, A QR-like Algorithm for the Palindromic Eigenvalue Problem, Technical report, [19] C. Schro Preprint 388, TU Berlin, Berlin, Germany, 2007. ¨ der, URV decomposition based structured methods for palindromic and even eigen[20] C. Schro value problems, Technical report, Preprint 375, TU Berlin, Berlin, Germany, 2007. [21] D. C. Sorensen, Implicit application of polynomial filters in a k-step Arnoldi method, SIAM J. Matrix Anal. Appl., 13 (1992), pp. 357–385. [22] G. W. Stewart, Matrix Algorithms, Volume II: Eigensystems, SIAM, Philadelphia, 2001. [23] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev., 43 (2001), pp. 235–286. [24] D. S. Watkins, The Matrix Eigenvalue Problem: GR and Krylov Subspace Methods, SIAM, Philadelphia, 2007. [25] Y. S. Wu, Y. B. Yang, and J. D. Yau, Three-dimensional analysis of train-rail-bridge interaction problems, Vehicle Syst. Dyn., 35 (2001), pp. 1–35.. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited..

(28)