(B.27)

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.

B.5. BLOCK AND “BORDERED” MATRICES 523

sums are taken over just two elements except for one row where the sum runs over all N

elements.

The second way is “matrix-bordering” (Faddeev & Faddeeva, 1963, Strang & Fix, 1973

). To solve

(B.28)

A x = f,

write A in the block form

A11 A12

(B.29)

A=

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)

(B.30)

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

(B.31)

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

(B.32)

=

L21 L22 0 U22 A21 A22

where both L & U are triangular and

(B.33)

L11 U11 = A11

(B.34)

L21 U11 = A21

(B.35)

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

T T T

(B.37)

U11 L21 = A21

APPENDIX B. DIRECT MATRIX-SOLVERS

524

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

reality.

B.7. PARTING SHOTS 525

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

accurate.

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

(C.1)

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

526

C.1. INTRODUCTION 527

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,

then

’1

a(i+1) = a(i) ’ J (C.3)

F

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

2

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

(C.6)

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:

2

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)

(i+1)

(C.8)

ux

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)

(x)

APPENDIX C. NEWTON ITERATION

528

Then (x) is the solution of

’ {Fu } = F ’ ux

(i)

(C.10)

x

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)

’0

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.