Computers them. Engng Vol. 20. No. 617. pp. Ml-655.

1996
Copyright 0 1996 Elscvicr Scicncc Ltd
Pergamon Printed in Great Britain. All rights rcscrvcd
,a /: 0098-1354(95)00199-9 009%1354/96
$15.00+0.00

A ROBUST PATH TRACKING ALGORITHM FOR

.i. - HOMOTOPY CONTINUATION ,.
S. H. CHOI,~ D. A. HARNEY and N. L. BOOKS
Department of Chemical Engineering, University of Missouri-Rolla, Rolla. MO 65401, U.S.A.

(Received 7 April 1995; received for publication II May 1995)

Abstract-Conventional path tracking algorithms used in homotopy continuation systems sometimes

miss roots owing to jumping from one segment of the homotopy path to another even if there exists
homotopy paths to the roots. A robust path tracking algorithm is proposed which loses significant
efficiency only on those portions of the path where segment jumping is likely to occur.
The method presented here basically performs the predictor-corrector procedure using the Euler
predictor and the Newton corrector. Any available algorithm can be used to control the step size.
Robustness is achieved by adding the following rule: control the step size so that each continuation step
causes a reasonably small change in the determinant of the augmented Jacobian. Case studies have
shown that allowing -50 to + 100% change virtually eliminates segment jumping in all the path-
tracking algorithms tested.
The determinant monitoring step size control algorithm can be applied to most of the currently
available path tracking algorithms so that extremely tangled homotopy paths can be traced, finding all
roots on them.

1. INTRODUCTION
._.

Homotopy continuation methods find roots of

systems of nonlinear equations by tracking the
homotopy path from one or multiple starting points.
There are now software packages such as CONSOL
(Morgan, 1987) or HOMPACK (Watson et al.,
1987) that can determine all real and complex roots
of systems of polynomial equations without requir-
ing any user-supplied initial guesses. These systems
use the multiple paths algorithm for polynomial
systems for which the existence of homotopy paths t=o t=1

that lead to all roots is guaranteed (Morgan, 1987). Fig. 1. Segment jumping in the multiple paths algorithm.
Given a known point on or near a homotopy path,
a path tracking algorithm follows the path by deter- than the distance between the two points. Figure 1
mining new points near the path. Most of the cur- schematically illustrates segment jumping in the
rently used path tracking algorithms are based on multiple paths algorithm (Morgan, 1987).
the predictor-corrector procedure, and use step size Generally, if segment jumping occurs, the complete
control algorithms based on the convergence pattern path may not be traced, and roots or other points of
of the corrector (Allgower, 1981) or the curvature of interest on the untraced segments may be missed.
the path currently being traced (den Heijer and Allgower and Georg (1980) indicate that the sign
Rheinboldt, 1981). As noticed by Allgower and of the determinant of the augmented Jacobian
Georg (1980) and by Morgan (1987), the path track- should be constant when a path without bifurcation
ing algorithms used in homotopy continuation is correctly traced. Allgower (1981) presented a path
systems suffer in practice from the problem of seg- tracking algorithm which monitors the sign of the
ment jumping, which can be roughly defined by the determinant of the augmented Jacobian in order to
arc length between the current point and the next locate a bifurcation point where the determinant is
point on the homotopy path being much greater zero. These properties of the determinant of the
augmented Jacobian are generalized in this work.
As a result, the segment jumping problem can be
t Currently at Department of Chemical Engineering,
Chonbuk National University, Korea. avoided by the determinant monitoring algorithm
$ To whom all correspondence should be addressed. presented here.
641
648 g. H. CHOI et ai.

2. HOMOTOF’~‘~~6NTINUATION METHODS 3. PREDICTOR-CORRECTOR PROCEDURE

Consider a system of equations

3. I. Euler predictor _ _ -~~
f(x) = 0 (1)
Hereafer, the following vector represents a point
where .f:..~D- c R”-+ R” is a continuous and smooth ,.
in the continuation space.
function. A homotopy is defined by a continuous

[I
function h(x, t) such that h(x, 0) = g(x) and h(x, 1) = X
f(x), where teR is the homotopy parameter, and WE
E R”f’.
t
g: E c R” -+ R” is a continuous and smooth function
such that a solution to g(x) = 0, x0 E D c E is known. A point near the homotopy path after k continua-
The linear homotopy function is defined by tion steps is represented by

h(x, t) = tf(x) + (1 - t)g(x) = 0 (2) Wk= w(sk) (4)

This function is usually defined for t E [0,11, and where sk(kaO) is the arc length traced from the
then it is called a convex homotopy, which means starting point to the current point, and for k> 1.
that h is a convex combination off and g. For convex ASk-,=Sk-Sk-,=~IWk-Wk-‘~~~. (5)
homotopies, each curve satisfying equation (2) and
connecting t = 0(x = x0) and t = 1(x = x*) is called a The Euler predictor is defined by the unit tangent
homotopy path. If we deal with the entire set of vector times the step size, i.e., the first order predic- --
curves of a homotopy system h(x, t) = 0 regardless of tor. The unit tangent vectors dw/d.r at s = Sk satisfy
the value of t, h is called a global homotopy and ‘< .-:
h-‘(O) is called global homotopy paths. Homotopy
continuation methods find a solution to f(x)=0 by
(6)
tracking the path from t = 0 [g(x) = 0; x=x”] to t = 1
[f(x) = 0; x=x*], if such a path exists.

2.5

2.0

1.5

w2 1.0

0.5

0.0

-0.5
-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0

Wl

Fig. 2. Results of curve tracking for Example 1.

 Path tracking algorithm 649
where en+’is a II + 1 dimensjopal vector of which the locally convergent. If the predicted point is suffi-
n + 1st component-equals one and the other compo- ciently close to the curve, the corrector iteration is
nents are zero. Furthermore, to keep the path track- guaranteed to converge to a point on the exact curve
ing direction consistent, the unit tangent vectors for (den Heijer and Rheinboldt, 1981). Note that if H’L
k > 1 must satisfy the following condition. is sufficiently close to the homotopy path, equation
(11) with a sufficiently small 6sk gives w[+’ which is .’
(7) sufficiently close to the homotopy path.
A point after i iterations of the corrector process
If the norm of the tangent vector is left to later is represented by w:+’ (ia0). The Newton correc-
adjustment, the bottom row of the matrix in equa- tor, Awf+‘, is obtained by solving the following
tion (6) can be replaced by a known row vector to system of equations.
make the system linear. Two kinds of replacements

.(12)
have been adopted in this work.
3.1 .I. Standard Base Vector (den Heijer and
Rheinboldt, 1981). The most popular method to
calculate a tangent is to solve
A corrected point is obtained by applying the full
step as follows.
vk,entl
$:I’ = w;+’ + Aw;+‘. (13)
This corrector process is repeated until the following
where ej(l c j< II + 1) is the jth coordinate standard
condition is satisfied.
base vector, which is a n+ 1 dimensional vector of
which the jth component is one and the other IIh(wfZ:f:)lb< a,, max(cSingIdet A(w”),
components are zero.
.. If [dwld.r(q)lT (ej)#O, this ]det A(wf:;)l”“) (14)
equation gives a tangent vector vk. In our case, the
where E,,, and Eying are small positive numbers, and A
above equation is used only for k = 0.
represents the matrix in equation (6). The final
3.1.2. Previous Tangent Vector (Watson et al.,
determinant is not raised to the power of l/n if it has
1987). The formula to be used for k 2 1 is as follows.
an absolute value less than 1. If the above condition
is satisfied, the latest corrected point is accepted as a
new point near the homotopy path as follows.
ticen+‘.
(9) Wk+r= w;2,‘. (15)
The errors of new points are independent of the
continuation history. Therefore, the parameter E,,,
If [dW/dY(S(Sk_,)]TdW/d.S(Sk)T#O, this equation
does not have to be extremely small during the
gives a tangent vector vk. The unit tangent vector is
continuation process, and can be reduced to the
obtained by
desired accuracy just near the points to be located.
k
dw Experimental results suggest the value of cCur= 10m5
,(sk,=&2 (10) during the continuation process and ccur= lo-’
during the search for a target point. The parameter
where ok = + 1 or - 1 is the direction factor. For
eSingis used to numerically define a singular point,
k b 1, the tangent vector in the forward direction,
and the experimental results suggest use of the value
which satisfies inequality (7), is automatically
of E.inp= lo-‘. These values work well unless the
obtained by using equations (9) and (10) with ok =
problem is poorly scaled.
+1.
The matrix in equation (6) is referred to as the
The next point near the homotopy path (predicted
augmented Jacobian of h (Allgower, 1981) at w=
point) is
wk, and represented by A(d), i.e.
dw
w^;+‘=w”+--(Sk)& (11)
ds
where dsk > 0 is the step size. A[w(s)] = (16)

3.2. Newton corrector

A more accurate point near the homotopy path and except at the starting point, is approximated as
can be obtained by Newton’s method, which is follows.
650 !j. H. CHOI et al.

increase the step size only after five successive con-

tinuation steps. Allgower and Georg (1980) pro-
(17) p0sed.a somewhat more sophisticated rule: for con-
traction factors a, and a2 such that 0 <d, < ai Ql , -_

(1) if 11
Aw:+‘112>a211Aw~+‘llz,then dsk~bsli/2
The augmented Jacobian has the following (2) or if I~Aw~+‘~~~<~,IIAw(;+‘~~~,
then r.&+, =26sk
properties. (3) otherwise, &+, = dsli .

(1) Reversing the path tracking direction changes Corvalan and Saita (1991), using an analysis ofthe
the sign of det A(w). Newton-Kantorovich theorem, attempt to keep the
(2) Rank A(W) = rank ahlaw + 1. Therefore number of corrector steps at each continuation point
det A(w) = 0 if and only if w is a singular point, constant. They claim that in order to achieve con-
where all II X II minors of ah/% are singular. vergence after a specified number of corrector steps
(3) If det A(w’) det A(w’+‘) < 0, a bifurcation point N, the error of the predicted point must lie between
is located on the arc between & and &+I upper and lower bounds given by
(Allgower, 1981).
II,e _ w;ll = [E;:W@/$ WW]
(1%
A corrector in a different direction can be
obtained by using a different bottom row of the for some 9 such that O< 0 s 1. & is a function of the
matrix in equation (12), for example, a standard successive errors of the Newton corrector steps and
base vector as in equation (8). Generally, it is E,, is the specified convergence tolerance. They use
efficient to use the same type of matrix for correc- this expression to yield an estimate of the error of
tion as for prediction. In our case, the matrix in the next predicted point which will require_&
equation (17), if is 1, is identical to the matrix in Newton corrector steps to achieve the specified
equation (12), and can be reused in LU-decomposed convergence and then calculate the necessary step
form in the next predictor process as the matrix in size.
equation (9). The third method of den Heijer and Rheinboldt
Another numerically adequate convergence cri- (1981) uses the local curvature and the ratio of
terion for the corrector process instead of inequality corrector size to step size to estimate the safe local
(14), adopted by Allgower and Georg (1980), is convergence radius, pk+‘, while requiring that
certain parameters lie between calculated maximal
Ilh(wk+L)IIc&Ilah/aw(w~+I)II (18) and minimal values.
where E is a small positive number. However, this
4.2. Prevention of segment jumping
criterion causes insufficient correction near a
singular or nearly singular point, and thus frequently Any quantity q [w(s)] which is a continuous func-
results in inefficient path tracking for bifurcating tion of arc length and which equals zero when the
homotopies. augmented Jacobian is rank deficient can be used to
detect bifurcation points or indicate the existence of
a near bifurcation point. Both the determinant of
4. STEP SIZE CONTROL the augmented Jacobian and the reciprocal of its
condition number have these properties, while the
determinant has the additional property that it will
4.1. Conventional algorithms change sign when a simple bifurcation point is
passed, or a higher order one where an odd number
A reasonable aim in the construction of the steps of eigenvalues of A[w(s)] have changed sign.
along a curve is to ensure that the number of Segment jumping is likely to occur when two
corrector iterations remains approximately constant segments of the homotopy path are close together.
(den Heijer and Rheinboldt, 1981). Most of the step Compare the real domain homotopy paths (solid
size control algorithms seek this aim. Some path lines) in Figs 3 and 5. Two segments of the hdmo-
tracking systems, such as HOMPACK (Watson et topy path that are connected (Fig. 3, solid line) or
al., 1987), PITCON (Rheinboldt and Burkardt. that are unconnected (Fig. 5, solid line) can be close
1983) and the algorithm of Corvalan and Saita to each other. If a starting point is selected which is
(1991) estimate appropriate step sizes analytically. near the starting point that gives real bifurcation,
Others, such as CONSOL (Morgan, 1987), use then two segments of the path will be close together
simple step size control strategies. For example, and the variable q[w(s)] will be close to zero. We
CONSOL (Morgan, 1987) uses the five rule: use the terminology, a near bifurcation point, to
 Path tracking algdrithm 651

\
\ . t=o
\ . t=1
\
\ . dUds=O-
\
0 t=*m
- Im (2) = 0
.. --- Im(z)#O

I
\
-5 , , \
-4 -3 -2 -1 0 1 2 3 4 5 6 I I
---.‘-3.-2 -1 0 1 2 3 4 5 6
Re (q)
Fig. 3. Homotopy paths from the starting point (2.336, Re (2,)
- 2.5). Fig. 5. Homotopy paths from the staring point (2.337,
- 2.5).

(Zb) If q(w”) = l/~‘[A(w~)], then

refer to a point on a segment of the homotopy path
at. which 1q [w(s)] ( is minimum and near zero. Note if IsW’)l>lqW)lv
that inequality (14) gives more accurate corrected lq(wk)<lq(wk+‘)land ’
points as 1q [w(s)] ) becomes smaller. lq(~)lC~~i~lq(w”)l,
then begina
A tentative new point wk+’which was determined search to locate the local minimum and
by the corrector process is examined by the follow- check if this is a bifurcation point.
ing algorithm.
(3) If ti+’ is still tentative, accept it only when the
(1) If lq(wk) I< ~,~~~1q(w”)I, assume that wk is a following condition is satisfied.
singular point, and to escape there, accept wk+’
o<(q(wk+l)l//q(w*)l<~ (20)
as a new point without further test. All bifur-
cating branches from wk can be located by a where a is a user specified number such that
search algorthm such as the one presented by O-Ca < 1. If this condition is not satisfied, dis-
Kearfott (1983). card we+‘, halve the step size, and try continua-
(Za) If q(wk) = det A(wk), then tion again from wk.
(4) If a user specified number of step halvings fails
if q (wk)q (wk+‘) C 0, accept wk+’as a new point
to locate an acceptable new point, accept the
without further test. Begin a search to locate a
current wk+’ and indicate the possibility of
bifurcation point or correct segment jumping.
segment jumping by warning that a large
change in the monitoring variable, q[w(s)],
occurred on this continuation step. This rule
prohibits too small a step size near a singular
point. Generally, if Iq(w’)I is close to zero,
inequality (20) can be satisfied only by a small
step size. As a result, the step size automati-
cally becomes small around the near bifurca-
0 t=*ao
- Im (2) = 0 tion points. The larger a is, the more robust
--- Irn(Z)#O but less efficient the continuation process is.
Experimental results show that a = 0.5 works
well.
-0.4 -
4.3. Examples
The following two examples represent extremely
I ill and well behaved curves in the 2-dimensional
-0.9 -0.8 -0.7 -0.6 -0.5 -0.4
space. Although they are not practical homotopy
Re (z,) paths, similar curves in higher dimensions can
Fig. 4. Near bifurcation in the homotopy paths in Fig. 3. appear in the homotopy continuation problems.
652 S. H. CHOI et al.

Example 1-A:; f2-type curve. The following and 3, and maximum number of corrector
equation represents an extremely ill-behaved curve iterations = 4.
shown inFig. 2. (lb) Algorithm of Corvalan and Saita plus
determinant monitoring with ti= 0.5-.
h(w,, wr)= wr[4w:+(w7- l)‘- l] - 10-6=o. (21) (Monitoring of a continuation variable,
Let us.tr_ac_ethis curve from a starting point near the when the step size is being controlled by
curve (-2,O) in the direction (1,O). Since the this algorithm is achieved by temporarily
traced curve is almost a straight line except around updating the parameters a,, and 8 when a
(0,0),any conventional algorithm will increase the step size is halved.)
step size as much as possible at every continuation Algorithm of Corvalan and Saita (1991)
step, and segment jump unless it is extremely lucky. with the reciprocal of the condition
The segment jumping problem can be avoided by number (based on the Pi norm) monitored,
using the proposed step size control algorithm. The with a = 0.5.
empty squares in Fig. 2 are the continuation points Algorithm of Corvalan and Saita with
without the q monitoring phase, and the dots are the l/x,,,(A) monitored, the estimate made
continuation points obtained by determinant moni- using the algorithm of Cline et al. (1979).
toring with .a = 0.5. (24 Algorithm of Allgower and Georg (1980)
Example 2-A circle. Consider the curve repre- with &,=0.05, a,=0.05, ar=0.5 and
sented by the following equation. maximum number of corrector iterations =
4.
h(w,, w,)=w:+w:-l=O. (21) (2b-d) Algorithm (2a) with the same parameters
It can be shown that if this circle is traced clockwise, and the variable q monitored as in (l$g:
det A= - 2 anywhere on the curve. Therefore, the (3a) PITCON Version 6.1 (Rheinboldt and
determinant control phase in any step size control Burkardt, 1983), with all parameters set to
algorithm in this case does not cause any loss in the the default values.
curve tracking efficiency. Generally, when a curve (3b) PITCON with default parameter values
without near bifurcation is traced, there is only a and either the determinant or conditon
small change in det A for a relatively large step size, number monitored. (PITCON augments
and therefore the loss in efficiency caused by deter- the Jacobian with a standard base vector of
minant control is also small. dimension n + 1 in the bottom row. Hence,
a continuation variable is monitored only
5. PERFORMANCE EVALUATION when the previous point was calculated
with the same base vector augmenting the
Jacobian.)
5.1. Program HOMES
(4a) Mkthod III of den Heijer and Rheinboldt
The algorithm presented here has been coded into (1981) with 6s,,= 0.05, adh= 0.05, 6’= lo-’
an interactive subroutine HOMCON, which is dri- and xdh= 3.0.
ven by a conversational program, HOMES (HOMo- (4d-b) Algorithm (4a) with the same parameters
topy Equation Solver). This program starts either and the variable q montored as in (lb-d).
from the real domain or from the complex domain.
If the real domain calcuiation mode is selected, Test problems presented here are global homotopy
complex domain calculation is performed only when paths in the complex space (with the exception of
it. is necessary, for example, when it is told to PITCON, where all calculations are in the real
bifurcate to the complex domain from a turning domain), with near bifurcations obtained from the
point with respect to t. This program locates roots, fixed point homotopy, i.e.,
turning points, bifurcation points, and starting u(z, t) = tf(z) + (1 - t) (z - z”) = 0 (23)
points using the inverse Hermite cubic interpolation
(Watson et al., 1987). for the following system of equations. .

5.2. Test problems f,(z,,zr)=zf+z;-1.5?=0 (24a)

The following step size control algorithms were f~(z,,zr)=zr-z;+2=0. (24b)
tested by tracking several types of ill-behaved homo-
The starting points (2.336, - 2.5) and (2.337, - 2.5)
topy paths using program HOMES.
are to be used. In Figs 3-6, the solid lines represent
(la) Algorithm of Corvalan and Saita (1991) the curves in the real domain (zero imaginary parts),
with 6s0=0.05, ~,~=10-“, 6=0.05, N=2 and the dotted lines represent the projection onto
 Path tracking algorithm 653

the real domain of the curves’in the complex domain PITCON, which tracked the homotopy path
(nonzero imaginary parts). Each point on the pro- wholly in the real domain, located the first two roots
jection of the complex curves corresponds to a and then made a large jump onto the negative t
conjugate pair of complex vectors. portion of the homotopy path and was thus unable
to locate any subsequent roots. However, when a
5.3. Results ‘: .~_ monitoring phase was implemented into PITCON’.
The positive t portion of the real homotopy path with a = 0.5,the correct path was followed with only
in Fig. 3 was traced in the complex space. The 17 additional step halvings necessary when the
determiant of the augmented Jacobian defined by determinant was monitored and seven additional
equation (17) is - 11.22 at the starting point and step halvings ‘when the reciprocal of the condition
about +0.09 at the near bifurcation points and number was monitored.
HIKE is 4.41x 10m2at the starting point and about As can be seen from the tables, most of the step
1.4 x 10m4 at the near bifurcation points. Tables halvings take place in the Rl-Bl and B3-R4 seg-
l(a-d) show the results of path tracking with the q ments where the near bifurcation points are located.
monitoring phase. This indicates that such monitoring decreases
Algorithms (2a) and (4a) with a=0 segment efficiency significantly only in the neighborhood of
jumped twice. As a result, all roots were found near bifurcation points.
(good fortune from poor path tracking) but in a The positive t portion of the ral homotopy path in
different sequence. The correct path was traced by Fig. 5 was traced in the complex space. There are
algorithms (2a-d), for ~~20.34 when the determi- two segments unconnected in the real domain. If the
nant is monitored, for aSO. when the condition path is followed correctly, only two real roots will be
number is monitored, and for aS0.63 when an located. Algorithms (la) with N=3, (2a) and (4a),
estimate of the condition number is monitored. for a = 0, jumped from the main path to the isola,
Algorithms (4a-d) required the respective values of circled the isola and jumped back to the main path,
a to be greater than or equal to 0.49, 0.37 and 0.60 hence finding all roots. The path was followed cor-
for the correct homotopy path to be followed with rectly with determinant monitoring, for a 5 0.027
the given parameters used. using algorithm (lb), a Z=0.24 using (2b), and for
Algorithm (la), both when the desired number of a 20.42 when (4b) is used. No monitoring of the
corrector iterations, N, is two or three, segment continuation variable q [w(s)] was necessary for
jumped over the second near bifurcation point and algorithm (la) with N=2.
then returned to the starting point. The traced path PITCON, with the parameters set to their default
is a fixed point homotopy which cannot return to the values, jumped from the main path to the isola after
(non-singular) starting point if the path is correctly locating one root on the main path, found one root
followed. For N = 2, the correct path is traced when on the isola, and then jumped onto the negative t
aa0.22, 0.27, and 0.31 for determinant, condition portion of the homotopy path, from which it could
number, and estimate of condition number monitor- locate no further roots. The path was followed
ing, respectively. When N=3, the correct path is correctly using the installed determinant monitoring
traced when as0.16, 0.14 and 0.51. with a = 0.5, which caused only eight additional step
halvings.
It should be noted that in all cases, the correct
path can be followed without the monitoring phase,
by altering the relevant parameters that will ensure
step-sizes consistently small enough. These values
. dt/dsrO can be found experimentally, and for this problem
0 t=*m lead to an overall increase in the number of LU
- Im (2) = 0
--- Irn(Z)#O decompositons required for each algorithm, com-
pared to the same algorithm used with less stringent
values and a monitoring phase.
-0.4-
6.CONCLUSION

Any path tracking algorithm can be made robust

I by adding any of the three monitoring phases pre-
-0.9 -0.8 -0.7 -0.6 -0.5 -0.4
sented here. Experimental results have shown that
Rc (q) a=OS works well when the determinant and con-
Fig. 6. Near bifurcation in hte homotopy paths in Fig. 5. dition number are used to monitor the step size, and
654 S; H. CHOI et al.

Table 1, Results of pftp tracking (a) Algorithm of Allgowcr and Gcorg (1980) with ds,, = 0.05. o, = 0.C15.or = 9.5 and ,max. COT.= 4

Dcicrminant monitored .. Condition number monitored Estimate of condition number monitored

a=O.S a=0.5 a=0.65
Point iUd Cont Search Det Point LUd Cont Search Cond Point LUd Cont Search _Condest-

s 0 0 0 0 s 0 0 0 0 s 0 0 0 0
RI 112 42 10 9 Rl 58 22 7 1 RI 103 37 9 i 10
BI ‘i/289 92 I6 38 Bl 240 71 I2 32 Bl 336 96 13;. 57
R2 369 120 24 45 R2 316 102 21 34 R2 470 137 26 80
82 444 147 30 55 82 386 132 29 38 52 594 168 21 111
R3 504 169 36 63 R3 443 I58 38 40 R3 694 202 28 127
83 657 215 44 83 B3 570 199 47 48 B3 919 257 32 173
R4 840 275 53 II4 R4 732 251 55 74 R4 1363 352 40 2%

(h) Algorithm of Corvalrln and Saita (1991) with N=2. &=6.05, E,,= 1.0X 10mJ, 8=0.05 and max. car. =4

Determinant monitored Condition number monitored Estimate of condition numhcr monitored

cl= 0.5 a=0.5 a=0.65
Point LUd Cont Search Det Point LUd Cont Search Cond Point LUd Cont Search Condest

S 0 0 0 0 s 0 0 0 0 s 0 0 0 0
RI 69 30 6 0 RI 69 30 6 0 RI 69 30 6 ,O
Bl 166 64 6 8 Bl 251 92 7 23 Bl 224 86 7 17
R2 254 100 I2 13 R2 319 121 I2 26 R2 306 120 I3 23
B2 321 123 I3 22 B2 386 144 13 35 B2 417 I59 I4 38
R3 389 154 19 25 R3 444 171 19 37. R3 509 196 20 .47
83 526 203 20 42 B3 562 213 20 50 B3 643 245 21 61 ..
R4 683 265 25 55 R4 716 275 576 63 R4 879 325 26 97

(c) Algorithm of Corvalan and Saita (1991) with N=3, dsa=O.O5, E,,= 1.0X lo-‘. 0=0.05 and max. car. =4 __ __

Determinant monitored Condition number monitored Estimate of condition number monitored

o=o.5 a=0.5 a=O.05
Point lSJd Cont Search Dct Point LUd Cont Search Cond Point LUd Cont Search Condest

S 0 0 0 s 0 0 0 0 s 0 0 0 0
RI 26 7 3 Rl 76 23 6 7 Rl’ 77 26 7 3
Bl 250 77 8 28 Bl 263 80 8 35 Bl 214 64 9 16
R2 341 103 14 35 R2 341 104 14 39 R2 314 97 I5 23
B2 422 128 15 46 B2 422 129 I5 50 B2 385 I19 16 32
R3 491 153 21 ‘53 R3 494 158 22 55 R3 458 146 23 39
B3 644 197 22 75 B3 599 1% 24 67 B3 559 II7 25 50
R4 829 256 28 97 R4 765 242 30 85 R4 729 227 31 74

(d) Algorithm III of den Heijer and Rheinboldt (1981) with b’=l.Ox lo-‘, adh=5.Ux lo-‘. cd,“= I.05 and ~~=3.0
Determinant monitored Condition number monitored Estimate of conditon number monitored
a=0.5 a=0.5 a=0.6
Point LUd Cont Search Det Point LUd Cont Search &nd Point LUd Cont Search Condest

S 0 0 0 0 s 0 0 0 0 s 0 0 0 0
Rl 95 36 8 2 Rl 77 29 6 I Rl 67 27 6 0
Bl 280 82 13 31 Bl 307 86 14 38 Bl 267 81 I4 26
R2 363 II2 21 35 R2 379 II4 21 40 R2 358 II2 21 35
B2 431 139 27 40 B2 447 I41 27 45 82 461 144 26 53
R3 492 165 35 _ 508 167 35 47 R3 525 169 33 58
B3 648 213 43 y8 ;; 647 211 43 59’ B3 716 220 39 93
R4 861 266 52 82 R4 735 250 55 59 R4’ 1436 333 47 176

Point = point located (S: I = 0. R: I= 1, B: dct A = 0)

Cont = continuation and search count
Det = determinant control count
Condest = condition number cstimatc control count
LUd = LU-decomposition count
Search = search count
Cond = condition number control count

that a = 0.667 prevents segment jumping when an condition number at each continuation point is equi-
estimate of the condition number is used. valent to that of an extra LU decomposition each
Monitoring the condition number results in a slight time. The method of estimating the condition
decrease in the number of LU decompositions number by Cline et al. (1979) makes it considerably
required but overall is much less efficient than deter- cheaper to calculate, but the uncertainty of its accur-
minant monitoring, as the cost of calculating the acy, which generally is such that 0.1 G
 Path tracking algorithm 655

K&A)/K,(A) < 1.0, gives riseto a higher LUd count, Ediv= measure of corrector divergence for den
Heijer and Rheinboldt algorithm
and a higher value of 6 necessary to prevent seg- E,~= upper bound of corrector-convergence for
ment jumping, than both the determinant and con- algorithm of Corvalan and Saita
dition number monitoring. 6 = step size parameter for algorithm of Coriralfin~
and Saita (O<fIsl)
It is proposed that determinant monitoring is the K, = condition number based on C, norm
,.
superior monitoring phase for preventing segment K~., = estrmate of K,

jumping in the continuation space, especially when K~,,. /K,,~ = maXhUm and minimum relative step size
change for den Heijer and Rheinboldt algor-
homotopy techniques are being used to solve ithm
medium to large scale problems, where the addi- rl = function of successive errors of the corrector
tional CPU time required for the other phases of steps used in algorithm of Corvdlan and Saita
u = direction factor (u = + 1)
monitoring makes them infeasible. Determinant
monitoring has also been shown to be a more Subscripts and superscripts

efficient way of preventing segment jumping than i= correction count

k= continuation count
simply tightening the parameters of the path- 0= starting point
tracking algorithm used. *= root
Another method of preventing segment jumping 11
./I2= Euclidian norm
being currently studied is the interval arithmetic
technique proposed by Kearfott (1990). For small REFERENCES
scale problems, where the cost of calculating inter-
val extensions of the Jacobian matrix is not too Allgower E. L., A survey of homotopy methods for
smooth mappings. Lecture Notes in Mathematics No.
prohibitive, this could provide an alternative
878, Numerical Solution of Nonlinear Equations (Edited
method of ensuring that the correct continuation by Allgower E. L., K. Glashoff and H.-O. Peitgen), pp.
path is followed. l-29. Springer-Verlag. New York (1981).
Allgower E. L. and K. Georg, Simplical and continuation
methods for approximating fixed points and solutions to
.
systems of equations. SIAM Reu. 22, 28-85 (1980).
Cline A. K., C. B. Moler, G. W. Stewart and _I. H.
NOMENCLATURE
Wilkinson, An estimate for the condition number of a
matrix. SIAM J. Numer. Anal. 16,363-375 (1979).
A = augmented Jacobian Corvahtn C. M. and F. A. Saita, Automatic stepsize
f = vector of equation functions control in continuation procedures. Computers them.
g = vector of starting functions Engng 15.729-739 (1991).
h = vector of real homotopy functions den Heijer C. and W. C. Rheinboldt, On steplength
N = desired number of corrector iterations algorithms for a class of continuation methods. SIAM J.
q = continuation control variable Numer. Anal. 18, 925-948 (1981).
s = arc length Kearfott R. B., Some general bifurcation techniques.
t = homotopy parameter SIAM J. Stat. Comma. 4. 52-68 (1983).
II = vector of complex homotopy functions Kearfott R. B., Inter&d arithmetic iechniques in the com-
v = tangent vector putational solution of nonlinear systems of equations.
w = continuation point Lectures in Applied Mathematics, Vol. 26, pp. 337-357
x = vector of real variables (1990).
z = vector of complex variables Morgan A. P., Solving Polynomial Systems Using
a = continuation control factor (0 <a < 1) Continuation for Scientific and Engineering Problems.
ai, a2 = contraction factors used in Allgower and Prentice-Hall, Englewood Cliffs, NJ (1987).
Georg algorithm Rheinboldt W. and J. Burkardt, Algorithm 596. A locally
adh = lower bound for value of angle of curvature parametrized continuation process. ACM Trans. Math.
of the homotopy path in den Heijer and Software 9, 236-241 (1983).
Rheinboldt algorithm Watson L. T., S. C. Billups and A. P. Morgan, Algorithm
d’= convergence tolerances (i = 1,4) for den 652 HOMPACK: a suite of codes for globally conver-
Heijer and Rheinboldt algorithm gent algorithms. ACM Trans. Math. Software 13, 281-
E = tolerance 310 (1987).