. 33
( 136 .)


usual Chebyshev-Lobatto grid:
N ’1
π(x) ≡ (x ’ xi ) (7.37)

After the usual spectral discretization using these basis functions to give the generalized
matrix eigenproblem Au = » Bu, Huang and Sloan multiply both matrices by the diagonal
matrices whose elements are

Djj = (1 ’ x2 )2 , j = 2, . . . , N ’ 1 (7.38)

Like Heinrichs™ improvement for the non-cardinal basis, Huang and Sloan™s method
reduces the condition number to its square root. For a differential equation of order p, the
condition number is reduced from O(N 2p ) to O(N p ). Thus, it is possible to solve differen-
tial equations of quite high order with large N .

Although we have chosen to discuss condition number in the midst of a chapter on
eigenproblems because condition number is measured as a ratio of eigenvalues, the Hein-
richs and Huang and Sloan basis sets can be equally well employed for high order BOUND-
ARY VALUE problems.

7.8 Alternatives for Costly Problems, I: The Power Method
When the QR/QZ algorithm is too expensive, one needs iteration schemes which compute
only a single eigenvalue at a time. The simplest is the “power method”, which can be
applied both to the matrix discretization of an eigenproblem and also to the original time-
dependent equations from which the eigenproblem was derived. A couple of examples
will make the idea clearer.
First, suppose the goal is to calculate hydrodynamic instabilities. The ¬‚ow is split into
a user-speci¬ed “basic state” plus a perturbation, which is assumed to be initially small
compared to the basic state but is otherwise unrestricted. The ¬‚ow is unstable if the per-
turbation ampli¬es in time.
A useful illustration is the barotropic vorticity equation in a “beta-plane” approxima-
tion on a rotating planet:

ψxxt + ψyyt + ψx (ψxxy + ψyyy ) ’ ψy (ψxxx + ψxyy ) + βψx = 0 (7.39)

where ψ(x, y, t) is the streamfunction for two-dimensional ¬‚ow and β is the latitudinal
derivative of the Coriolis parameter. For simplicity, we shall consider only a basic state
which is independent of x, that is, denoting the basic state by upper case letters,
Ψ(y) = ’ (7.40)
U (z)dz

where U (y) is the “mean ¬‚ow”, a jet moving parallel to the x-axis.
The time-dependent variant of the power method is to integrate the equations of motion
from an arbitrary initial condition until the fastest-growing mode completely dominates the
solution, at which point the solution is the desired eigenmode and its growth rate and
spatial translation rate are the imaginary and real parts of the complex-valued phase speed
c, which is the eigenvalue.
There are two variants. The ¬rst is to linearize the equations of motion about the basic
ψ ∼ O( ), (7.41)
ψ=ψ + U (z)dz, << 1

Since the perturbation is assumed to be very small, we can always expand the equations
of motion in a power series in the amplitude of the perturbation; the linearized equation
of motion is simply the O( ) term in this expansion. The linearized barotropic vorticity
equation, for example, is

ψxxt + (β ’ Uyy )ψx + U (y)(ψxxx + ψxyy ) = 0 (7.42)

where the prime on ψ denotes the perturbative part of the streamfunction.
In the linearized equation, feedback from the perturbation to the basic state, as happens
for the fully nonlinear equations of motion, is suppressed. This is a detriment because it
restricts the model to the early stages of the instability when the perturbation is very small.
However, it is also a virtue. The linearized stability problem is completely speci¬ed by the
speci¬cation of the basic state.

The second variant of the time-dependent power method is to solve the full, nonlinear
equations of motion by starting from a perturbation of amplitude and integrating until
the amplitude is O(a) where a >> , but a is still small compared to the basic state. The
result will be the same as that of solving the linearized problem (on the same time interval)
to within an absolute error of O(a2 ) [relative error O(a)] because this is the leading order of
the terms neglected in the linearization. One virtue of this method is that it is unnecessary
to create a second computer program, with slightly different (i. e., linearized) equations of
motion, to explore linearized stability. Another is that one can easily extend the integration
further in time to explore how nonlinear feedback alters the growth and structure of the
amplifying disturbances. The disadvantage is that if there are unstable modes of similar
growth rates, one may be forced to stop the integration (to keep a << 1) before the fastest
growing mode has become large compared to all other modes.
In either variant, the growth rate and amplitude are estimated by ¬tting the approxi-
ψ(x, y, t2 ) ’ ψ(x, y, t1 ) ≈ Ξ(y) exp(ikx) {’ikct2 ) ’ exp(’ikct1 )}
+complex conjugate (7.43)
where t1 and t2 are two large but otherwise arbitrary times and Ξ(y) gives the complex-
valued latitudinal structure of the instability while c is the complex phase speed and k the
x-wavenumber, which must also be determined from the perturbation for large time.
Because the ¬tting is a little tricky, and it sometimes requires a huge number of time
steps before the fastest-growing mode dominates its slower-growing brethren, it is often
more convenient to convert the time-dependent equations into a matrix eigenproblem. The
justi¬cation for this step is that as long as the linearized equations have coef¬cients which
are independent of time, all discretizations of the spatial dependence will collapse the lin-
earized PDE into a system of ordinary differential equations in time of the form du/dt = Au
where A is a square matrix whose elements are independent of time. The usual theory of
constant coef¬cient ODEs, found in most undergraduate differential equation texts, then
asserts that the general solution is a superposition of the eigenmodes of A, each oscillating
as an exponential of time with a complex frequency which is the corresponding eigenvalue
of A, assuming these are all distinct. The matrix formulation of linearized stability theory
is simply the computation of these eigenvectors and eigenvalues of A.
When the coef¬cients of the linearized PDE do not vary with all the spatial coef¬cients,
then one can separate variables still further. For the barotropic vorticity equation, one can
expand the solution in a Fourier series in x; if linearized with respect to a mean current
which varies only with y, the x-wavenumbers are uncoupled. A large two-dimensional
eigenvalue problem therefore collapses into a set of differential eigenproblems in y only,
one for each wavenumber in x. Denoting the wavenumber by k, the eigenproblem is
β ’ Uyy
’ k2 ψ = 0 (7.44)
ψyy +
U ’c
(In a minor abuse of notation that is almost ubiquitous in the literature, ψ(y) is used for
the eigenfunction even though the same symbol has already been employed for the total
time-dependent streamfunction.)
The power method can be applied directly to the matrix form of an eigenvalue problem
by ¬rst discretizing the differential operator and then multiplying through by the appro-
priate inverse matrix to convert the problem into a standard (as opposed to a generalized)
matrix eigenproblem:

Au = »u

Note that if we assume that A has as complete set of eigenvectors ej and expand a vector

u0 = aj ej

then the result of multiplying u0 by the square matrix A will be a vector whose eigencoef-
¬cients will be the same except that aj is multiplied by »j , the eigenvalue associated with
ej :

Au0 = »j aj ej

If we iterate by repeatedly multiplying by powers of A, it follows that as the iteration
number k increases, the eigencoef¬cients will be multiplied by »k and therefore the eigen-
function associated with the eigenvalue of largest magnitude will more and more dominate
the k-th power of A. As it does so, the ratio ||Au||/||u|| will tend more and more towards
the magnitude of the largest eigenvalue where || || denotes any reasonable matrix/vector
This suggests the following algorithm for computing the eigenfunction and eigenvalue
for the mode of largest (absolute value) of eigenvalue. First, choose an arbitrary starting
vector u0 and then divide it by its norm so as to rescale it to unit norm. (The algorithm will
fail if u0 is orthogonal to the target eigenmode, but if u0 is a column of numbers from a
random number generator, the probability of such failure is negligibly small.) Then, repeat
the following loop over the iteration number k until the eigenfunction and eigenvalue have
converged to within the user-chosen tolerance:

uk+1 = Auk , k = 0, 1, . . .
= ||uk+1 || (7.48)
= uk+1 / ||uk+1 ||

This is simply the matrix form of the power method: iterate until the fastest-growing mode
dominates. The third line of the loop is meant in the sense of a computational assignment
statement rather than mathematical equality ” uk+1 is replaced by the same vector divided
by its norm to rescale it to unit norm. This renormalization avoids over¬‚ow problems.
When » is complex-valued, the second line must be replaced by the average of the element-
by-element division of the new eigenvector by the old:
»k+1 = (1/N ) uk+1,j /uk,j

The matrix power method can be applied to any matrix, whether connected with sta-
bility problems or differential eigenproblems or not. It is very fast and cheap per iteration
because the most costly step is merely a matrix-vector multiply, which is a highly vectoriz-
able operation. However, the power method has the disadvantage that it can ¬nd only the
mode of largest absolute value of eigenvalue, which may not be the only mode of interest.
In stability problems, it is quite common to have two unstable modes whose growth rates
switch dominance at some point in parameter space. In such transition regions where two
modes have the same or almost the same growth rates, the power method will fail.
In the next section, we therefore describe another algorithm which ¬xes many of the
defects of the power method, albeit at a much higher cost per iteration. First, though, we

First Guess
-30 -20 -10 0 10 20 30

Dispersive Soliton



-30 -20 -10 0 10 20 30

Figure 7.10: Illustration of the power method applied to the KdV equation. The initializa-
tion (“¬rst guess”) is u = 4 exp(’x2 /2). At t = 20, this has split into a dispersing transient
(which moves left) and a single soliton (which moves rightward).

must note that the power method is not limited only to instability problems or only to
linear eigenvalue problems.
For example, the power method is a good way to compute solitary waves of the Korteweg-
deVries (KdV) equation:

ut + uux + uxxx = 0

The solitary waves are the solutions to the nonlinear eigenvalue problem

(u ’ c)uX + uXXX = 0 (7.51)

where X ≡ x ’ ct is the coordinate in a frame of reference which is travelling with the
wave and c is the phase speed. By the power method, the ODE nonlinear eigenproblem
can be bypassed and solitons computed directly from the time-dependent equation: Inte-
grating forward in time from an arbitrary initial condition until the dominant eigenmode
has separated itself from the transients.
In stability problems, the separation is by separation-by-amplitude: the dominant mode
is fastest-growing and therefore eventually becomes larger than the other modes every-
where in space. For the KdV equation, the separation is spatial: the dispersing transients
travel leftward (because the group velocity for small amplitude waves of all wavenumbers
is negative) whereas the velocities of all solitary waves are positive. Because the amplitude
of a solitary wave decreases exponentially fast with distance from the center of the soliton,
the overlap between a given solitary wave and the rest of the solution decreases exponen-
tially with time. Thus, even though there is no instability and the eigenvalue problem is
nonlinear, the power method converges geometrically for the KdV equation.

7.9 Alternatives for Costly Problems, II: Inverse Power Method
The matrix variant of the power method has the disadvantage that it can ¬nd only the mode
of largest eigenvalue. The inverse power method, also sometimes called simply “inverse
iteration”, lacks this limitation. The algorithm is to apply the ordinary power method to
the matrix

M ≡ (A ’ ΛI)’1 (7.52)

where I is the identity matrix and Λ is a user-speci¬ed constant. The eigenvalues of M are

µj =
»j ’ Λ

where the »j are the eigenvalues of the original matrix A, so the largest eigenvalue of the
shifted-and-inverted matrix M is the reciprocal of whichever eigenvalue of A is closest to
the shift Λ. It follows that the inverse power method will converge geometrically to any
eigenvalue of the matrix A, provided we have a suf¬ciently good ¬rst guess (and set Λ
equal to it).
As usual, the inverse is not explicitly computed. Instead, one iterates

(A ’ Λk I)uk+1 = uk , k = 0, 1, . . .


. 33
( 136 .)