- •1 Introduction
- •1.1 What makes eigenvalues interesting?
- •1.2 Example 1: The vibrating string
- •1.2.1 Problem setting
- •1.2.2 The method of separation of variables
- •1.3.3 Global functions
- •1.3.4 A numerical comparison
- •1.4 Example 2: The heat equation
- •1.5 Example 3: The wave equation
- •1.6 The 2D Laplace eigenvalue problem
- •1.6.3 A numerical example
- •1.7 Cavity resonances in particle accelerators
- •1.8 Spectral clustering
- •1.8.1 The graph Laplacian
- •1.8.2 Spectral clustering
- •1.8.3 Normalized graph Laplacians
- •1.9 Other sources of eigenvalue problems
- •Bibliography
- •2 Basics
- •2.1 Notation
- •2.2 Statement of the problem
- •2.3 Similarity transformations
- •2.4 Schur decomposition
- •2.5 The real Schur decomposition
- •2.6 Normal matrices
- •2.7 Hermitian matrices
- •2.8 Cholesky factorization
- •2.9 The singular value decomposition (SVD)
- •2.10 Projections
- •2.11 Angles between vectors and subspaces
- •Bibliography
- •3 The QR Algorithm
- •3.1 The basic QR algorithm
- •3.1.1 Numerical experiments
- •3.2 The Hessenberg QR algorithm
- •3.2.1 A numerical experiment
- •3.2.2 Complexity
- •3.3 The Householder reduction to Hessenberg form
- •3.3.2 Reduction to Hessenberg form
- •3.4 Improving the convergence of the QR algorithm
- •3.4.1 A numerical example
- •3.4.2 QR algorithm with shifts
- •3.4.3 A numerical example
- •3.5 The double shift QR algorithm
- •3.5.1 A numerical example
- •3.5.2 The complexity
- •3.6 The symmetric tridiagonal QR algorithm
- •3.6.1 Reduction to tridiagonal form
- •3.6.2 The tridiagonal QR algorithm
- •3.7 Research
- •3.8 Summary
- •Bibliography
- •4.1 The divide and conquer idea
- •4.2 Partitioning the tridiagonal matrix
- •4.3 Solving the small systems
- •4.4 Deflation
- •4.4.1 Numerical examples
- •4.6 Solving the secular equation
- •4.7 A first algorithm
- •4.7.1 A numerical example
- •4.8 The algorithm of Gu and Eisenstat
- •4.8.1 A numerical example [continued]
- •Bibliography
- •5 LAPACK and the BLAS
- •5.1 LAPACK
- •5.2 BLAS
- •5.2.1 Typical performance numbers for the BLAS
- •5.3 Blocking
- •5.4 LAPACK solvers for the symmetric eigenproblems
- •5.6 An example of a LAPACK routines
- •Bibliography
- •6 Vector iteration (power method)
- •6.1 Simple vector iteration
- •6.2 Convergence analysis
- •6.3 A numerical example
- •6.4 The symmetric case
- •6.5 Inverse vector iteration
- •6.6 The generalized eigenvalue problem
- •6.7 Computing higher eigenvalues
- •6.8 Rayleigh quotient iteration
- •6.8.1 A numerical example
- •Bibliography
- •7 Simultaneous vector or subspace iterations
- •7.1 Basic subspace iteration
- •7.2 Convergence of basic subspace iteration
- •7.3 Accelerating subspace iteration
- •7.4 Relation between subspace iteration and QR algorithm
- •7.5 Addendum
- •Bibliography
- •8 Krylov subspaces
- •8.1 Introduction
- •8.3 Polynomial representation of Krylov subspaces
- •8.4 Error bounds of Saad
- •Bibliography
- •9 Arnoldi and Lanczos algorithms
- •9.2 Arnoldi algorithm with explicit restarts
- •9.3 The Lanczos basis
- •9.4 The Lanczos process as an iterative method
- •9.6 Partial reorthogonalization
- •9.5 An error analysis of the unmodified Lanczos algorithm
- •9.7 Block Lanczos
- •9.8 External selective reorthogonalization
- •Bibliography
- •10 Restarting Arnoldi and Lanczos algorithms
- •10.2 Implicit restart
- •10.3 Convergence criterion
- •10.4 The generalized eigenvalue problem
- •10.5 A numerical example
- •10.6 Another numerical example
- •10.7 The Lanczos algorithm with thick restarts
- •10.8 Krylov–Schur algorithm
- •10.9 The rational Krylov space method
- •Bibliography
- •11 The Jacobi-Davidson Method
- •11.1 The Davidson algorithm
- •11.2 The Jacobi orthogonal component correction
- •11.2.1 Restarts
- •11.2.2 The computation of several eigenvalues
- •11.2.3 Spectral shifts
- •11.3 The generalized Hermitian eigenvalue problem
- •11.4 A numerical example
- •11.6 Harmonic Ritz values and vectors
- •11.7 Refined Ritz vectors
- •11.8 The generalized Schur decomposition
- •11.9.1 Restart
- •11.9.3 Algorithm
- •Bibliography
- •12 Rayleigh quotient and trace minimization
- •12.1 Introduction
- •12.2 The method of steepest descent
- •12.3 The conjugate gradient algorithm
- •12.4 Locally optimal PCG (LOPCG)
- •12.5 The block Rayleigh quotient minimization algorithm (BRQMIN)
- •12.7 A numerical example
- •12.8 Trace minimization
- •Bibliography
2 |
CHAPTER 1. INTRODUCTION |
1.2Example 1: The vibrating string
1.2.1Problem setting
Let us consider a string as displayed in Fig. 1.1. The string is fixed at both ends, at x = 0
u |
|
|
u(x,t) |
|
x |
0 |
L |
Figure 1.1: A vibrating string fixed at both ends.
and x = L. The x-axis coincides with the string’s equilibrium position. The displacement of the rest position at x, 0 < x < L, and time t is denoted by u(x, t).
We will assume that the spatial derivatives of u are not very large:
|
∂u |
|
|
|
|
|
|
|
∂x |
|
is small. |
|
|
|
|
This assumption entails that we may neglect terms of higher order.
Let v(x, t) be the velocity of the string at position x and at time t. Then the kinetic energy of a string section ds of mass dm = ρ ds is given by
|
|
|
|
|
|
|
1 |
|
|
|
1 |
|
|
|
∂u |
2 |
||
(1.1) |
|
|
dT = |
|
dm v2 = |
|
|
ρ ds |
|
|
. |
|||||||
2 |
2 |
∂t |
||||||||||||||||
From Fig. 1.2 we see that ds2 = dx2 + ∂u∂x 2 dx2 and thus |
|
|||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
ds |
= s1 + |
∂u |
|
2 |
1 |
|
|
∂u |
2 |
|
|
|
|||||
|
|
|
|
|
= 1 + |
|
|
|
|
|
+ higher order terms. |
|||||||
|
dx |
∂x |
|
2 |
∂x |
|||||||||||||
Plugging this into (1.1) and omitting also the second order term (leaving just the number 1)
gives
dT = ρ dx ∂u 2 . 2 ∂t
The kinetic energy of the whole string is obtained by integrating over its length,
T = |
Z0 |
dT (x) = 2 |
Z0 |
ρ(x) |
∂t |
dx |
|
L |
1 |
L |
|
∂u |
2 |
The potential energy of the string has two components
1.2. EXAMPLE 1: THE VIBRATING STRING |
3 |
ds |
dx |
Figure 1.2: A vibrating string, local picture.
1. the stretching times the exerted strain τ, |
|
|
|
|
|
|
|||||||||||
τ Z0 |
ds − τ Z0 |
dx = τ |
Z0 |
s |
|
|
|
|
− 1 dx |
||||||||
1 + ∂x |
2 |
||||||||||||||||
L |
L |
|
|
L |
1 |
|
|
∂u |
|
||||||||
|
|
|
|
L |
∂u |
|
2 |
|
|||||||||
|
|
= τ Z0 |
|
|
|
|
|
|
+ |
higher order terms! dx |
|||||||
|
|
2 |
∂x |
||||||||||||||
2. exterior forces of density f, |
|
|
− Z0L fudx. |
|
|
||||||||||||
|
|
|
|
|
|
|
|||||||||||
Summing up, the potential energy of the string becomes |
|
|
|||||||||||||||
|
|
L |
|
τ |
∂u |
|
2 |
|
|
|
|
||||||
(1.2) |
V = Z0 |
|
|
|
|
|
|
− fu! dx. |
|||||||||
2 |
|
∂x |
|||||||||||||||
To consider the motion (vibration) of the string in a certain time interval t1 ≤ t ≤ t2 we form the integral
(1.3) |
I(u) = |
Zt1 |
(T − V ) dt = |
2 |
Zt1 |
Z0 |
"ρ(x) |
∂t |
− τ |
∂x |
+ 2fu# dx dt |
|
|
t2 |
|
1 |
t2 |
|
L |
∂u |
2 |
∂u |
2 |
Here functions u(x, t) are admitted that are di erentiable with respect to x and t and satisfy the boundary conditions (BC) that correspond to the fixing,
(1.4) |
u(0, t) = u(L, t) = 0, |
t1 ≤ t ≤ t2, |
||
as well as given initial conditions and end conditions, |
||||
(1.5) |
u(x, t1) = u1 |
(x), |
0 < x < L. |
|
u(x, t2) = u2 |
(x), |
|||
|
|
|||
4 |
CHAPTER 1. INTRODUCTION |
According to the principle of Hamilton a mechanical system with kinetic energy T and potential energy V behaves in a time interval t1 ≤ t ≤ t2 for given initial and end positions
such that
Z t2
I = L dt, L = T − V,
t1
is minimized.
Let u(x, t) be such that I(u) ≤ I(w) for all w, that satisfy the initial, end, and boundary conditions. Let w = u + ε v with
v(0, t) = v(L, t) = 0, v(x, t1) = v(x, t2) = 0.
v is called a variation. We now consider I(u + ε v) as a function of ε. Then we have the equivalence
|
|
|
|
|
|
|
|
dI |
|
|
|
|
|
|
|
|||||||
I(u) minimal |
|
|
dε (u) = 0 for all admitted v. |
|
||||||||||||||||||
Plugging u + ε v into eq. (1.3) we obtain |
|
|
|
− τ ∂(u∂x |
+ 2f(u + ε v)# dx dt |
|||||||||||||||||
I(u + ε v) = 2 Z |
Z "ρ(x) |
|
|
∂t |
|
|
|
|||||||||||||||
1 |
|
t2 |
L |
|
∂(u + ε v) 2 |
|
|
+ ε v) 2 |
||||||||||||||
(1.6) |
|
|
t1 |
0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Z ρ(x) ∂t ∂t − τ ∂x ∂x + 2fv dx dt + O(ε2). |
||||||||||||||||||
= I(u) + ε Z |
||||||||||||||||||||||
|
|
|
|
t2 |
L |
|
|
|
∂u ∂v |
|
|
∂u ∂v |
|
|
||||||||
|
|
|
|
t1 |
0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Thus, |
= Zt1 |
Z0 |
ρ ∂t2 |
|
|
− τ ∂x2 + 2 f v dx dt = 0 |
||||||||||||||||
|
|
|
∂ε |
|
|
|||||||||||||||||
|
|
|
∂I |
t2 |
L |
|
|
∂2u |
|
|
|
∂2u |
|
|
||||||||
for all admissible v. Therefore, the bracketed expression must vanish, |
||||||||||||||||||||||
|
|
|
|
|
|
∂2u |
|
|
|
|
∂2u |
|
|
|
|
|
||||||
(1.7) |
|
|
|
|
− ρ |
|
|
|
+ τ |
|
|
|
= 2 f. |
|
|
|||||||
|
|
|
|
∂t2 |
|
∂x2 |
|
|
|
|||||||||||||
This last di erential equation is named Euler-Lagrange equation.
If the force is proportional to the displacement u(x, t) (like, e.g., in Hooke’s law) then we get a di erential equation of the form
−ρ(x)∂2u + ∂ p(x)∂u + q(x)u(x, t) = 0. (1.8) ∂t2 ∂x ∂x
u(0, t) = u(1, t) = 0
which is a special case of the Euler-Lagrange equation (1.7). Here, ρ(x) plays the role of a mass density, p(x) of a locally varying elasticity module. We do not specify initial and end conditions for the moment. Note that there are no external forces present in (1.7).
From physics we know that ρ(x) > 0 and p(x) > 0 for all x. These properties are of importance also from a mathematical view point! For simplicity, we assume that ρ(x) = 1.
1.2. EXAMPLE 1: THE VIBRATING STRING |
5 |
1.2.2The method of separation of variables
For the solution u in (1.8) we make the ansatz
(1.9) u(x, t) = v(t)w(x).
Here, v is a function that depends only on the time t, while w depends only on the spacial variable x. With this ansatz (1.8) becomes
(1.10) v′′(t)w(x) − v(t)(p(x)w′(x))′ − q(x)v(t)w(x) = 0.
Now we separate the variables depending on t from those depending on x,
v′′(t) |
= |
|
1 |
(p(x)w′(x))′ + q(x). |
|
v(t) |
w(x) |
||||
|
|
||||
This equation holds for any t and x. We can vary t and x independently of each other without changing the value on each side of the equation. Therefore, each side of the equation must be equal to a constant value. We denote this value by −λ. Thus, from the left side we obtain the equation
(1.11) |
− v′′(t) = λv(t). |
√ |
|
√ |
|
|
|
This equation has the well-known solution v(t) = a · cos( |
λt) + b · sin( λt) where λ > 0 |
||||||
is assumed. The right side of (1.10) gives a so-called Sturm-Liouville problem |
|||||||
|
|
|
|
||||
(1.12) |
−(p(x)w′(x))′ + q(x)w(x) = λw(x), |
w(0) = w(1) = 0. |
|
||||
A value λ for which (1.12) has a non-trivial (i.e. nonzero) solution w is called an eigenvalue; w is a corresponding eigenfunction. It is known that all eigenvalues of (1.12) are
positive. By means of our ansatz (1.9) we get
h √ √ i u(x, t) = w(x) a · cos( λt) + b · sin( λt)
as a solution of (1.8). It is known that (1.12) has infinitely many real positive eigenvalues
0 < λ1 ≤ λ2 |
≤ · · · , (λkk→∞∞). (1.12) has a non-zero solution, say wk(x), only for these |
||||||
|
−→ |
|
|
|
|
|
|
particular values λk. Therefore, the general solution of (1.8) has the form |
|||||||
(1.13) |
u(x, t) = k=0 wk(x) hak · cos( |
|
λk t) + bk · sin( |
λk t)i . |
|||
|
∞ |
|
|
|
|
|
|
|
X |
p |
|
|
p |
|
|
The coe cients ak and bk require that
are determined by initial and end conditions. We could, e.g.,
u(x, 0) =
∂u∂t (x, 0) =
where u0 and u1 are given functions. It is known that the wk form an orthogonal basis in the space of square integrable functions L2(0, 1). Therefore, it is not di cult to compute the coe cients ak and bk.
In concluding, we see that the di cult problem to solve is the eigenvalue problem (1.12). Knowing the eigenvalues and eigenfunctions the general solution of the time-dependent problem (1.8) is easy to form.
Eq. (1.12) can be solved analytically only in very special situation, e.g., if all coe cients are constants. In general a numerical method is needed to solve the Sturm-Liouville problem (1.12).
6 |
CHAPTER 1. INTRODUCTION |
1.3Numerical methods for solving 1-dimensional problems
In this section we consider three methods to solve the Sturm-Liouville problem.
1.3.1Finite di erences
We approximate w(x) by its values at the discrete points xi = ih, h = 1/(n + 1), i = 1, . . . , n.
|
|
|
|
|
|
|
x |
|
|
|
|
|
|
|
|
0 |
xi−1 xi |
x |
i+1 |
L |
|
||
Figure 1.3: Grid points in the interval (0, L).
At point xi we approximate the derivatives by finite di erences. We proceed as
follows. First we write |
|
|
|
|
|
|
|
|
g(xi+ 21 ) − g(xi− 12 ) |
|
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
|
|
d |
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
g(xi) ≈ |
|
|
|
|
|
|
. |
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
dx |
|
|
|
|
h |
|
|
|
|
|
|
||||||
For g = pdw |
we get |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
dx |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
g(x |
i+ |
1 ) = p(x |
i+ |
1 ) |
w(xi+1) − w(xi) |
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
2 |
|
|
2 |
|
h |
|
|
|
|
|
|
|
|||
and finally, for i = 1, . . . , n, |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
d |
dw |
1 |
|
|
|
|
|
w(xi+1) |
w(xi) |
|
|
w(xi) |
w(xi |
|
1) |
|
||||||||
− |
|
p |
|
(xi) |
≈ − |
|
p(xi+ 21 ) |
|
|
|
|
− |
|
− p(xi− 21 ) |
|
− |
− |
|
|
||||||
dx |
dx |
h |
|
|
|
|
h |
|
|
h |
|||||||||||||||
= h2 hp(xi− |
21 )wi−1 + (p(xi− 12 ) + p(xi+ 21 ))wi − p(xi+ 21 )wi+1i . |
|
1 |
|
|
Note that at the interval endpoints w0 = wn+1 = 0.
We can collect all equations in a matrix equation,
|
p(x1 ) + p(x3 ) |
||||
2 |
2 |
|
+ q(x1) |
||
|
h2 |
||||
|
|
3 |
|
|
|
|
|
p(x2 ) |
|
||
|
2 |
||||
|
|
− h |
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
or, briefly,
|
p(x3 ) |
|||
− h2 |
||||
|
2 |
|
|
|
p(x3 ) + p(x5 ) |
||||
2 |
2 |
+ q(x2) |
||
h2 |
|
|||
p(x5 ) |
||||
− |
||||
2 |
|
|
||
h2 |
|
|||
p(x5 )
− h22
...
|
|
w2 |
|
|
w2 |
|
|
|||
|
|
|
|
w1 |
|
|
|
w1 |
|
|
|
|
|
. |
|
. |
|
||||
|
|
|
. |
|
|
|
. |
|
|
|
.. |
|
|
. |
|
. |
|
||||
|
|
|
w3 |
|
= λ |
|
w3 |
|
, |
|
|
. |
|
|
|
|
|
||||
|
|
wn |
|
wn |
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(1.14) |
Aw = λw. |
By construction, A is symmetric and tridiagonal. One can show that it is positive definite as well. Note that this matrix has just a few nonzeros: out of the n2 elements of A only 3n − 2 are nonzero. This is a first example of a sparse matrix.
1.3.2The finite element method
We write (1.12) in the form
1.3. NUMERICAL METHODS FOR SOLVING 1-DIMENSIONAL PROBLEMS |
7 |
Find a twice di erentiable function w with w(0) = w(1) = 0 such that
Z 1
−(p(x)w′(x))′ + q(x)w(x) − λw(x) φ(x)dx = 0
0
for all smooth functions φ that satisfy φ(0) = φ(1) = 0.
To relax the requirements on w we integrate by parts and get the new so-called weak form of the problem:
Find a di erentiable function w with w(0) = w(1) = 0 such that
Z 1
(1.15) p(x)w(x)′φ′(x) + q(x)w(x)φ(x) − λw(x)φ(x) dx = 0
0
for all di erentiable functions φ that satisfy φ(0) = φ(1) = 0.
Remark: Requiring di erentiability is too strong and does not lead to a mathematically suitable formulation. In particular, the test functions that will be used below are not di erentiable in the classical sense. It is more appropriate to require w and φ to be weakly di erentiable. In terms of Sobolev spaces: w, φ H01([0, 1]). An introduction to Sobolev spaces is, however, beyond the scope of these notes.
|
Ψi |
|
|
|
||||
|
1 |
|
|
|
|
|
||
|
|
|
|
|
|
|
|
x |
|
|
|
|
|
|
|
|
|
0 |
|
xi−1 xi |
x |
i+1 |
L |
|
||
Figure 1.4: A basis function of the finite element space: a hat function. We now write w as the linear combination
|
|
|
|
|
|
n |
|
|
|
|
|
(1.16) |
|
|
|
|
w(x) = |
ξi Ψi(x), |
|
|
|
||
|
|
|
|
|
|
=1 |
|
|
|
|
|
|
|
|
|
|
|
Xi |
|
|
|
|
|
where |
|
|
|
|
|
| + |
|
|
|
|
|
(1.17) |
Ψ |
(x) = |
1 |
− | |
x − xi |
= max |
0, 1 |
− |
|x − xi| |
, |
|
|
i |
|
|
h |
{ |
|
h |
} |
|||
is the function that is linear in each interval (xi, xi+1) and satisfies |
|
||||||||||
|
|
|
|
|
|
|
1, |
i = k, |
|
|
|
|
|
|
Ψi(xk) = δik := 0, |
i 6= k. |
|
|
|||||
An example of such a basis function, a so-called hat function, is given in Fig. 1.4.
We now replace w in (1.15) by the linear combination (1.16), and replace testing ‘against all φ’ by testing against all Ψj . In this way (1.15) becomes
0 |
1 |
−p(x)( |
n |
ξi Ψi′ (x))Ψj′ (x) + (q(x) − λ) |
n |
ξi Ψi(x)Ψj (x)! dx, for all j, |
Z |
|
|
X |
|
X |
|
|
|
i=1 |
|
i=1 |
|
