instead of monotonic for some or all of the eigenmodes.

Alas, even squaring the matrix won™t help when the dif¬culty, as for the Helmholtz

equation (15.99), is an eigenvalue which is approximately zero. (Squaring the matrix makes

the corresponding eigenvalue even smaller!) Still, the conclusion is that iterative methods

can be applied ” with a little care ” to the vast majority of real-life problems including

those with eigenvalues of both signs.

15.14 Nonlinear Iterations & the Preconditioned Newton Flow

The usual strategy to solve a nonlinear boundary value or eigenvalue problem is, after

discretization, to apply Newton™s iteration. When the linear matrix problems of each New-

ton iteration are themselves solved by iteration, it is often desirable to avoid an iteration-

within-an-iteration by combining Newton™s iteration and the preconditioned Richardson™s

iteration into a single step. To do so, it is useful to employ a conceptual device that we used

to motivate Richardson™s method: interpreting an iteration as the temporal discretization

of a time-dependent problem8 .

Newton™s iteration for the system of N equations

(15.106)

r (u) = 0

is

J(un+1 ’ un ) = r(un ) (15.107)

where J is the usual Jacobian matrix whose elements are

‚ri

(15.108)

Jij =

‚uj u=un

This is identical with the result of integrating the system of differential equations

du 1

=’ [“Newton ¬‚ow” equation] (15.109)

r

dT J

by the Euler forward method with unit step in the pseudotime, T .

If we expand r(u) in an N-dimensional Taylor series centered on a root u0 , then (15.109)

becomes locally

du 1

=’ r(u0 ) + J(u ’ u0 ) + O (u ’ u0 )2 (15.110)

dT J

= ’(u ’ u0 ) + O (u ’ u0 )2 (15.111)

T

7 Tominimize the condition number, one should solve (15.105) as the system Λy = f, Λ x = y.

8 Chu(1988) is a good review of other examples of “Continuous realizations of Iterative Processes”, to quote

his title.

15.14. PRECONDITIONED NEWTON FLOW 321

In the language of the theory of differential equations, (15.111) is a linearization of the

nonlinear system about a critical point. The linearization shows that every root of r(u) = 0

is a local attractor. The exact “Newton ¬‚ow” converges geometrically on the root since one

can prove without approximation that

r(u[T ]) = e’T r(u[0]) (15.112)

Newton™s iteration actually converges much faster than (15.112); the number of correct

digits roughly doubles with each iteration. The reason for the faster convergence is that

if we ignore the quadratic terms in (15.111), Euler™s method with unit time step reduces

(u ’ u0 ) to 0 instead of to 36%(= exp(’1)) of its initial value. Thus, poor accuracy in

time integration brings us to the root faster than solving the Newton ¬‚ow accurately in

pseudotime.

Unfortunately, Newton™s method requires inverting the Jacobian matrix at each time

step. When r(u) = 0 is the pseudospectral discretization of a differential equation, the

Jacobian matrix is dense and costly to invert.

Fortunately, if the Newton ¬‚ow is modi¬ed by replacing the inverse-Jacobian matrix by

a matrix which is nonsingular at the roots of r(u), but is otherwise arbitrary, the roots of

r(u) are still critical points of the differential equation.

Not all approximations of the Jacobian matrix will succeed because the modi¬cations

may turn some roots from attractors into repellors. Suppose, however, the inverse-Jacobian

is replaced by the inverse of a preconditioning matrix:

du 1

= ’ r(u) [“Pre-conditioned Newton ¬‚ow”] (15.113)

dT H

Linearization near a root gives

’1

du

= ’H J(u ’ u0 ) + O (u ’ u0 )2 (15.114)

dT

Orszag (1980) showed that if H is the ¬nite difference approximation to J, then all the

eigenvalues of (1/H)J will lie on the interval ± ∈ [1, 2.47] (as already explained in Sec. 4). If

the error (u ’ u0 ) is decomposed into the eigenvectors of (1/H)J, then the error in different

components must decay at rates ranging between exp(’T ) and exp(’2.47 T ).

Because of these varying rates of decay, it is not possible to ¬nd a magic time step that

will keep the quadratic, digit-doubling convergence of Newton™s iteration. However, the

Euler forward scheme does give very rapid convergence: with a time step of 4/7, one can

guarantee that each component of the error will decrease by at least a factor of 2.33 at each

iteration ” just as fast as for a linear problem. Since the ¬nite difference matrix can be

inverted in O(1/N 2 ) of the cost for inverting the pseudospectral Jacobian where N is the

number of basis functions in each dimension, the preconditioned Newton ¬‚ow is much

cheaper than the classical Newton™s iteration.

The iterative form of (15.114) is

H(un+1 ’ un ) = ’„ r(un ) [Nonlinear Richardson Iteration] (15.115)

where H is any acceptable preconditioning matrix for the Jacobian matrix and „ is the pseu-

dotime step. It is appropriate to dub (15.115) the “Nonlinear Richardson Iteration” because

it reduces to the standard Richardson iteration whenever r is a linear function of u.

This “preconditioned Newton ¬‚ow” derivation of (15.115) is taken from Boyd (1989a),

but the same algorithm was independently invented through different reasoning in the

CHAPTER 15. MATRIX-SOLVING METHODS

322

technical report by Streett and Zang (1984)9 . Eq. (15.115) has probably been re-invented

half a dozen times.

Simple though it is, the Nonlinear Richarson™s Iteration (NRI) is very useful. It is silly

and expensive to apply Richardson™s iteration to solve the linearized differential equations

of a strict Newton method; one is then wasting many iterations to converge on an interme-

diate Newton iterate which is of no physical interest. With (15.115), every iteration reduces

the nonlinear residual and brings us closer to the root.

Numerical examples are given in Boyd (1989a).

15.15 Summary & Proverbs

MORAL PRINCIPLE #1: Never solve a matrix using anything but Gaussian elimination

unless you really have to.

This may seem strange advice to conclude a chapter mostly on iterative methods. Nev-

ertheless, it is more than a little important to remember that there are no style points in

engineering. Iterative methods, like mules, are useful, reliable, balky, and often infuriat-

ing. When N is small enough so that Gaussian elimination is practical, don™t fool around.

Be direct!

Note further that for constant coef¬cient ordinary differential equations and for sepa-

rable partial differential equations, Galerkin/recurrence relation methods yield sparse ma-

trices, which are cheap to solve using Guassian elimination.

MORAL PRINCIPLE #2: Never apply an iteration without ¬rst preconditioning it.

Well, almost never. However, un-preconditioned iterations converge very slowly and

are much less reliable because a single eigenvalue of the wrong sign will destroy them.

MORAL PRINCIPLE #3: Never use multigrid as a purely iterative method; apply Gaussian

elimination on the coarsest grid.

Besides being the difference between convergence and divergence if the matrix is indef-

inite, elimination on the coarsest grid seems to improve the rate of convergence even for

positive de¬nite matrices.

Good, reliable iterative methods for pseudospectral matrix problems are now available,

and they have enormously extended the range of spectral algorithms. Iterations make it

feasible to apply Chebyshev methods to dif¬cult multi-dimensional boundary value prob-

lems; iterations make it possible to use semi-implicit time-stepping methods to computed

wall-bounded ¬‚ows. Nevertheless, this is still a frontier of intense research ” in part be-

cause iterative methods are so very important to spectral algorithms ” and the jury is still

out on the relative merits of many competing strategies. This chapter is not the last word

on iterations, but only a summary of the beginning.

9 See also Canuto et al. (1988, pg. 390).

Chapter 16

The Many Uses of Coordinate

Transformations

“There are nine and sixty ways of constructing tribal lays”

” R. Kipling, In the Neolithic Age

16.1 Introduction

A change of coordinates is a popular tool in elementary calculus and physics. Such trans-

formations have an equally important role in spectral methods. In the next section, the

Chebyshev polynomial-to-cosine change-of-variable greatly simpli¬es computer programs

for solving differential equations.

However, there are many other uses of mapping. In this chapter, we concentrate on

one-dimensional transformations, whose mechanics is the theme of Sec. 3. In Sec. 4, we

show how in¬nite and semi-in¬nite intervals can be mapped to [-1, 1] so that Chebyshev

polynomials can then be applied. When the ¬‚ow has regions of very rapid change ”

near-singularities, internal boundary layers, and so on ” maps that give high resolution

where the gradients are large can tremendously improve ef¬ciency. We discuss a number of

speci¬c cases in Secs. 5-7. In Sec. 8, we discuss two-dimensional mappings and singularity

subtraction for coping with “corner” singularities. Finally, in the last part of the chapter, we

give a very brief description of the new frontier of adaptive-grid pseudospectral methods.

16.2 Programming Chebyshev Methods

The Chebyshev polynomials can be computed via the 3-term recurrence relation

T0 ≡ 1 T1 ≡ x

;

(16.1)

Tn+1 = 2 x Tn ’ Tn’1 n = 1, . . . .

This is numerically stable, and the cost is O(N ) to evaluate all the polynomials through

TN at a given x. Since the derivatives of the Chebyshev polynomials are Gegenbauer poly-

nomials, it follows that we can evaluate these derivatives via three-term recurrences, too

323

CHAPTER 16. COORDINATE TRANSFORMATIONS

324

(Appendix A). Indeed, through a similar three-term recurrence derived from (16.1), (Fox &

Parker, 1968 and Boyd, 1978a), we can sum a Chebyshev series in O(N ) operations without

evaluating the individual Tn .

As simple and effective as these recurrences are, however, it is often easier to exploit

the transformation

(16.2)

x = cos(t)

which converts the Chebyshev series into a Fourier cosine series:

Tn (x) ≡ cos(nt) (16.3)

By using the tables given in Appendix E, we can easily evaluate the derivatives of the

Chebyshev polynomials directly in terms of derivatives of the cosine and sine. For example

dTn (x) n sin(nt)

(16.4)

=

dx sin(t)

This is just the trigonometric representation of n Un’1 (x) where Un (x) is the Chebyshev

polynomial of the second kind de¬ned in Appendix A. Similarly

d2 Tn (x) 1

sin(t) ’n2 cos(nt) ’ cos(t) [’n sin(nt)] (16.5)

= 3

2

dx sin (t)

Table 16.1 is a short FORTRAN subroutine that computes the n-th Chebyshev polynomial

and its ¬rst four derivatives. Its simplicity is remarkable.

One modest defect of (16.4) and (16.5) is that the denominators are singular at the end-

points. However, this is irrelevant to the Chebyshev “roots” grid because all grid points

are on the interior of the interval. The problem also disappears for the “extrema” (Gauss-

Lobatto) grid if the boundary conditions are Dirichlet, that is, do not involve derivatives.

When boundary derivatives are needed, one must modify the subroutine in Table 16.1

by adding an IF statement to switch to an analytical formula for the endpoint values of the

Chebyshev derivatives of all orders:

d

(i) (16.6)

Tn (x) = n Un’1 (x)

dx

p’1

n2 ’ k 2

d p Tn n+p

(ii) (16.7)

= (±1)

dxp 2k + 1

x=±1 k=0

Instead of using formulas like (16.3) to (16.5) to compute in the original coordinate x,

one may alternatively convert the original problem into an equivalent differential equa-

tions on t ∈ [0, π] and then solve the result using a Fourier cosine series. For example,

uxx ’ q u = f (x) x ∈ [’1, 1] (16.8)

u(’1) = u(1) = 0

becomes

sin(t) utt ’ cos(t) ut ’ q sin3 (t) u = sin3 (t) f (cos[t]) (16.9)

t ∈ [0, π]

u(0) = u(π) = 0

Which of these two forms is simpler is a matter of taste. My personal preference is to

solve the problem in x, burying the trigonometric formulas in the subroutines that evaluate

derivatives.

16.3. THEORY OF 1-D TRANSFORMATIONS 325

Table 16.1: A sample subroutine to compute the Chebyshev polynomials and their deriva-

tives via the derivatives of the cosine functions.

SUBROUTINE BASIS(N,X,TN,TNX,TNXX,TNXXX,TNXXXX)

C INPUT: X (on interval [-1, 1] )

C OUTPUT: TN, TNX, TNXX, TNXXX, TNXXXX are the N-th Chebyshev

C polynomials and their first four derivatives.