. 67
( 136 .)


(MRR) variant. The imaginary parts of the »j simply mean that the decay will be oscillatory
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
r (u) = 0
J(un+1 ’ un ) = r(un ) (15.107)

where J is the usual Jacobian matrix whose elements are
Jij =
‚uj u=un
This is identical with the result of integrating the system of differential equations

du 1
=’ [“Newton ¬‚ow” equation] (15.109)
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)
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.

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
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
= ’H J(u ’ u0 ) + O (u ’ u0 )2 (15.114)

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

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

“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
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


(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

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)
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
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:

(i) (16.6)
Tn (x) = n Un’1 (x)
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


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

Table 16.1: A sample subroutine to compute the Chebyshev polynomials and their deriva-
tives via the derivatives of the cosine functions.

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.


. 67
( 136 .)