Rational Krylov for Real Pencils with Complex Eigenvalues Axel Ruhe School of Computer Science and Communications, Royal Institute of Technology, SE-10044 Stockholm, Sweden ( ruhe@kth.se ). Work with Henrik Olsson, Zenon Matthews and Henrik Holst. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 1/33
Introduction We consider a linear time-invariant dynamical system in state space form E ˙ x ( t ) = Ax ( t ) + bu ( t ) , (1) y ( t ) = cx ( t ) , connecting input u ( t ) and output y ( t ) via the state x ( t ) . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 2/33
Transfer function Apply Laplace transformation, assuming zero initial conditions, and get Y ( s ) = c ( sE − A ) − 1 bU ( s ) ≡ H ( s ) U ( s ) . (2) The transfer function H ( s ) gives the output Y ( s ) as a function of the input U ( s ) at the complex frequency s . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 3/33
Eigenvalue Problem The poles of the transfer function H ( s ) are eigenvalues of the linear matrix pencil ( A − λE ) x = 0 . (3) Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 4/33
Rational Krylov The rational Krylov algorithm is a shift invert Arnoldi decomposition ( A − σE ) − 1 EV j = V j +1 H j +1 ,j , where the shift σ is changed in some of the steps j . The resulting matrix H is of Hessenberg form. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 5/33
Basic recursion, one step Consider step j , ( A − σ j E ) − 1 Ev j = V j +1 h j , Multiply by ( A − σ j E ) and get Ev j = ( A − σ j E ) V j +1 h j , and separate A and E terms EV j +1 ( e j + h j σ j ) = AV j +1 h j , Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 6/33
Basic recursion, two matrices Put earlier columns up front and get, AV j +1 H j +1 ,j = EV j +1 ( I + H j +1 ,j diag ( σ i )) with two Hessenberg matrices, H consisting of the Gram Schmidt coefficients, and K = I + H diag ( σ i ) . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 7/33
Arnoldi formulation We may now formulate this as an Arnoldi iteration with another starting vector w 1 chosen from the Krylov space spanned by V j +1 by QR factorizing either H or K . Choose H , direct iteration. AV j +1 Q j +1 R j +1 ,j = EV j +1 Q j +1 ( Q ′ j +1 K j +1 ,j ) AW j R j,j = EW j +1 M j +1 ,j Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 8/33
Arnoldi formulation, cont The matrix M is full, but can be brought to Hessenberg form by a backwards Householder algorithm. That will modify the starting vector w 1 but keep the residuals in the direction w j +1 . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 9/33
Arnoldi formulation, cont Apply the QZ algorithm to the pencil ( M j,j , R j,j ) . A Ritz pair ( θ, y ) is given by M j,j s j = R j,j s j θ y = W j R j,j s j Its residual is, Ay − Eyθ = Ew j +1 m j +1 ,. s j in the direction of the next basis vector Ew j +1 and length given by last row of M and the eigenvector s j . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 10/33
Real matrices, double vector We have developed a double shift variant, that gives a real Hessenberg pair when the original A, E are real, but the shifts σ and the eigenvalues λ are complex. In each step add two vectors, given by the real and imaginary parts of r = ( A − σE ) − 1 Ev j The matrix diag ( σ i ) in K = I + H diag ( σ i ) will have 2 × 2 blocks. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 11/33
Real matrices, single vector 1 We are now testing the Parlett and Saad original variant, where only one vector is added each step. In each step, take real or imaginary part of r = ( A − σE ) − 1 Ev j adding one column to H each time. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 12/33
Real matrices, single vector 2 We get BV j = V j +1 H j +1 ,j where B = B + = real (( A − σE ) − 1 E ) or B = B − = imag (( A − σE ) − 1 E ) . Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 13/33
Real matrices, single vector 3 Note that B + = 1 ( A − σE ) − 1 ) + ( A − σE ) − 1 ) � � E 2 and similarly for B − . It will get the same eigenvectors but other complex conjugate eigenvalues. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 14/33
Real matrices, single vector 4 Note that if we partition σ = ρ + iθ µ + = λ ( B + ) = 1 1 1 2( λ − σ + λ − σ ) = 1 2( λ − σ + λ − σ ( λ − σ )( λ − σ )) λ − ρ = ( λ 2 − 2 ρλ + | σ | 2 ) θ µ − = λ ( B − ) = ( λ 2 − 2 ρλ + | σ | 2 ) Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 15/33
Numerical Example, Hopf bifurcation: This small problem was reported by Parlett and Saad. • Reaction and transport in a tubular reactor. • Two one dim equations n = 200 • Most eigenvalues real and negative. • Seek eigenvalues close to imaginary axis. Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 16/33
Hopf, eigenvalue approximations: Eigenvalue approximations 15 approximations close appr converged comparison 10 shift 5 0 −5 −10 −15 −30 −25 −20 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 17/33
Hopf, Transformed eigenvalues: Transformed spectrum 1 0.8 Pencil B+ B− 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −2 −1.5 −1 −0.5 0 0.5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 18/33
Hopf, residuals j = 20 : 4 10 estimated residuals computed residuals 2 10 0 10 −2 10 −4 10 −6 10 −8 10 −10 10 −12 10 −14 10 −16 10 0 2 4 6 8 10 12 14 16 18 20 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 19/33
Numerical Examples, Inlet: • Supersonic engine inlet . • Two dimensional Euler equations n = 11730 • Seek eigenvalues close to imaginary axis Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 20/33
Inlet, transfer function: Inlet: Transfer functions, jmax=12 2 full 1.5 reduced |H(i ω )| 1 0.5 0 0 5 10 15 20 ω Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 21/33
Inlet, eigenvalue approximations: Inlet: Eigenvalues, jmax=12 25 20 15 10 5 Im{ λ } 0 −5 full −10 reduced −15 pole −20 −25 −10 −8 −6 −4 −2 0 Re{ λ } Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 22/33
Inlet, more eigenvalues: 150 all λ RKeig: λ RKeig: σ 100 Im{ λ } 50 0 −20 −15 −10 −5 0 Re{ λ } Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 23/33
Inlet, real algorithm, zq residuals: −2 10 estimated residuals computed residuals −4 10 −6 10 −8 10 −10 10 −12 10 −14 10 −16 10 −18 10 −20 10 0 5 10 15 20 25 30 35 40 45 50 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 24/33
Inlet, real algoritm σ = 0 eigenvalues: Eigenvalue approximations 80 approximations close appr converged comparison 60 shift 40 20 0 −20 −40 −60 −50 0 50 100 150 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 25/33
Inlet, real algorithm eigenvalues: Eigenvalue approximations 10 8 6 4 2 0 −2 −4 approximations close appr −6 converged comparison shift −8 −10 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 26/33
Inlet, real algorithm σ = 6 i eigenvalues Eigenvalue approximations 20 approximations 18 close appr converged comparison 16 shift 14 12 10 8 6 4 2 0 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 27/33
Inlet, σ = 5 + 5 i eigenvalues of H : 0.25 0.2 0.15 0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 −0.25 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 28/33
Inlet, real algorithm σ = 20 i eigenvalue Eigenvalue approximations 40 approximations close appr converged 35 comparison shift 30 25 20 15 10 −25 −20 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 29/33
Inlet, real algorithm σ = 40 i eigenvalue Eigenvalue approximations 55 approximations close appr 50 converged comparison shift 45 40 35 30 25 −25 −20 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 30/33
Inlet, real algorithm σ = 75 i eigenvalue Eigenvalue approximations 85 approximations close appr 80 converged comparison shift 75 70 65 60 55 −25 −20 −15 −10 −5 0 5 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 31/33
Inlet, real algorithm qz residuals: −2 10 estimated residuals computed residuals −4 10 −6 10 −8 10 −10 10 −12 10 −14 10 −16 10 −18 10 −20 10 −22 10 0 20 40 60 80 100 120 Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 32/33
Further developments: • Follow nonlinear iteration Joakim Möller Special question: How nonlinear is the iteration? Estimate R in A ( λ ) = λ − µ σ − µA ( σ )+ λ − σ µ − σA ( µ )+( λ − σ )( λ − µ ) R ( λ ) • Iteration instead of factorization Sonia Gupta • Model reduction Henrik Olsson Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 33/33
Recommend
More recommend
