. 107
( 136 .)


uxx + K u = f (x)

where K is constant, the result (after a little algebra) is almost a tridiagonal matrix (Chap-
ter 15, Sec. 10). (Actually, two tridiagonal matrices since we would normally split f (x) into
symmetric and antisymmetric parts (Chapter 8) because the alternative is solving a pen-
tadiagonal matrix of twice the dimension and half its elements equal to zero.) Although
most rows of one of these two Galerkin matrices do contain only three non-zero elements,
we must use one row of each matrix to impose the boundary condition, and this row has
all non-zero elements.
There are two ways to deal with this. First, one can use the formulas for the general LU
decomposition to show that the property of having just one full row can be preserved by
the LU factorization. Consequently, one can apply (B.5)“(B.13) in such a way that all the
2 At least for geophysical problems with so-called “inertial” or “critical” latitudes.

sums are taken over just two elements except for one row where the sum runs over all N
The second way is “matrix-bordering” (Faddeev & Faddeeva, 1963, Strang & Fix, 1973
). To solve

A x = f,

write A in the block form

A11 A12
A21 A22

where the Aij are themselves matrices. The blocks can be of different sizes provided they
follow the pattern shown below where m is an integer such that 1 ¤ m < N

A11 : m — m A12 : m — (N ’ m)
A21 : (N ’ m) — m A22 : (N ’ m) — (N ’ m)

For example, if A is tridiagonal except for a full row at the bottom, we can pick A11 to be
the (N ’ 1) — (N ’ 1) tridiagonal matrix. Then the dimensions of the blocks are

A11 : (N ’ 1) — (N ’ 1) A12 : (N ’ 1) — 1
A21 : 1 — (N ’ 1) A22 : 1 — 1 (scalar!)

A can be written as the product of matrices L and U of block form, too:

L11 0 U11 U12 A11 A12
L21 L22 0 U22 A21 A22

where both L & U are triangular and

L11 U11 = A11

L21 U11 = A21

L11 U12 = A12

L22 U22 = A22 ’ L21 U12 (B.36)

Optionally, one can force the diagonal blocks of U to be identity matrices, but this would
save no labor. To solve (B.35), we must compute the LU factorization of A11 anyway. Con-
sequently, the program is simpler and no less ef¬cient if we choose U11 to be one of the
usual LU factors of A11 . Taking the transpose of (B.34) gives
U11 L21 = A21

which can be easily solved for L21 since the transpose of U11 is a lower triangular matrix.
We then can solve (B.33)“(B.36) in sequence to obtain the block factorization of A.
The crucial point is that if A11 is banded, then L11 and U11 will be banded also. Thus, if
A11 is tridiagonal, as in Orszag and Kells(1980), we can ¬nd its factors in O(N ) operations.
The remaining steps cost O(N ) operations for the same reason. The fact that A21 is a full
row implies that L21 will be full also, but the extra expense is only O(N ).
L. M. Delves has developed spectral iterations which require inverting matrices with an
m — m full block for A11 but with zeros in all elements outside this block except along the
main diagonal (Delves and Freeman, 1981). The work is O((2/3)m3 +kN log2 N ) operations
per iteration; if m N and the number of iterations is not too large, this is much faster
than solving the full N — N matrix problem via Gaussian elimination. At present, Delves™
spectral/block iteration method is the fastest for solving integral equations via Chebyshev
polynomials. It can be applied to differential equations, too. Details are given in the mono-
graph by Delves and Freeman (1981) and Chapter 15.

B.6 Cyclic Banded Matrices (Periodic Boundary Conditions)
For periodic boundary conditions, ¬nite difference methods ” useful as pseudospectral
pre-conditioners ” generate so-called “cyclic” banded matrices. Outside the diagonal
band, the only non-zero elements are those in a small triangular region in the upper right
corner. Temperton (1975) and Navon (1987) describe matrix factorizations for cyclic ma-
trices. However, these are unnecessary for periodic problems with parity: the symmetries
give rise to ordinary banded banded matrices ” very helpful when using ¬nite differences
to precondition an iteration with a Fourier cosine or Chebyshev pseudospectral

B.7 Parting shots
• Always use a direct, rather than an iterative method, if you can afford it. Direct
methods almost always work; iterative schemes sometimes fail.

• Iterative methods are considerably less robust for real-world problems than for the
idealized problems in the textbooks of mathematicians.

• Hotellier proved in 1944 that Gaussian elimination (LU factorization) could magnify
roundoff error by 4N where N is the size of the matrix. This would seem to make it
impossible to solve matrix problems with N > 40 even in seventeen decimal place
arithmetic! However, this worst case is exceedingly rare in practice. In numerical
analysis, there is sometimes an enormous gap between what can be proven and what
actually happens.

• Clever matrix-of-a-matrix methods have lowered the asymptotic cost of solving a
dense matrix from O(N 3 ) operations to O(N 2.376 ). L. N. Trefethen has a $100 bet with
P. Alfeld that the exponent will lowered to O(N 2+ ) for arbitrarily small by 2006.
However, these new algorithms have huge proportionality constants. Consequently,
Gaussian elimination is the best unless N is really, really huge. A numerical analyst™s
concept of ef¬ciency is the N ’ ∞ limit, which may have little to do with engineering

Trefethen and Bau(1997) give more information about roundoff and asymptotically
cheap algorithms.
Karageorghis and Paprzycki(1996, 1998), and earlier papers by Karageorghis and col-
laborators, describe ef¬cient direct methods even for the somewhat complicated matrices
generated by spectral domain decomposition.
Appendix C

The Newton-Kantorovich Method
for Nonlinear Boundary &
Eigenvalue Problems

“Nature and Nature™s laws were hid in night. God said: Let Newton be! and all was light.”
” Alexander Pope (1688“1744)

“ . . . with his prism and silent face . . . a mind forever voyaging through strange seas
of thought, alone.”
” William Wordsworth (of Sir Isaac Newton)

C.1 Introduction
Most algorithms for solving nonlinear equations or systems of equations are variations of
Newton™s method. When we apply spectral, ¬nite difference, or ¬nite element approxima-
tions to a nonlinear differential equation, the result is a nonlinear system of algebraic equations
equal in number to the total number of degrees of freedom, regardless of whether the de-
grees of freedom are the grid point values of the solution, {u(xi )}, or the coef¬cients of
the spectral series for u(x), {an }. We must then apply Newton™s method ” or something
equivalent to it ” to iterate from a ¬rst guess until the numerical solution is suf¬ciently
There are several practical dif¬culties. First, we need a ¬rst guess ” this is a universal
problem, and no choice of iteration scheme can eliminate it. Second, Newton™s method
requires the elements of the “Jacobian matrix” of the system of equations. If we write the
system as

Fi (a0 , a1 , . . . , aN ) = 0 i = 0, . . . , N

then the elements of the Jacobian matrix J are

Jij ≡ ‚Fi /‚aj (C.2)
i, j = 0, . . . , N


It is usually messy and time-consuming to calculate the Jacobian, which must (in a strict
Newton™s method) be recomputed at each iteration. Third, each iteration requires inverting
the Jacobian. If we let a(i) denote the vector of the spectral coef¬cients at the i-th iteration,
a(i+1) = a(i) ’ J (C.3)

where F is the vector whose elements are the equations of the system and where both F
and J are evaluated using a = a(i) . This may be the most expensive step of all since the
Jacobian matrix ” at least with pseudospectral methods ” is usually a full matrix, so it
costs O(N 3 ) operations to determine each iterate from its predecessor.
Resolving the ¬rst dif¬culty requires the “continuation” method discussed below, but
we can often reduce the magnitude of the second and third problems by applying New-
ton™s method directly to the original differential equation, and then applying the spectral
method to convert the resulting iterative sequence of linear differential equations into the
corresponding linear matrix problem. When Newton™s method is applied directly to a
differential equation ” rather than to algebraic equations ” it is often called the “Newton-
Kantorovich” method.
The justi¬cation for the ordinary Newton™s method is a Taylor expansion followed by
linearization. If we wish to solve f (x) = 0, then assuming we have a good guess for a root,
x ≈ x(i) , we can Taylor expand f (x) to obtain
f (x) = f (x(i) ) + fx (x(i) ) x ’ x(i) + O x ’ x(i) (C.4)

If our guess is suf¬ciently good, we can ignore the quadratic terms and treat (C.4) as a linear
equation for x,

f (x(i) )
x(i+1) = x(i) ’ (C.5)
fx x(i)

This is identical with (C.3) except for the extension from one unknown to many.
To illustrate the generalization known as the Newton-Kantorovich method, consider
the ¬rst order differential equation

ux = F (x, u[x])

Suppose we have an iterate, u(i) (x), which is a good approximation to the true solution
u(x). Taylor expand (C.6) about this function:
u(x) ’ u(i) (x) + O u ’ u(i)
ux = F x, u(i) [x] + Fu x, u(i) [x] (C.7)

where Fu denotes the derivative of the function F (x, u) with respect to u. This implies that
the linear differential equation for the next iterate is

’ {Fu }u(i+1) = F ’ Fu u(i)

where F and Fu are both evaluated with u = u(i) . Equivalently, we can solve a differential
equation for the correction to the iterate by de¬ning

u(i+1) (x) ≡ u(i) (x) + (C.9)

Then (x) is the solution of

’ {Fu } = F ’ ux

Since Taylor expansions about a function instead of a point may be unfamiliar, a brief ex-
planation is in order. Technically, operations on functionals, which we may loosely de¬ne
as “functions of a function”, are described by a branch of mathematics called “functional
analysis”. The most elementary and readable introduction is a little paperback by W. W.
Sawyer (1978).
The punchline is this: the linear term in a generalized Taylor series is actually a linear op-
erator. In (C.3), the operator is multiplication by the Jacobian matrix. In (C.10), the operator
is the differential operator that is applied to (x) to generate the L. H. S. of the equation. In
either case, this linear operator is known as a “Frechet derivative”. The general de¬nition
is the following.

De¬nition 49 (FRECHET DERIVATIVE) Given a general nonlinear operator N (u), which may
include derivatives or partial derivatives, its FRECHET DIFFERENTIAL is de¬ned by

N (u + (x)) ’ N (u)
≡ lim
Nu (C.11)

which is equivalent to

‚N (u + )
Nu ≡ (C.12)
‚ =0

The numerical result, Nu , is the FRECHET DIFFERENTIAL of the operator N in the direction
while Nu is the FRECHET DERIV ATIVE and is a linear operator.


. 107
( 136 .)