544

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

whenever

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.

D.5. BIFURCATION POINTS 545

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

±’1

du

(D.32)

=

u’1

d±

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

APPENDIX D. THE CONTINUATION METHOD

546

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

in¬nities.

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.

(D.33)

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±

(D.35)

= +

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

(D.37a)

S1 (u, ±; s) = 0

du(s0 ) d±(s0 )

S1 ≡ [u(s0 + δs) ’ u(s0 )] + [±(s0 + δs) ’ ±(s0 )] ’ δs (D.37b)

ds ds

D.6. PSEUDOARCLENGTH CONTINUATION 547

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

2

d±

(D.39)

=1 [θ = 0],

ds

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; ±)

u

(D.40)

=

’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

scalar.

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

(D.41)

=

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

via

T

(du/dt, d±/dt)T

du d±

(D.42)

, =

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:

APPENDIX D. THE CONTINUATION METHOD

548

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

curve.

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

continuation.

D.6. PSEUDOARCLENGTH CONTINUATION 549

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

Transformations

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

(E.1)

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