. 109
( 136 .)


the number of parameters by two.
For some differential equations, solutions exist only for discrete values of the eigenpa-
rameter. An example is van der Pol™s equation:
utt + u = “ ut (1 ’ u2 ) [Van der Pol™s Eq.] (C.36)
If we replaced the term in the () by a constant, then (C.36) would describe a harmonic
oscillator which was damped or excited depending on the sign of the constant. Because this
“constant” in (C.36) is actually a function of u(x), we see that the solution of the nonlinear
equation is self-excited when u is small and self-damping when u is large. When (C.36) is
given an arbitary initial condition, the solution eventually settles down to a nonlinear limit
cycle in which the damping and self-excitation balance out over the course of one period of
the cycle to create a permanent, self-sustaining oscillation.
We can solve for the period and structure of the limit cycle as a nonlinear eigenvalue
problem by adopting the following device. Let ω denote the (unknown) frequency of the
limit cycle. If we make the change of variable
x ≡ ωt, (C.37)

then the limit cycle will automatically be periodic with a period of 2π and we can compute
a solution by assuming a Fourier series in x1 . The problem is then
ω 2 uxx + u = “ ω ux (1 ’ u2 ) (C.38)
In contrast to the undamped Korteweg-deVries equation, it is not possible to pick an
arbitrary value of ω and ¬nd a solution. The reason is that the tug-of-war between the
damping and the self-excitation normally has only a single solution for a given value of
the parameter “. It is not possible to have arbitrarily large amplitudes for u, for example,
because the equation would then be damped over almost the whole cycle. If we started
the calculation with a very large value of u(t = 0), we would ¬nd that u(t) would quickly
decay into a limit cycle of moderate amplitude and a unique frequency ω. Consequently,
just as in a linear eigenvalue problem, we must vary both u(x) and the eigenparameter ω
until we have obtained a unique solution.
Van der Pol™s equation happens to also have the property of translational invariance. In
this respect, it is like the Korteweg-deVries equation: its solution is not unique, but rather
is a continuous family, and the linearized differential equation of the Newton-Kantorovich
has an eigenvalue with zero eigenfunction. Thus, the matrix which is the pseudospectral
discretization of the linearized ODE will be almost singular [exactly singular in the limit
N ’ ∞] and the iteration will diverge unless we add additional constraints or conditions
to restrict the solution to a unique member of the continuous family of solutions.
Translational invariance may be destroyed by imposing a phase condition such as
[Phase Condition] (C.39)
ux (x = 0) = 0
which demands that the origin be either a maximum or a minimum of u(x). For the KdV
solutions, it is possible to impose (C.39) implicitly by restricting the basis set to functions
which are symmetric about x = 0.
However, because of the presence of both even and odd order derivatives in van der
Pol™s equation, its solution is not symmetric about any point and we cannot impose a phase
condition by restricting the basis set; we have to use both sines and cosines. To obtain a
unique solution, we must therefore impose the phase condition (C.39) explicitly.
Since a solution does not exist for arbitary ω, but only for one particular frequency, we
must retain ω as an additional unknown. Thus, if we apply the pseudospectral method
with N coef¬cients and N interpolation points, we must solve an (N + 1) — (N + 1) matrix
equation at each iteration. The (N +1)-st row of the matrix is the spectral form of (C.39);
the (N +1)-st unknown is δω, the correction to the frequency at the i-th iteration.
The Newton-Kantorovich iteration is a modest generalization of that used for earlier
examples. We make the simultaneous replacements
’’ (C.40a)
u(x) u(x) + (x)
’’ (C.40b)
ω ω + (δ ω)
in (C.38), collect powers of , and then the iteration, as for the other examples, consists of
equating the ¬rst order terms (divided by ) to the negative of the zeroth order terms. In
particular, if we write (C.38) in the symbolic form
uxx = F (x, u, ux , ω)
then the linearized differential equation is

’ Fux ’ Fu ’ Fω (δ ω) = ’ uxx ’ F
xx x

1 Note that time-periodic solutions, usually with the frequency as an eigenparameter, are an exception to the
usual rule that spectral methods are applied only to spatial coordinates.

where F and all its derivatives are evaluated with u = u(i) and ω = ω (i) . The only extra
term is the derivative with respect to ω. We expand (x) as a Fourier series with a total of
N terms, demand that (C.42) is satis¬ed at N interpolation points, add the spectral form of
the phase condition, and then solve the linear matrix equation for the column vector whose
¬rst N elements are the N spectral coef¬cients of (x) and whose (N+1)-st element is (δω).
We then add the corrections to u(i) and ω (i) and repeat the iteration until it has converged.
It should be clear, as stressed at the end of the preceding chapter, that the Newton-
Kantorovich method is always applicable and the iteration can always be derived by very
elementary methods. (“Elementary” is a euphemism for freshman calculus! Functional
analysis is not needed). However, the differences between the Korteweg-deVries and van
der Pol cases, the need to impose phase conditions in some problems but not others, the
possibility that the eigenparameter may assume a continuous range of values or just one ”
all in all, they show that the subject of nonlinear differential equations is very complicated.
It is important to thoroughly understand the particular problem you are trying to solve
because the range of possibilities ” no solution, a solution only for a discrete value of a
parameter, solution for a continuous range of values of the parameter ” is rather large.
The only sure statement is this: given a suf¬ciently good ¬rst guess, the Newton-
Kantorovich method will always compute a solution of the indicated form if it exists. If
the solution does not exist, or you have left out phase conditions so that the solution of the
numerical problem is not unique, the iteration will (probably!) diverge, or the pseudospec-
tral matrix for the iteration differential equation will be a singular matrix.

C.4 Summary
SUCCESS: Given a suf¬ciently good ¬rst guess, the Newton-Kantorovich method will con-
verge. Furthermore, the asymptotic convergence will be quadratic in the sense that the
magnitude of the correction (i+1) will be roughly the square of the magnitude of (i) for
suf¬ciently large i.

Notes: (i) The quadratic convergence does not occur until the iteration is close to the ¬nal
solution (ii) The corrections will not decrease inde¬nitely, but will “stall out” at some very
small magnitude which is determined by the roundoff of the computer ” typically O(10’12
to 10’14 ). (iii) The iteration is self-correcting, so it is possible to have small errors in solving
the linearized ODE while still obtaining the correct answer. The signature of such errors
is geometric rather than quadratic convergence. (iv) The iteration does not converge to the
exact solution to the differential equation, but only to the exact solution of the N -basis
function discretization of the differential equation. Remember to vary N and compare!

FAILURE: Newton™s method may diverge because of the following:
(i) Insuf¬ciently accurate ¬rst guess.
(ii) No solution exists for the nonlinear problem.
(iii) A family of solutions, parameterized by a continuous variable, exists, rather than a
single, unique solution.
(iv) Programming error.

Notes: (i) Close to a “limit” or “bifurcation” point as discussed in Appendix D, a very
accurate ¬rst guess is needed. (ii) Small perturbations may sometimes drastically modify
nonlinear solutions; perturbing the KdV equation (C.28) so as to make the coef¬cients vary
with x will destroy the translational invariance, and reduce an in¬nite number of solutions
to either one or none, depending upon whether one can ¬nd a soliton such that the ¬rst
C.4. SUMMARY 535

Table C.1: Nonlinear Boundary Value Problems Solved by Spectral Methods: A Select Bib-

References Comments
Boyd(1986c) FKdV bicnoidal; Fourier pseudospectral
Yang&Keller(1986) Fourier/¬nite difference algorithm for pipe ¬‚ow; good use of
pseudoarclength continuation to follow multiple branches
Van Dooren&Janssen(1989) period-doubling for forced Duf¬ng equation
Fourier Galerkin
Haupt&Boyd(1988, 1991) FKdV cnoidal waves; KdV double cnoidal wave
Boyd(1989b) review; many examples
Boyd & Haupt(1991) Fourier pseudospectral
Heinrichs(1992d) quasi-linear elliptic problems
Boyd(1995h) multiple precision; ¬nite difference preconditioned
nonlinear Richardson™s iteration; nonlocal solitons
Mamun&Tuckerman(1995) good discussion of how a time-marching code can
be easily modi¬ed to a Nonlinear Richardson™s iteration
nonlinear steady ¬‚ow between concentric spheres
Yang&Akylas(1996) weakly nonlocal capillary-gravity solitons
Boyd (1997b) Camassa-Holm “peakons”; T B pseudospectral
Boyd (1998c) monograph on solitons; 3 chapters on spectral methods
Nicholls(1998) water waves; Fourier in 2 horizontal coords.
Schultz&Vanden-Broeck & spectral boundary integral; Newton continuation;
Jiang&Perlin(1998) water waves with surface tension

order perturbation is orthogonal to the eigenfunctions. (iii) Too many solutions is as bad as
too few: one must add “side” conditions like (C.39) until a unique solution is obtained. (iv)
Because Newton™s iteration is self-correcting, a code may pass tests for some parameters ”
and then fail in a different parameter range. Geometric convergence, that is, /
∼ ρ for some constant ρ < 1, indicates that the Jacobian matrix is inconsistent with

the N equations which are being solved. (This is true if the iteration is a “honest” Newton™s
method with recomputation of the linearized differential equation at each step; if the same
linearized ODE is used at every iteration to save money, then geometric convergence is the
expected price for this modi¬cation of Newton™s method.)
Appendix D

The Continuation Method

“Creeps in this petty pace from day to day”
” W. Shakespeare, Macbeth

D.1 Introduction
Iterative methods have the vice that they must be initialized with an approximate solution.
The “continuation” method is a strategy for systematically constructing a ¬rst guess and
tracing a branch of solutions.
The basic idea is very simple. Suppose that a problem depends upon a parameter ±.
Suppose further that we know the solution u0 for ± = 0. We can then compute u(± = δ) by
applying the Newton or Newton-Kantorovich iteration with u(0) (δ) = u0 , provided that δ
is suf¬ciently small. Similarly, we can set u(0) (2 δ) = u(δ), and so on: at each step, we use
the solution for u(n δ) as the ¬rst guess for calculating u([n + 1] δ).
The better our guess, the faster the iteration will converge. The basic scheme can be
re¬ned by using the linear extrapolation:

u(0) ([n + 1] δ) = u(n δ) + δ {u(n δ) ’ u([n ’ 1] δ)} (D.1)

In summary,

’ ’ ’’
[Analytical solution] (D.2)
u(± = 0)

Newton iteration
’ ’ ’’
’’ ’ ’ ’ ’ ’ ’’
’ ’ ’ ’’
u(0) (δ) (D.3)
u(0) u(δ)

u([n ’ 1] δ) Linear Extrapolation Newton™s iter.
’ ’ ’ ’ ’ ’’ u(0) ([n + 1] δ) ’ ’ ’ ’ ’’ u([n + 1] δ)
’ ’ ’’ (D.4)
u(n δ)

In words, we have the following.


De¬nition 50 (CONTINUATION METHOD) This procedure solves the equation

N (u; ±) = 0 (D.5)

for u by marching in small steps in ±. N is the nonlinear operator and ± is the parameter. The
march begins with a known solution ” typically, a linear solution which can be found analytically
or numerically ” and applies an iteration to compute u at each value of ±, USING PREVIOUSLY
N may include derivatives or integrations, and may represent a system of nonlinear equations.
The unknown u may be the solution to either an algebraic equation or system of equations, or the
solution to a differential equation or system of such equations.

Strictly speaking, “continuation method” should be “continuation methods” in De¬-
nition 50 because the basic idea of parameter-marching to obtain the ¬rst guess can be
applied with a wide variety of iteration and extrapolation schemes. We can always replace
Newton™s method by the secant method, for example; since the latter has a smaller radius
of convergence than Newton™s, the price we must pay is a smaller maximum step in ±.
Similarly, to de¬ne better extrapolation formulas, we can use either ODE integration
schemes or higher-order Lagrangian interpolation in ± ” very carefully since interpolation
may diverge as the degree of the interpolating polynomial becomes large, especially when
The ODE methods are based on differentiating N (u[±], ±) with respect to ± to obtain

Nu (u, ±) u± + N± = 0 (D.6)

For the special case that N (u, ±) is a single algebraic equation, (D.6) is the ¬rst order ordi-
nary differential equation

du ‚N (u, ±)/‚±
=’ [“Davidenko equation”] (D.7)
d± ‚N (u, ±)/‚u

Given a single initial condition such as u(± = 0) = u0 , we can apply our favorite Runge-
Kutta or predictor-corrector method in ± to trace the solution branch.
This differential-equation-in-the-parameter is named after Davidenko who ¬rst observed
that nonlinear algebraic equations could be solved by integrating an ODE or system of
ODE™s in the parameter. When extended to systems of M nonlinear equations in M un-
knowns, u and N become M -dimensional vectors and ‚N /‚u becomes the Jacobian ma-
trix whose elements are given by ‚Ni /‚uj . Integrating (D.7) is then relatively expensive
because one must compute du/d± by Gaussian elimination.
In principle, Newton™s method can be completely replaced by the integration of the
Davidenko equation. However, the errors in a Runge-Kutta scheme will gradually accu-
mulate as we march, so it is a good idea to periodically apply Newton™s method so as to
reinitialize the march with a u(±) that is accurate to full machine precision.

D.2 Examples
Let the goal be to solve the polynomial equation


. 109
( 136 .)