[PPT] - Rational Krylov for Real Pencils with Complex Eigenvalues Axel Ruhe PowerPoint Presentation

SLIDE 1

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

SLIDE 2

Introduction

We consider a linear time-invariant dynamical system in state space form E ˙ x(t) = Ax(t) + bu(t), y(t) = cx(t),

(1)

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

SLIDE 3

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

SLIDE 4

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

SLIDE 5

Rational Krylov

The rational Krylov algorithm is a shift invert Arnoldi decomposition (A − σE)−1EVj = Vj+1Hj+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

SLIDE 6

Basic recursion, one step

Consider step j, (A − σjE)−1Evj = Vj+1hj , Multiply by (A − σjE) and get Evj = (A − σjE)Vj+1hj , and separate A and E terms EVj+1(ej + hjσj) = AVj+1hj ,

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 6/33

SLIDE 7

Basic recursion, two matrices

Put earlier columns up front and get, AVj+1Hj+1,j = EVj+1(I + Hj+1,jdiag(σi)) with two Hessenberg matrices, H consisting of the Gram Schmidt coefficients, and K = I + Hdiag(σi).

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 7/33

SLIDE 8

Arnoldi formulation

We may now formulate this as an Arnoldi iteration with another starting vector w1 chosen from the Krylov space spanned by Vj+1 by QR factorizing either H or K. Choose H, direct iteration. AVj+1Qj+1Rj+1,j = EVj+1Qj+1(Q′

j+1Kj+1,j)

AWjRj,j = EWj+1Mj+1,j

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 8/33

SLIDE 9

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 w1

but keep the residuals in the direction wj+1.

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 9/33

SLIDE 10

Arnoldi formulation, cont

Apply the QZ algorithm to the pencil (Mj,j, Rj,j). A Ritz pair (θ, y) is given by Mj,jsj = Rj,jsjθ y = WjRj,jsj Its residual is, Ay − Eyθ = Ewj+1mj+1,.sj in the direction of the next basis vector Ewj+1 and length given by last row of M and the eigenvector sj.

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 10/33

SLIDE 11

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)−1Evj The matrix diag(σi) in K = I + Hdiag(σi) will have 2 × 2 blocks.

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 11/33

SLIDE 12

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)−1Evj adding one column to H each time.

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 12/33

SLIDE 13

Real matrices, single vector 2

We get BVj = Vj+1Hj+1,j where B = B+ = real((A − σE)−1E)

r

B = B− = imag((A − σE)−1E) .

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 13/33

SLIDE 14

Real matrices, single vector 3

Note that B+ = 1 2

(A − σE)−1) + (A − σE)−1)
E

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

SLIDE 15

Real matrices, single vector 4

Note that if we partition σ = ρ + iθ µ+ = λ(B+) = 1 2( 1 λ − σ + 1 λ − σ) = 1 2( λ − σ + λ − σ (λ − σ)(λ − σ)) = λ − ρ (λ2 − 2ρλ + |σ|2) µ− = λ(B−) = θ (λ2 − 2ρλ + |σ|2)

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 15/33

SLIDE 16

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

SLIDE 17

Hopf, eigenvalue approximations:

−30 −25 −20 −15 −10 −5 5 −15 −10 −5 5 10 15 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 17/33

SLIDE 18

Hopf, Transformed eigenvalues:

−2 −1.5 −1 −0.5 0.5 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 Transformed spectrum Pencil B+ B−

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 18/33

SLIDE 19

Hopf, residuals j = 20:

2 4 6 8 10 12 14 16 18 20 10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 10

2

10

4

estimated residuals computed residuals

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 19/33

SLIDE 20

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

SLIDE 21

Inlet, transfer function:

5 10 15 20 0.5 1 1.5 2 ω |H(iω)| Inlet: Transfer functions, jmax=12 full reduced

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 21/33

SLIDE 22

Inlet, eigenvalue approximations:

−10 −8 −6 −4 −2 −25 −20 −15 −10 −5 5 10 15 20 25 Re{λ} Im{λ} Inlet: Eigenvalues, jmax=12 full reduced pole

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 22/33

SLIDE 23

Inlet, more eigenvalues:

−20 −15 −10 −5 50 100 150 Re{λ} Im{λ} all λ RKeig: λ RKeig: σ

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 23/33

SLIDE 24

Inlet, real algorithm, zq residuals:

5 10 15 20 25 30 35 40 45 50 10

−20

10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

estimated residuals computed residuals

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 24/33

SLIDE 25

Inlet, real algoritm σ = 0 eigenvalues:

−50 50 100 150 −60 −40 −20 20 40 60 80 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 25/33

SLIDE 26

Inlet, real algorithm eigenvalues:

−15 −10 −5 5 −10 −8 −6 −4 −2 2 4 6 8 10 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 26/33

SLIDE 27

Inlet, real algorithm σ = 6i eigenvalues

−15 −10 −5 5 2 4 6 8 10 12 14 16 18 20 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 27/33

SLIDE 28

Inlet, σ = 5 + 5i eigenvalues of H:

−0.25 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 0.25 −0.25 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 0.25

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 28/33

SLIDE 29

Inlet, real algorithm σ = 20i eigenvalue

−25 −20 −15 −10 −5 5 10 15 20 25 30 35 40 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 29/33

SLIDE 30

Inlet, real algorithm σ = 40i eigenvalue

−25 −20 −15 −10 −5 5 25 30 35 40 45 50 55 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 30/33

SLIDE 31

Inlet, real algorithm σ = 75i eigenvalue

−25 −20 −15 −10 −5 5 55 60 65 70 75 80 85 Eigenvalue approximations approximations close appr converged comparison shift

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 31/33

SLIDE 32

Inlet, real algorithm qz residuals:

20 40 60 80 100 120 10

−22

10

−20

10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

estimated residuals computed residuals

Talk at RANMEP2008 Taipei, January 8 2008, Axel Ruhe – p. 32/33

SLIDE 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