Finding Rightmost Eigenvalues of Large, Sparse, Nonsymmetric Parameterized Eigenvalue Problems

1. Finding Rightmost Eigenvalues of Large, Sparse, Nonsymmetric Parameterized Eigenvalue Problems Minghao Wu AMSC Program mwu@math.umd.edu Advisor: Dr. Howard Elman Department of Computer Science elman@cs.umd.edu

2. Motivation • To determine the stability of the linearized system of the form: • The steady state solution x* is • - stable, if all the eigenvalues of Ax = λBx have • negative real parts; • - unstable, otherwise.

3. Problem Statement To find the rightmost eigenvalues of: • where matrices A and B are • Real N-by-N • Large • Sparse • Nonsymmetric • Depend on one or several parameters

4. Method • Eigensolver: Arnoldi Algorithm - Iterative method - Based on Krylov subspace: - ComputesArnoldi decomposition: where: [Uk uk+1]: an orthonormal basis of Hk: k-by-k upper Hessenberg matrix, k « N βk: scalar, ek: k-by-1 vector [0 0 … 0 1]’

5. Arnoldi Algorithm (continue) • Eigenvaluesof Hkapproximate eigenvalues of A Premultiply previous equation by transpose(Uk): Let (λ,z) be an eigenpair of Hk, then: Residual As k increases, ||Residual|| will decrease. When k = N, Residual = 0.

6. Matrix Transformation • Motivation - Arnoldi Algorithm cannot solve generalized eigenvalue problem - It converges to well – separated extremal eigenvalues, not rightmost eigenvalues • Shift – Invert Transformation

7. Matlab Code of Arnoldi Method • Arnoldi Algorithm (with Shift – Invert matrix transformation) routine: [v,X,U,H]=SI_Arnoldi(A,B,k,sigma) Input: A, B: matrix A and B in Ax = λBx k: number of eigenpairs wanted sigma: the shift σin shift – invert matrix transformation Output: v: a vector of k computed eigenvalues X: k eigenvectors associated with the eigenvalues U: the Krylov basis Uk+1 H: the upper Hessenberg matrix Hk

8. Test Problem Olmstead Model (see Olmstead et al (1986)): with boundary conditions: This model represents the flow of a layer of viscoelastic fluid heated from below. u: the speed of the fluid S: related to viscoelastic forces b,c: scalars, R: scalar, Rayleigh number

9. Test Problem (continue) • Discretize the model with finite differences - grid – size: h = 1/(N/2) - (4) can be written as dy/dt = f(y) with • Evaluate the Jacobian matrix A = df/dy at steady state • solution y* • - N = 1000, b = 2, c = 0.1, R = 0.6 • - y* = 0 • - A = df/dy(y*) is a nonsymmetric sparse matrix with bandwidth 6

10. Test Problem (continue) Computational Result: • Rightmost eigenvalues: λ1,2 = 0 ± 0.4472i • Residual: ||Axi - λi xi|| = 8.4504e-012, i=1,2 • The result agrees with the literature.

11. Implicitly Restarted Arnoldi (IRA) • Motivation - Large k is not practical Example: • When B is singular, Arnoldi algorithm may give rise to spurious • eigenvalues

12. IRA (continue) • Basic idea of Implicitly Restarted Arnoldi Filter out the unwanted eigendirections from the starting vector by using the most recent spectrum information and a clever filtering technique • IRA steps 1. Compute m eigenpairs (k<m«N) by Arnoldi method with starting vector u1 2. Filter out the m-k unwanted eigendirections from u1 (Key Technique: shifted QR algorithm) 3. Restart the process with filtered starting vector till the k eigenvalues of interest converge

13. Test Problem K: 200-by-200 matrix, full rank; C: 200-by-100 matrix, full rank; M: 200-by-200 matrix, full rank. Eigenvalue problem with this kind of block structure appears in the stability analysis of steady state solution of Navier – Stokes equations for incompressible flow.

14. Test Problem (continue) • Use Matlab function “rand” to generate K, C, M • -2.7377 ≤ Re(λ) ≤ 49.9129 • Find out 10 rightmost eigenvalues • Use the IRA code written by Fei Xue

15. Test Problem (continue) Computational Result: (shift σ= 60)

16. Future Work (AMSC 664) • Solve the third test problem • Implement iterative solvers for linear systems