Secular Equations: Solving and Understanding the Eigenvalue Problem, Study notes of Latin

Download Secular Equations: Solving and Understanding the Eigenvalue Problem and more Study notes Latin in PDF only on Docsity!

Chapter 12

Solving secular equations

G´erard MEURANT

January-February, 2012

(^1) Examples of Secular Equations

(^2) Secular equation solvers

(^3) Numerical experiments

Eigenvalues of a Tridiagonal Matrix

We look for an eigenvalue λ and an eigenvector x =

y ζ

)T

of Jk+1 where y is a vector of dimension k and ζ is a real number. Then Jk y + ηk ζek^ = λy ηk yk + αk+1ζ = λζ where yk is the last component of y , αj , j = 1,... , k + 1 are the diagonal entries of Jk+1 and ηj , j = 1,... , k are the subdiagonal entries By eliminating the vector y from these two equations we have

(αk+1 − η^2 k (ek^ )T^ (Jk − λI )−^1 ek^ )ζ = λζ

αk+1 − η k^2

∑k

j=

(ξj )^2 θj − λ

= λ

where ξj = zkj is the kth (i.e., last) component of the jth eigenvector of Jk the θj ’s are the eigenvalues of Jk , that is, the Ritz values Too obtain the eigenvalues of Jk+1 from those of Jk we have to solve

f (λ) = λ − αk+1 + η^2 k

∑k

j=

ξ^2 j θj − λ

The secular function f has poles at the eigenvalues (Ritz values) of Jk for λ = θj = θ( j k), j = 1... , k

Modification by a Rank-One Matrix

Assume that we know the eigenvalues of a matrix A and we would like to compute the eigenvalues of a rank-one modification of A We have Ax = λx where we know the eigenvalues λ and we want to compute μ such that (A + ccT^ )y = μy where c is a given vector (not orthogonal to an eigenvector of A) Clearly μ is not an eigenvalue of A Therefore A − μI is nonsingular and we obtain an equation for μ

y = −(A − μI )−^1 ccT^ y

Multiplying by cT^ to the left,

cT^ y = −cT^ (A − μI )−^1 ccT^ y

The secular equation is

1 + cT^ (A − μI )−^1 c = 0

Using the spectral decomposition of A = QΛQT^ with Q orthogonal and Λ diagonal and z = QT^ c

∑^ n

j=

(zj )^2 λj − μ

where λj are the eigenvalues of A

Using the spectral decomposition of A = QΛQT^ and d = QT^ c

f (λ) =

∑^ n

j=

d j^2 λj − λ

There are n − 1 solutions to the secular equation When we have the values of λ that are solutions, we use the constraint xT^ x = 1 to remark that

xT^ x = μ^2 cT^ (A − λI )−^2 c = 1

Therefore, μ^2 =

cT^ (A − λI )−^2 c and x = −μ(A − λI )−^1 c

− (^10) − 1 0 1 2 3 4 5 6 7

− 8

− 6

− 4

− 2

0

2

4

6

8

10

Example of secular function

Assume that all the eigenvalues λi of A are simple If λ is equal to one of the eigenvalues, say λj , we must have dj = 0 For all i 6 = j we have

yi =

di λi − λ Then, there is a solution or not, whether we have ∑

i 6 =j

di λi − λj

= α^2

or not If λ is not an eigenvalue of A, the inverse of A − λI exists and we obtain the secular equation

cT^ (A − λI )−^2 c = α^2

which is written as

f (λ) =

∑^ n

i=

di λi − λ

− α^2 = 0

− (^10) − 1 0 1 2 3 4 5 6 7

− 8

− 6

− 4

− 2

0

2

4

6

8

10

Example of secular function

The solution is sought in the interval ]0, δi+1[ with δi+1 > 0

Let us denote δ = δi+

In ]0, δ[ we have ψ(t) < 0 and φ(t) > 0

Assume we know all the poles of the secular function f and we are able to compute values of ψ and φ and their derivatives

BNS methods

Bunch, Nielsen and Sorensen interpolated ψ to first order by a rational function p/(q − t) and φ by r + s/(δ − t) This is called osculatory interpolation The parameters p, q, r , s are determined by matching the exact values of the function and the first derivative of ψ or φ at some given point ¯t (to the right of the exact solution) where f has a negative value q = ¯t + ψ(¯t)/ψ′(¯t) p = ψ(¯t)^2 /ψ′(¯t) r = φ(¯t) − (δ − ¯t)φ′(¯t) s = (δ − ¯t)^2 φ′(¯t)

− (^10) − 2 − 1 0 1 2 3 4 5

− 8

− 6

− 4

− 2

0

2

4

6

8

10

Functions ψ (solid) and φ (dots) as function of t

To obtain an efficient method it remains to find a good initial guess

Melman’s strategy is to look for a zero of an interpolant of the function tf (t) When seeking the ith root, this function is written as

tf (t) = t

c i^2 t

c i^2 + δ − t

• h(t)

The function h is written in two parts h = h 1 + h 2

h 1 (t) =

∑^ i−^1

j=

c^2 j δj − t

, h 2 (t) =

∑^ n

j=i+

c j^2 δj − t