. 111
( 136 .)



Keller and Reinhart, 1985, and Keller, 1992) have developed a “pseudo-arclength continu-
ation method” which can follow a solution smoothly around a limit point. The “globally
convergent homotopy” methods reviewed by Li (1987) are very similar. We postpone a
discussion until Sec. 6.

D.5 Bifurcation points
As in the previous section, the de¬nition below and the examples are restricted to a single
equation in a single unknown.

De¬nition 52 (BIFURCATION POINT) This is a point in parameter space where two branches
cross. The formal de¬nition for a single equation N (u; ±) = 0 is that a bifurcation point occurs
Nu (u[±]; ±) = 0 N± (u[±]; ±) = 0, (D.29)
which differs from the de¬nition of a limit point in that both derivatives are 0 instead of just one.
Newton™s method fails in the vicinity of a bifurcation point; the iteration literally cannot decide
which branch of the two branches to converge to.
For a system of equations N (u) = 0 where u is a vector of unknowns, (D.29) become the
conditions that (i) the Jacobian matrix Nu has a vanishing determinant and (ii) N± is within the
range of Nu so that the matrix equation Nu du/d± = ’N± has bounded solution even though the
Jacobian matrix is singular.
An example is the quadratic equation
u2 ’ 2 u ’ ±(± ’ 2) = 0 (D.30)

Figure D.4: Same as Fig. D.3 except that the quadratic has been altered so that its solution
branches have a bifurcation point instead of a pair of limit points. At a bifurcation point,
two branches cross. As at a limit point, Newton™s iteration cannot converge to either one
but instead diverges. Unlike a limit point, however, all branches have a ¬nite slope du/d±
at a bifurcation point. A limit point is one where branches merge; a bifurcation point is one
where branches cross.

Figure D.5: Schematic of “shooting a bifurcation point”. Newton™s iteration converges
when the ¬rst guess is within the cross-hatched region. Unfortunately, this convergence
region shrinks to zero width at the bifurcation point. However, linear extrapolation from
the left “x” will give a ¬rst guess as indicated by the right “x”. This lies within the conver-
gence zone on the far side of the bifurcation point. Newton™s iteration will converge from
a guess of the point marked by the right “x”, and one may then apply continuation to trace
out the whole solution branch to the right of the bifurcation point. The only complication
is the issue of how far to proceed along the straight line when extrapolating beyond the
bifurcation point; some trial-and-error may be needed.

whose two roots are
u2 (±) = 2 ’ ± (D.31)
u1 (±) = ± ;
The bifurcation point is located at ± = 1 where u1 = u2 and the two branches cross. The
ODE in ± for (D.30) is

At the limit point, the denominator of (D.32) tends to 0 ” but so does the numerator.
The slope du/d± is always ¬nite. Consequently, although we cannot compute u(±) at the
bifurcation point, we can calculate either branch without using a small step-size even at
points very near ± = ±bifurcation .
Fig. D.5 shows how it is possible to, in Keller™s words, “shoot the bifurcation point”.
The radius of convergence of Newton™s method (in u) shrinks to 0 at the bifurcation point,
but expands on either side. Consequently, if we use a rather large step in ±, we can jump
from the region of convergence on one side of the bifurcation point to the corresponding
region on the other side. We do not need tricks like “pseudoarclength continuation” near
a bifurcation point.
A harder task is to shift to the other branch, that is, to jump from u1 (±) to u2 (±). By
using a local Taylor approximation near the bifurcation point (all the needed terms can

be computed numerically), we can obtain a ¬rst guess for points on the other branch very
close to the bifurcation. All one needs is a single point on u2 (±); one may then trace the rest
of the branch by ordinary continuation. Seydel (1988) and Keller (1992) are good references.

D.6 Pseudoarclength Continuation
As noted in Sec. 4, a major problem for continuation in a physical parameter ± is the ex-
istence of limit points where du/d± ’ ∞. In pseudoarclength continuation, both u and ±
are taken as functions of a new parameter s which approximates arclength along the so-
lution curve in the u-± plane. The gradient (du/ds, d±/ds) is always ¬nite so there are no
The arclength is de¬ned by approximating the solution curve u(±) by many short line
segments. In the limit that the segments become in¬nitely short, the sum of the length of
the segments between two points on the curve is the arclength between those points.
Mathematically, suppose that the solution curve is speci¬ed parametrically, i. e.
u = u(t) & ± = ±(t)
Let t denote the change in t between the endpoints of one of the in¬nitesimal line seg-
ments that we use to measure arclength. Let u and ± denote changes in u and ± be-
tween the ends of the segment. By the Pythagorean theorem,
( s)2 = ( u)2 + ( ±)2 (D.34)
which, dividing by t and then taking the limit, is
2 2 2
ds du d±
= +
dt dt dt
Our goal is to parameterize the solution curve in terms of the arclength so that t ≡ s.
This implies that the L. H. S. of (D.35) is one, and this in turn gives a constraint which can
be applied to “in¬‚ate” our system of M unknowns, u(±), into a system of M + 1 unknowns
(u(s), ±(s)) by giving us the (M +1)-st nonlinear equation:
2 2
du d±
[“Arclength Constraint”] (D.36)
+ =1
ds ds
For simplicity, we have written u as a scalar, but the only change if u is a vector is that
(du/ds)2 = (du1 /ds)2 + (du2 /ds)2 + . . . , i. e. is the usual scalar product of the vector du/ds
with itself.
There is one modest complication: (D.36) involves the s-derivatives of u and ± rather
than these quantitities themselves. Strictly speaking, when u is a vector of dimension M ,
the system N (u; ±) = 0 plus (D.36) is a system of (M +1) equations in (2M +2) unknowns!
However, H. Keller, who invented this technique, pointed out that we are not interested
in the arclength for its own sake, but rather as a computational device for evading limit
points. It follows that it is legitimate to approximate the s-derivatives in (D.36) in any way
such that u(s) & ±(s) become the unknowns. When such approximations are made, s is
no longer exactly equal to the arclength, but is merely an approximation to it. It is for this
reason that this family of algorithms is called “pseudoarclength” continuation.
Many approximations are possible. The simplest is a forward ¬nite difference approxi-
mation to (some of) the derivatives in (D.36), which gives
S1 (u, ±; s) = 0
du(s0 ) d±(s0 )
S1 ≡ [u(s0 + δs) ’ u(s0 )] + [±(s0 + δs) ’ ±(s0 )] ’ δs (D.37b)
ds ds

where s0 is a constant which is normally the value of the arclength at the previous step of
the continuation, and δs ≡ s ’ s0 .
An alternative is useful in overcoming a drawback of standard pseudoarclength con-
tinuation: Because s is the marching parameter and ± is an unknown, there is no simple
way to compute u(±) at a given “target” value of ±. The problem is that we do not know
which value of s corresponds to the target ±. However, if we generalize (D.37) to
du d±
Sθ ≡ θ [u(s0 + δs) ’ u(s0 )] + (2 ’ θ) [±(s0 + δs) ’ ±(s0 )] ’ δs (D.38)
ds ds
then Sθ (θ = 1) = S1 . However, when θ = 0, (D.38) becomes an approximation to

=1 [θ = 0],
that is, s becomes identical with ±. We have switched back to ± as the marching parameter,
and it is trivial to march to a desired “target” value of ±.
Thus, it is possible to switch from ordinary continuation (in ±) to pseudosarclength
continuation in s and back again merely by changing θ. In the vicinity of a limit point, of
course, taking θ = 0 will lead to disaster, but very close to a limit point, we have no hope
of marching to a target ± anyway. There may not even be a solution if ±target > ±limit . Away
from a limit point, however, we can apply pseudosarclength continuation to march close
to a target value of ± and then switch to θ = 0 to march directly to ±target on the ¬nal step.
In any event, Newton™s iteration for the arclength continuation requires solving the
(M + 1) — (M + 1) matrix equation which in block form is
Nu N± ’N (u; ±)
’S(u; ±)
Su S± ±
When u is an M -dimensional vector, Nu is the M — M Jacobian matrix whose elements are
(Nu )ij ≡ ‚Ni /‚uj . Similarly, N± is the M -dimensional column vector whose elements are
(N± )i ≡ ‚Ni /‚±, Su is the M -dimensional row vector with (Su )j ≡ ‚S/‚uj , and S± is a
The one minor technical complication is that elements of S require the s-derivatives of
u and ±. Fortunately, we can obtain them in two steps. First, solve the system
Nu N± du/dt 0
0 1 d±/dt 1
which modi¬es only the R. H. S. and the last row of (D.40). Second, rescale du/dt and d±/dt
(du/dt, d±/dt)T
du d±
, =
ds ds [du/dt]2 + 1
The justi¬cation for (D.41)“(D.42) is that along the solution curves, N (u[t], ±[t]) = 0
where t is any parameter that parameterizes the solution curve. It follows that dN /dt ≡ 0,
but this is merely the top M rows of (D.41). This condition, dN /dt = 0, gives only M con-
straints on (M +1) unknowns, so to obtain a system with a unique solution, we arbitrarily
demand that d±/dt = 1, which is expressed by the bottom row of (D.41). For small dis-
placements along the solution curve, the parameter t must be proportional to s so that the
t-derivatives of u and ± are proportional to the desired s-derivatives. We can determine
the proportionality constant by recalling that the length of the vector (du/ds, d±/ds)T is
unity, which implies (D.42).
Thus, the overall algorithm may be summarized as:

1. Generate a ¬rst guess for u & ±. [Predictor step]
2. Apply Newton™s method to compute u & ± by repeatedly solving the (M +1)—(M +1)
matrix problem (D.40). [Corrector step]
3. Modify the last row of the square matrix and the R. H. S. of (D.40) to give (D.41).
Solve (D.41) and apply the rescaling (D.42) to compute du/ds and d±/ds.
4. Increase s ’ s + δs and return to 1.

When s = 0, the ¬rst guess is usually ± = ±0 , u = u0 (±0 ) where u0 is a known analytical
solution (often a linear solution). For larger s, the simplest predictor is
du(s0 ) d±(s0 )
u(s0 + δs) ≈ u(s0 ) + & ±(s0 + δs) ≈ ±(s0 ) + (D.43)
δs δs
ds ds
but many alternatives are possible. Bank & Chan (1986) point out that a good predictor
can greatly reduce the computation time by reducing the number of Newton iterations and
apply a rather ad hoc but effective procedure that requires solving two simultaneous equa-
tions in two unknowns. Generating a good ¬rst guess is still a subject of active research,
but (D.43) is safe and reliable.
If the ¬rst guess for s = 0 is exact, one may skip the Newton™s iteration for this one
parameter value. If, however, the march is initialized with an approximate solution, then
(D.40) cannot be applied because it requires s-derivatives which we do not yet know. The
remedy is to modify (D.40) by replacing the last row by of the square matrix by (0 1)T
and setting the last element of the R. H. S. equal to 0. This modi¬ed matrix problem is
then the classical Newton™s iteration with u as the unknown and ± ¬xed; once u(± = ±0 )
has been re¬ned to the desired degree of accuracy, one may compute the corresponding s-
derivatives in step 3 and solve (D.40) without modi¬cation at all later points on the solution
The extra computation required by pseudoarclength continuation versus the classic
Newton™s method is very minor, especially when M 1. By solving (D.40) and (D.41)
via the 2 — 2 block decomposition which is described in Appendix B, one may solve (D.41)
for the s-derivatives without having to recompute the full LU decomposition of a large
matrix; the factorization of the Jacobian matrix Nu , which is the costly step, can be applied
to both (D.40) and (D.41).
The block-factorization fails near limit points because the Jacobian matrix Nu is singular
at a limit point. This is not a serious dif¬culty in practice, however. Near the limit point,
the solution varies with s proportional to the singular (zero eigenvalue) eigenfunction of
the Jacobian anyway, so the 2 — 2 block solution of (D.41) will give the correct s-derivatives
even very near the limit point.
Unfortunately the accuracy of the Newton™s step (D.40) is degraded if the block fac-
torization is used too near a limit point. The remedy is to switch to standard Gaussian
elimination for the entire in¬‚ated system. The beauty of the pseudoarclength continuation
method is precisely that the in¬‚ated, (M +1)-dimensional matrix system is not singular even
when the M -dimensional Jacobian matrix is. Bank & Chan (1986) describe a method for
obtaining accurate solutions via block-factorization for (D.40), but it would take us too far
a¬eld to describe their technique in detail.
It is important to note that pseudoarclength continuation eliminates only limit points.
At bifurcation points, continuation still has dif¬culties. However, we can “shoot” the bifur-
cation points as described in the previous section. It is comforting, however, that Keller has
found that “shooting” tends to be easier with arclength continuation than with standard

Seydel (1988) is a good treatment of continuation, limit points, and bifurcation points
which assumes only an undergraduate mathematical background: ordinary differential
equations, partial derivatives, and matrix algebra. Keller (1992) is an outstanding mix of
review and research papers at a higher mathematical level.
Appendix E

Change-of-Coordinate Derivative

“A Chebyshev polynomial is a Fourier cosine in disguise”
” J. Boyd

In many problems, it is useful to replace the original coordinate y by a new variable x:

y = f (x)

One common application is to convert a Chebyshev series in y into a Fourier cosine series
in x, in which case f (x) ≡ cos(x). In this appendix, we collect (mostly unpublished) tables


. 111
( 136 .)