Existence Verification for Singular Zeros of Complex Nonlinear Systems.
Computational fixed point theorems can be used to automatically verify existence and uniqueness of a solution to a nonlinear system of n equations in n variables ranging within a given region of n-space. Such computations succeed, however, only when the Jacobi matrix is nonsingular everywhere in this region. However, in problems such as bifurcation problems or surface intersection problems, the Jacobi matrix can be singular, or nearly so, at the solution. For n real variables, when the Jacobi matrix is singular, tiny perturbations of the problem can result in problems either with no solution in the region, or with more than one; thus no general computational technique can prove existence and uniqueness. However, for systems of n complex variables, the multiplicity of such a solution  can be verified.  That is the subject of this paper.Such verification is possible by computing the topological degree, but such computations heretofore have required a global search on the (n-1)-dimensional boundary of an n-dimensional region. Here it is observed that preconditioning leads to a system of equations whose topological degree can be computed with a much lower-dimensional search. Formulas are given for this computation, and the special case of rank-defect one is studied, both theoretically and empirically.Verification is possible for certain subcases of the real case. That will be the subject of a companion paper.
Introduction
. Given an approximate solution -
x to a nonlinear system of
equations F
is useful in various contexts to construct bounds
around -
x in which it is proven that there exists a unique solution x # , F
continuously di#erentiable F for which the Jacobian det(F # (x #= 0 and for which
that Jacobian is well conditioned, interval computations have no trouble proving that
there is a unique solution within small boxes with x # reasonably near the center; see
[8], [16], [23]. However, if F # conditioned or singular, such computations
necessarily must fail. In the singular case, for some classes of systems F
arbitrarily small perturbations of the problem can lead to no solutions
or an even number of solutions, so multiplicity verification is not logical. In contrast,
verification is always possible if F maps C n into C n . Here, algorithms are developed
for the multiplicity of such solutions for F
The algorithms are presented in the context of solutions that lie near the real line
of complex extensions of real systems. (Such solutions arise, for example, in bifurcation
problems.) However, the algorithms can be generalized to arbitrary solutions
z # C n with z not necessarily near the real line.
Also, verification is possible for singular solutions of particular general classes of
We will cover this in a separate paper.
# Received by the editors September 10, 1999; accepted for publication (in revised form) February
21, 2000; published electronically July 19, 2000. This work was supported by National Science
Foundation grant DMS-9701540.
http://www.siam.org/journals/sinum/38-2/36107.html
Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504 (rbk@
louisiana.edu, dian@louisiana.edu).
# Institut f?r Mathematik, Universit?t Wien, Strudhofgasse 4, A-1050 Wien, Austria (neum@
cma.univie.ac.at).
SINGULAR COMPLEX ZEROS VERIFICATION 361
1.1. Previous work, related material, and references. The emphasis in
this paper is on rigorous verification of existence of a zero of a system of nonlinear
equations in a small region containing an approximate, numerically computed solu-
tion. Verification for F : R n
with the Jacobi matrix of F nonsingular at points
x with F done with computational fixed point theorems based on interval
Newton methods. Such methods are introduced, for example, in the books [2], [8],
[11], [16], [21], and [23].
The techniques in this paper for handling singularities are based on the topological
degree. Introductions to degree theory include parts of [3] (in German) or [20]. A
basic computational procedure for the degree over large regions appears in Stenger
[27]. Stynes [28], [29] and Kearfott [12], [13], [14] derived additional formulas and
algorithms based on Stenger's results. These degree computation procedures, however,
involved heuristics, and the result was not guaranteed to be correct. Aberth [1] based
a verified degree computation method on interval Newton methods and a recursive
degree-computation formula such as Theorem 2.2 below. The work here di#ers from
this previous work in two important aspects:
. The algorithms here execute in polynomial time with respect to the number
of variables and equations, 1 and
. the algorithms here assume at least second-order smoothness, and are meant
to compute the degree over small regions containing the solution, over which
certain asymptotic approximations are valid.
The treatment of verified existence represented in this paper involves computation
of the topological degree in n-dimensional complex space. In loosely related work,
Vrahatis et al. develop an algorithm for computing complex zeros of a function of a
complex variable in [31].
Finally, most of the literature we know on specialized methods for finding complex
zeros, verified or otherwise, of equations and systems of equations deals with
polynomial systems. Along these lines, continuation methods, as introduced in [6]
and [22], figure prominently. The article [4] contains methods for determining the
complex zeros of a single polynomial, while [7] and [9] contain verified methods for
determining the complex zeros of a single polynomial.
1.2. Notation. We assume familiarity with the fundamentals of interval arith-
metic; see [16, 23] for an introduction in the present context. (The works [2], [8], [24]
also contain introductory material.)
Throughout, scalars and vectors will be denoted by lower case, while matrices
will be denoted by upper case. Intervals, interval vectors (also called "boxes"), and
interval matrices will be denoted by boldface. For instance,
an interval vector, denotes an
interval matrix. Real n-space will be denoted by R n , while the set of n-dimensional
interval matrices will be denoted by IR n-n . Similarly, complex n-space will be denoted
by C n . The midpoint of an interval or interval vector x will be denoted by m(x). The
nonoriented boundary of a box x will be denoted by #x while its oriented boundary
will be denoted by b(x) (see section 2).
1.3. Traditional computational existence and uniqueness. Computational
existence and uniqueness verification rests on interval versions of Newton's method.
Typically, such computations can be described as evaluation of a related interval
1 The general degree computation problem is NP-complete; see [26].
362 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
operator implies existence and uniqueness of the solution of
To describe these, we review the following definition.
Definition 1.1 (see [23, p. 174], etc. Let F : R n
. The matrix A is said
to be a Lipschitz matrix for F over x provided for every x # x and y # x, F (x) -
A.
Most interval Newton methods for F : R n
abstractly, are of the general
where v is computed to contain the solution set to the interval linear system
and where, for initial uniqueness verification, A is generally a Lipschitz matrix 2 for F
over the box (interval vector) x and - x # x is a guess point. We sometimes write F # (x)
in place of A, since the matrix can be an interval extension of the Jacobi matrix of
F . Uniqueness verification traditionally depends on regularity of the matrix A. We
have the following lemma.
Lemma 1.2. (see [16], [23]). Suppose -
is the image under the interval
Newton method (formula (1.1)), where v is computed by any method that bounds the
solution set to the interval linear system (1.2), and -
x # x. Then A is regular.
The method of bounding the solution set of (1.2) to be considered here is the
interval Gauss-Seidel method, defined by the following definition.
Definition 1.3. The preconditioned interval Gauss-Seidel image GS(F ; x, - x) of
a box x is defined as GS(F ; x, -
x i is defined sequentially for
to n by
where
and where -
is an initial guess point, Y A # IR n-n and Y F (-x) are the
matrix and right-hand-side vector for the preconditioned interval system Y A(x-
-Y F (-x), Y # R n-n is a point preconditioning matrix, Y i denotes the ith row of Y ,
and A j denotes the jth column of A.
Lemma 1.2 applies when N(F
x), provided we specify that
x) be in the interior 3 int(x) of x. In particular, we have the following
theorem.
Theorem 1.4 (see [16], [23]). Suppose F : x # R n
A is a Lipschitz
matrix such as an interval extension F # (x) of the Jacobi matrix. If -
x is the image
under an interval Newton method as in formula (1.1) and -
x # int(x), then there is a
Various authors have proven Theorem 1.4; see [16], [23]. In particular, Miranda's
theorem can be used to easily prove Theorem 1.4 for
see
[19], [30], or [16, p. 60]. For worked-out examples, see [18, p. 3] or [17].
However, see [16, 25] for techniques for using slope matrices.
3 We must specify the interior because of the intersection step in Definition 1.3.
SINGULAR COMPLEX ZEROS VERIFICATION 363
Inclusion in the interval Gauss-Seidel method is possible because the inverse
midpoint preconditioner reduces the interval Jacobi matrix to approximately a diagonal
matrix. In the singular case, an incomplete factorization for the preconditioner
leads to an approximate diagonal matrix in the upper (n-1)- (n-1) submatrix, but
with approximate zeros in the last row. We discovered the methods in this paper by
viewing the interval Gauss-Seidel method on this submatrix, then applying special
techniques to the preconditioned nth function.
1.4. A simple singular example. Consider the following example.
Example 1. Take
and
Even though there is a unique root x
F is as in Example 1, the interval Gauss-Seidel method cannot prove this, since the
In fact, the interval Jacobi matrix is computed to
be
,
and the midpoint matrix is m(F #
). The midpoint matrix, often used
as the preconditioner Y , is singular. 4
Symbolic methods can be used to show that Example 1 has a unique solution at
arbitrarily small perturbations of the problem result in
either no solutions or two solutions. Consider the following example.
Example 2. Take
and
. Here, |#| is very small.
The system in Example 2 has two solutions for # < 0 and no solutions for # > 0.
Roundout in computer arithmetic and, perhaps, uncertainties in the system itself due
to modelling or measurement uncertainties, however, make it impossible to distinguish
systems such as in Example 2 for di#erent #, especially when computer arithmetic is
used as part of the verification process. In such instances, no verification is possi-
ble. However, if F is viewed as a complex function of two variables, then, for all #
su#ciently small, F has two solutions in a small box in C 2 containing the real point
(0, 0).
More generally, we can extend an n-dimensional box in R n to an n-dimensional
box in C n by adding a small imaginary part to each variable. If the system can
be extended to an analytic function in complex n-space (or if it can be extended to
a function that can be approximated by an analytic function), then the topological
degree gives the number of solutions, counting multiplicities, within the small region
in complex space. (See section 2 for an explanation of multiplicity.) For example,
Alternate preconditioners can nonetheless be computed; see [16]. However, it can be shown that
uniqueness cannot be proven in this case; see [16], [23].
364 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
the degree of the system in Example 2 within an extended box in complex space
can be computed to be 2, regardless of whether # is negative, positive, or zero. (See
the numerical results in section 8.) The topological degree corresponds roughly to
algebraic degree in one dimension; for example, the degree of z n in a small region in
containing 0 is n.
1.5. Organization of this paper. A review of properties of the topological
degree, to be used later, appears in section 2. The issue of preconditioning appears
in section 3. Construction of the box in the complex space appears in section 4.
Several algorithms have previously been proposed for computing the topological
degree [1], [12], [28], but these require computational e#ort equivalent to finding all
solutions to 4n (2n-1)-dimensional nonlinear systems within a given box, or worse. In
section 5, a reduction is proposed that allows computation of the topological degree
with a search in a space of dimension equal to the rank defect of the Jacobian matrix.
A theorem is proven that further simplifies the search.
In section 6, the actual algorithm is presented and its computational complexity is
given. Test problems and the test environment are described in section 7. Numerical
results appear in section 8. Future directions appear in section 9.
2. Review of some elements of degree theory. The topological degree or
Brouwer degree, well known within algebraic topology and nonlinear functional anal-
ysis, is both a generalization of the concept of a sign change of a one-dimensional
continuous function and of the winding number for analytic functions. It can be used
to generalize the concept of multiplicity of a root. The fundamentals will not be
reviewed here, but we refer to [3], [5], [12]. We present only the material we need.
Here we explain what we mean by "multiplicity." Actually, there is a more general
concept index (see [5, Chapter I]) for an isolated zero. The topological degree is equal
to the sum of the indices of zeros in the domain. The index is always positive in
our context. For this reason, we use the more suggestive term multiplicity as an
alternative term for index.
Suppose that F : D # C n
# C n is analytic. Then the real and imaginary
components of F and its argument z # C n may be viewed as real components in R 2n .
by -
R 2n by -
we have the following property of topological
degree
D, 0), and relationships between
D, 0) and the solutions of the system
Theorem 2.1 (see [5], [20], etc. Suppose F : D # C n
# C n is analytic, with
F (z) #= 0 for any z #D, and suppose -
D and -
D # R 2n are defined as above.
Then
D,
D, only if there is a solution z # D, F (z #
D, is equal to the number of solutions z # D, F (z # counting
multiplicities.
(4) If the Jacobi matrix F # (z # ) is nonsingular at every z # D with F (z #
then
D, is equal to the number of solutions z # D, F (z #
The following three theorems lead to the degree computation formula in Theorem
5.1 in section 5, the formula used in our computational scheme.
Theorem 2.2. (see [27, section 4.2]). Let D be an n-dimensional connected,
oriented region in R n and are continuous
SINGULAR COMPLEX ZEROS VERIFICATION 365
functions defined in D. Assume F #= 0 on the oriented boundary b(D) of D, b(D) can
be subdivided into a finite number of closed, connected (n - 1)-dimensional oriented
subsets # k
and there is a
on the oriented boundary b(# k
n-1 ) of
has the same sign at all solutions of there are any, on # k
Choose s # {-1, +1} and let K 0 denote the subset of the integers k # {1, . , r}
such that has solutions on # k
n-1 and sgn(f p at each of those solutions.
Then
The formula in Theorem 2.2 is a combination of formulas (4.15) and (4.16) in
[27]. The orientation of D is positive and the orientations of # k
positive
or negative, are induced by the orientation of D. If we assume that the Jacobi matrices
of F-p are nonsingular at all solutions of
depending on whether # k
n-1 has positive orientation or
negative orientation, and JF-p (x) is the determinant of the Jacobi matrix of F-p at
x. (See Theorem 5.2 and Theorem 7.2 in Chapter I of [5].) Thus we can simplify the
formula in Theorem 2.2 as follows.
Theorem 2.3. Suppose the conditions of Theorem 2.2 are satisfied and, addi-
tionally, the Jacobi matrix of F-p is nonsingular at each solution of
for each k # K 0 (s). Then
depending on whether # k
n-1 has positive orientation or
negative orientation, and JF-p (x) is the determinant of the Jacobi matrix of F-p at
x.
In our context, the region D is an n-dimensional box
The boundary #x of x consists of 2n (n - 1)-dimensional
boxes
The following theorem, necessary for the main characterization used in our algo-
rithm, is a basic property of oriented domains in n-space and follows from definitions
such as in [3]. See [18, pp. 7-8] for a detailed derivation in terms of oriented simplices.
Theorem 2.4. If x is positively oriented, then the induced orientation of x k is
and the induced orientation of x k
is
The oriented boundary b(x) can be divided into x k and x k
the associated orientations. Also, F #= 0 on b(x) is the same as F #= 0 on #x.
366 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
z }| {
Fig. 3.1. A singular system of rank n - p preconditioned with an incomplete LU factorization,
where "#" represents a nonzero element.
Now fix a p between 1 and n. Then F-p
) is the same
as F-p
. For this fixed p, let K 0 denote the subset of the
integers k # {1, . , n} such that has solutions on x k and sgn(f p at these
solutions, and let K 0 denote the subset of the integers k # {1, . , n} such that
has solutions on x k
and sgn(f p at these solutions, where s # {-1, +1}.
Then, by Theorem 2.3, we have the following theorem.
Theorem 2.5. Suppose F #= 0 on #x, and suppose there is
that
(1) F-p #= 0 on #x k or #x k ,
has the same sign at all solutions of there are any, on x k or
(3) the Jacobi matrices of F-p are nonsingular at all solutions of
Then
sgn
#F-p
x#x
sgn
#F-p
3. On preconditioning. The inverse midpoint preconditioner approximately
diagonalizes the interval Jacobi matrix when F # (and well enough
conditioned). This preconditioner can be computed with Gaussian elimination with
partial pivoting. We can compute (to within a series of row permutations) an LU
factorization of the midpoint matrix m F # (x)  . The factors L and U may then be
applied to actually precondition the interval linear system.
When the rank of F # Gaussian elimination with
full pivoting can be used to reduce F # (x) to approximately the pattern shown in

Figure

3.1. Actually, an incomplete factorization based on full pivoting will put the
system into a pattern that resembles a permutation of the columns of the pattern
in

Figure

3.1. However, for notational simplicity, there is no loss here in assuming
exactly the form in Figure 3.1.
SINGULAR COMPLEX ZEROS VERIFICATION 367
In the analysis to follow, we assume that the system has already been precondi-
tioned, so that it is, to within second-order terms with respect to w(x), of the form
in

Figure

3.1. Here we concentrate on the case p=1, although the idea can be applied
to the general case.
4. The complex setting and system form. Below, we assume
R n can be extended to an analytic function in C n .
small box that will be constructed
centered at the approximate solution -
x is near a point x # with F such that #-x - x # is much smaller than
the width of the box x, and width of the box x is small enough so that mean
value interval extensions lead, after preconditioning, to a system like Figure
3.1, with small intervals replacing the zeros.
(4) F has been preconditioned as in Figure 3.1, and F # space of
dimension 1.
The following representation is appropriate under these assumptions:
.
to complex space: x iy, with y in a small box
is centered at (0, . , 0). Define
z #
iy)), and v k (x, y) #(f k
Then, if preconditioning based on complete factorization of the midpoint matrix for
F # (x) is used, the first-order terms are eliminated in the pattern of Figure 3.1, and,
,
#xn (-x)y n +O  #(x -
,
and
(-x)y k y l +O  #(x -
,
.
5. Simplification of a degree computation procedure. To use Theorem 2.5
to compute the topological degree
directly in a verification algorithm would
require a global search of the 4n (2n-1)-dimensional faces of the 2n-dimensional box
z for zeros of -
F-p . This is an inordinate amount of work for a verification process
368 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
that would normally require only a single step of an interval Newton method in the
nonsingular case. However, if the system is preconditioned and in the form described
in section 3 and section 4, the verification can be reduced to 4n-4 interval evaluations
and four one-dimensional searches.
To describe the simplification, define
Similarly define y k and y k
. Also define
F-un
To compute the degree
F , z, 0), we will consider -
F-un on the boundary of z. The
boundary of z consists of the 4n faces x 1 , x 1
Observe that, for
F-un
#xn
|#fk /#xn (-x)| #
w(x n ). Similarly, -
F-un
F-un
F-un
on y k
implies w(y k )/|#f k /#x n (-x)| # w(y n ). Thus if x n is chosen so that
min
,
then it is unlikely that u k (x,
. Similarly, if y n is chosen so
that
min
,
then it is unlikely that v k (x, on either y k or y k . Here, the coe#cient 1
2 is
to take into consideration the fact that u k (x, y) #
#xn
#xn
(-x)y n are only approximate equalities. (When #f k /#x n
there is no restriction on w(x n ) or w(y n ) due to w(x k ) or w(y k ).)
By constructing the box z in this way, we can eliminate search of 4n - 4 of the
4n faces of the boundary of z, since we have arranged to verify -
F-un (x, y) #= 0 on
each of these faces. Elimination of these 4n - 4 faces needs only 4n - 4 interval eval-
uations. Then, we need only to search the four faces x n , x n , y n , and y n for solutions
of -
F-un regardless of how large n is. This reduces total computational cost
dramatically, since searching a face is expensive. Based on this, the following theorem
underlies our algorithm in section 6.1.
Theorem 5.1. Suppose
, and v k #= 0 on y k and y k
has a unique solution on x n and x n with y n in the interior of y n ,
and -
has a unique solution on y n and y n with x n in the interior of
(3) at the four solutions of -
condition 2; and
(4) the Jacobi matrices of -
F-un are nonsingular at the four solutions of -
in condition 2.
SINGULAR COMPLEX ZEROS VERIFICATION 369
Then
F , z,
F-un (x,y)=0
un (x,y)>0
sgn
F-un
xn=xn
F-un (x,y)=0
un (x,y)>0
sgn
F-un
F-un (x,y)=0
un (x,y)>0
sgn
F-un
F-un (x,y)=0
un (x,y)>0
sgn
F-un
Proof. Condition 1 implies -
conditions 2 and 3 imply -
F #= 0 on #z.
Condition 1 implies -
F-un #= 0 on #x k , #x k
consists of 2(n - 1) (2n - 2)-dimensional boxes, each of which is either embedded
in some x k , x k
or is embedded in y n or y n . Thus, by
conditions 2 and 3, -
F-un #= 0 on #x n . Similarly, -
F-un #= 0 on #x n , #y n and #y n .
Thus condition 1 in Theorem 2.5 is satisfied.
Condition 2 in Theorem 2.5 is automatically satisfied since either has no
solutions or a unique solution on x k , x k
Then, with condition 4, the conditions of Theorem 2.5 are satisfied. The formula
is thus obtained with
The conditions of Theorem 5.1 will be satisfied when the system is that as described
in section 3 and section 4, the box z is constructed as in (5.1) and (5.2), and
the quadratic model is accurate. (See Theorem 5.2 and its proof of the results when
all the approximations are exact.)
In Theorem 5.1, the degree consists of contributions of the four faces we search.
We can compute the degree contribution of each of the four faces, then add them to
get the degree.
In Theorem 5.1 we choose We can also choose s = -1. That doesn't
make any di#erence in our context if we ignore higher order terms in the values of u n
at the solutions of -
on the four faces x n , x n , y n , and y n . To be specific, the
four values of u n are
,
,
,
,
370 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
respectively, where # is defined in (5.3). When we choose w(y k ) the same (or roughly
the same) as w(x k ), the values of u n as a function of y n
(or y n ) will be the same (or
roughly the same) as the values of u n as a function of x n -
we ignore higher order terms, the cost of verifying u n < 0 and searching for solutions
of -
is approximately the same as the cost of verifying u n > 0
and searching for solutions of -
Next we will give a theorem that will further reduce the search cost by telling us
how we should search. Define
Theorem 5.2. If the approximations of (4.1) and (4.2) are exact, if we construct
the box z as in (5.1) and (5.2), and if #= 0, then
F , z, 2.
Proof. Under the assumptions,
Due to the construction of the box z, u
1. Next we locate
the solutions of -
(1) On x n ,
Plugging (5.8) and (5.9) into (5.6) and (5.7), we get
Then
SINGULAR COMPLEX ZEROS VERIFICATION 371
since #= 0. Thus by (5.9)
Therefore -
has a unique solution (-x, -
on x n . Plugging (5.12) into (5.10), we get the u n value at this solution, which
is
Next we compute the determinant of the Jacobi matrix of -
F-un at this solu-
tion.
Noting (5.4), (5.5), and (5.7), we have
F-un
(2) Similarly, on x n , -
has a unique solution (-x, -
y) on x n . The u n value
at this solution is
The determinant of the Jacobi matrix of -
F-un at this solution is
F-un
(3) On y n ,
Plugging (5.18) and (5.19) into (5.6) and (5.7), we get
372 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
Then
since #= 0. Thus by (5.18),
Therefore -
has a unique solution (-x,
on y n . Plugging (5.22) into (5.20), we get the u n value at this solution,
which is
Next, as in (5.15), we compute the determinant of the Jacobi matrix of -
F-un
at this solution. Noting (5.4), (5.5), and (5.7), we have
F-un
#.
Similarly, -
has a unique solution (-x, - y) on y n . The u n value at this
solution is
n .
The determinant of the Jacobi matrix of -
F-un at this solution is
F-un
Finally, we can use the formula in Theorem 5.1 to compute the topological degree
F , z, 0). If # > 0, then we know from (5.14), (5.16), (5.24), and (5.26) that
at the solutions of -
We also know the signs of
the determinants of the Jacobi matrices at the two solutions from (5.15) and (5.17).
Therefore,
F , z, 2. If # < 0, then we know from (5.14),
(5.16), (5.24), and (5.26) that u n > 0 at the solutions of -
y n . We also know the signs of the determinants of the Jacobi matrices at the two
solutions from (5.25) and (5.27). Therefore
F , z, also in this
case.
The proof of Theorem 5.2 tells us approximately where we can expect to find the
solutions of -
on the four faces we search and the value of the degree we can
expect when the approximations (4.1) and (4.2) are accurate.
From (4.1), we know that if x n is known precisely, formally solving u k (x,
for x k gives sharper bounds -
larly, if y n is known precisely, formally solving v k (x, sharper bounds
SINGULAR COMPLEX ZEROS VERIFICATION 373
n- 1. Thus when we search x n (or x n ) for
solutions of -
0, we can first get sharper bounds for x k , 1 # k # n-1, since x n is
known precisely. Then, for a small subinterval y 0
n of y n , we can solve v k (x,
y k to get sharper bounds -


Thus we get a small subface of x n (or x n ) over which we can either use an interval
Newton method to verify the existence and uniqueness of the zero of -
F-un or use
mean-value extensions to verify that -
F-un has no zeros, depending on whether y 0
n is in
the middle of y n or not. Thus the process reduces to searching over a one-dimensional
interval y n . This further reduces the search cost. We can similarly search y n or y n .
6. The algorithm and its computational complexity.
6.1. Algorithm. The algorithm consists of three phases. In the box-setting
phase, we set the box z. In the elimination phase, we verify that u k #= 0 on x k and
, and v k #= 0 on y k and y k
1. In the search phase, we verify
the unique solution of -
in the interior of y n , and on
y n and y n with x n in the interior of x n , compute the signs of u n and the signs of
the Jacobi matrices of -
F-un at the four solutions of -
compute the degree
contributions of the 4 faces x n , x n , y n , and y n according to the formula in Theorem
5.1, and finally add the contributions to get the degree.
Algorithm
Box-setting phase
1. Compute the preconditioner of the original system, using Gaussian elimination
with full pivoting.
2. Set the widths of x k and y k (see explanation below), for 1 # k # n - 1.
3. Set the widths of x n and y n as in (5.1) and (5.2).
Elimination phase
1. Do for
(a) Do for x k and x k
i. Compute the mean-value extension of u k over that face.
ii. If 0 # u k , then stop and signal failure.
(b) Do for y k and y k
i. Compute the mean-value extension of v k over that face.
ii. If 0 # v k , then stop and signal failure.
Search phase
1. Do for x n and x n
(a) i. Use mean-value extensions for u k (x, to solve for x k to get
sharper bounds -
x k with width O  #(x - x, y)# 2
ii. If -
return the degree contribution of that face as
iii. Update x k .
(b) i. Compute the mean-value extension u n over that face.
ii. If u n < 0, then return the degree contribution of that face as 0.
(c) Construct a small subinterval y 0
n of y n which is centered at 0.
(d) i. Use mean-value extensions for v k (x, to solve for y k to get
sharper bounds -
y k with width O  max(#(x -


thus getting a subface x 0
n ) of x n (or x n .)
ii. If -
#, then stop and signal failure.
374 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
up an interval Newton method for -
F-un to verify existence and
uniqueness of a zero in the subface x 0
ii. If the zero cannot be verified, then stop and signal failure.
(f) Inflate y 0
n as much as possible subject to verification of existence and
uniqueness of the zero of -
F-un over the corresponding subface, and thus
get a subinterval y 1
n of y n .
(g) In this step, we verify -
solutions when y
n .
n has two separate parts; we denote the lower part by y l
n and the
upper part by y u
n . We present only the processing of the lower part.
The upper part can be processed similarly.
A. Use mean-value extensions for v k (x, to solve for y k to
get sharper bounds for y k , 1 # k # n - 1, and thus to get a
subface of x n (or x n ).
B. Compute the mean-value extensions -
F -un over the subface obtained
in the last step.
F -un , then bisect y l
update the lower part as a new y l
and cycle.
F -un , then exit the loop.
ii. Do
A. If y 1
exit the loop.
B. y l
C. Use mean-value extensions for v k (x, to solve for y k to
get sharper bounds for y k , 1 # k # n - 1, and thus to get a
subface of x n (or x n ).
D. Compute the mean-value extensions -
F -un over the subface obtained
in the last step.
F -un , then cycle.
F -un , then y l
, mid(y l
Compute the mean-value extension of u n over x 0
.)
ii. If u n < 0, then return the degree contribution of that face as 0.
F-un
F-un
ii. If 0 # -
F-un
F-un
stop and signal failure.
(j) Use the formula in Theorem 5.1 to compute the degree contribution of
that face.
2. Do for y n and y n
(a) Same as step 1(a) except change x k to y k , -
x k to -
to v k .
(b) Same as step 1(b).
(c) Same as step 1(c) except change y 0
n to x 0
n , y n to x n , and 0 to -
x n .
(d) Same as step 1(d) except change y k to x k , -
y k to -
n to y 0
n to y 0
to y n , and x n to y n .
Same as step 1(e) except change x 0
n to y 0
n and x 0
n to y 0
n .
(f) Same as step 1(f) except change y 0
n to x 0
n to x 1
n , and y n to x n .
(g) Same as step 1(g) except change y n \ y 1
n to x n \ x 1
n .
Same as step 1(h) except change x 0
n to y 0
n and x 0
n to y 0
SINGULAR COMPLEX ZEROS VERIFICATION 375
(i) Same as step 1(i) except change
F-un
F-un
F-un
F-un
(j) Same as step 1(j).
3. Add the degree contributions of the four faces obtained in steps 1 and 2 to
get the degree.
End of algorithm
An explanation of the algorithm
1. In the box-setting phase, in step 2, the width w(x k ) of x k depends on the
accuracy of the approximate solution -
x of the system F should
be much larger than |-x k - x # k |. At the same time, w(x k ) should not be too
large, since the quadratic model needs to be accurate over the box.
2. In the search phase, in step 1(b) (or 2(b)), we check the sign of u n on that
face and discard that face at the earliest possible time if u n < 0 on that
face, since we know the degree contribution of that face is 0 according to the
formula in Theorem 5.1. This will save time significantly if it happens that
on that face. It did happen for all the test problems. (See section 8
for the test results.)
3. In the search phase, in step 1(e) (or 2(e)), we precondition the system -
F-un
before we use an interval Newton method, so that the method will succeed
(see section 1.3 and section 3). The system -
F-un is nonsingular over the
subfaces under consideration.
4. In the search phase, in step 1(f) (or 2(f)), we first expand the subinterval
n  at both ends. If existence and uniqueness of the
zero of -
F-un can be verified over the corresponding subface, then we expand
the subinterval by 2# at both ends, then 4# and so on until existence and
uniqueness verification fails.
5. In the search phase, in step 1(g) (or 2(g)), the underlying idea is that the
farther away the interval y l
n is from the interval y 0
whose corresponding
subface of x n (or x n ) contains a unique solution of -
or the narrower
the interval y l
is, the more probable it is that we can verify that -
F-un
over the subface of x n (or x n ) corresponding to y l
n .
6.2. Computational complexity.
Derivation of the computational complexity
Box-setting phase: Step 1 is of order O n 3
. Step 2 is of order O (n). Step 3 is of
order O n 2
. Thus, the order of this phase is O n 3
.
Elimination phase: Step 1(a)i and 1(b)i are of order O n 2
. Step 1(a)ii and 1(b)ii
are of order O (1). Thus, the order of this phase is O n 3
.
Search phase: Step 1(a) and 2(a) are of order O n 3
. Step 1(b) and 2(b) are of
order O n 2
. Step 1(c) and 2(c) are of order O (1). Step 1(d) and 2(d) are
of order O n 3
. Step 1(e) and 2(e) are of order O n 3
. Step 1(f) and 2(f)
are of order N inf l *O n 3
. (See explanation below.) Step 1(g) and 2(g) are
of order N proc *O n 3
. (See explanation below.) Step 1(h) and 2(h) are of
order O n 2
. Step 1(i) and 2(i) are of order O n 3
. Step 1(j) and 2(j) are
of order O (1) . The last step of this phase is of order O (1) too. Thus, the
order of this phase is O n 3
376 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
The order of the overall algorithm is thus O n 3
.
Remark. The order of the algorithm cannot be improved, since computing preconditioners
of the original system and the system -
F-un is necessary and computing
each preconditioner is of order O n 3
.
7. Test problems and test environment.
7.1. Test problems. Before describing the test set, we introduce one more
problem. Motivated by [10, Lemma 2.4], we considered systems of the following form.
Example 3. Set
the matrix corresponding to central di#erence discretization of the boundary value
problem
. The parameter t was
chosen to be equal to t is the largest eigenvalue of A.
The homotopy h in Example 3 has a simple bifurcation point at where the
two paths cross obliquely. That is, there are two solutions to
all t near t 1 and on either side of t 1 . Furthermore, the quadratic terms in the Taylor
expansion for f do not vanish at
The test set consists of Example 1, Example 2 with
and Example 3 with For all the test problems, we used
(0, 0, . , 0) as a good approximate solution to the problem F Actually, it's the
exact solution in Example 1 and Example 3. w(x k ) and w(y k ) were set to 10 -3 for 1 #
computed automatically by the algorithm. In fact,
can also be computed automatically by the algorithm,
depending on the accuracy of the approximate solution. At present, we used the known
true solutions to Example 1 and Example 3 and the known approximate solution to
Example 2 to test the algorithm and set the widths apparently small but otherwise
arbitrary.
For all the problems, the algorithm succeeded and returned a degree of 2.
7.2. Test environment. The algorithm in section 6.1 was programmed in the
Fortran 90 environment developed and described in [15], [16]. Similarly, all the functions
in the test problems were programmed using the same Fortran 90 system, and
internal symbolic representations of the functions were generated prior to execution of
the numerical tests. In the actual tests, generic routines then interpreted the internal
representations to obtain both floating point and internal values.
The LINPACK routines DGECO and DGESL were used in step 1 of the box-setting
phase, and in step 1(e), 2(e), 1(f), and 2(f) of the search phase to compute the
preconditioners. (See the algorithm and its explanation in section 6.1.)
The Sun Fortran 90 compiler version 1.2 was used on a Sparc Ultra model 140
with optimization level 0. Execution times were measured using the routine DSECND.
All times are given in CPU seconds.
8. Numerical results. We present the numerical results in Table 8.1 and some
statistical data in Table 8.2.
The column labels of Table 8.1 are as follows.
Problem: names of the problems identified in section 7.1.
n: number of independent variables.
Success: whether the algorithm was successful.
Degree: topological degree returned by the algorithm.
CPU time: CPU time in seconds of the algorithm.
Time ratio: This applies only to Example 3. It's the ratio of two successive CPU
times.
SINGULAR COMPLEX ZEROS VERIFICATION 377

Table
Numerical results.
Problem n Success Degree CPU time Time ratio
Example 2
Example 2
Example 3 5
Example 3 20
Example 3
Example 3 160
Example 3 320

Table
Statistical data.
Problem
Example
Example 2
Example 2
Example
Example
Example 3
Example 3 160
Example 3
The column labels of Table 8.2 are as follows.
Problem: names of the problems identified in section 7.1.
n: number of independent variables.
number of inflations the algorithm did in step 1(f) or 2(f) for the indicated
face x n , x n , y n , or x y .
number of subintervals of y n \ y 1
n the algorithm processed in step 1(g) or
subintervals of x n \ x 1
n the algorithm processed in step 2(g), i.e., the number
of y l
's plus number of y u
n 's in step 1(g) or number of x l
's plus number of
x u
n 's in step 2(g) for the indicated face x n , x n , y n , or x y .
We can see from Table 8.1 that the algorithm was successful on each problem in the
test set. The overall algorithm is O n 3
, but the are many O n 3
and O n 2
steps.
Some steps have many O n 3
and O n 2
substeps, and some of the substeps still
have many O n 2
structures. Thus, when n was small, those lower order structures
had significant influence on the CPU time. However, for the larger n in the examples
tried, the O n 3
terms dominated. We can see this from the time ratios of Example 3
in

Table

8.1.
In

Table

8.2, in each problem there were two faces of x n , x n , y n , and y n for
which N inf l = 0. This is because the algorithm verified that u n < 0 on each of those
two faces in step 1(b) or 2(b), and returned a degree contribution of each of those
378 R. BAKER KEARFOTT, JIANWEI DIAN, AND A. NEUMAIER
two faces as 0. Thus, the algorithm didn't proceed to step 1(f) or 2(f). For the same
reason, those two faces. For the remaining two faces for which the
algorithm did proceed to step 1(f) or 2(f), N inf l is small.
In step 1(g) or 2(g), which immediately follows the inflations, N
Example 1 and Example 2. This is because the inflations had covered the whole
interval y n . More significant is that N proc = 2 in Example 3 regardless of small n
or large n. This is because only one interval was processed to verify that -
has no solutions when x n # x l
n and only one interval was processed to verify that
solutions when x
n . This means that the algorithm was quite
e#cient.
9. Conclusions and future work. When we tested the algorithm, we took
advantage of knowing the true solutions (see section 7.1. For this reason, we set
arbitrarily. But we plan to have the algorithm
eventually compute these, based on the accuracy of the approximate solution
obtained by a floating point algorithm and the accuracy of the quadratic model.
We presented an algorithm which was designed to work for the case that the rank
deficiency of the Jacobian matrix at the singular solution is one. But the analysis in
section 5 and the algorithm in section 6.1 can be generalized to general rank deficiency.
Also, at present, it is assumed that the second derivatives # 2 fn
#xk#x l
don't vanish simultaneously at the singular solution. In fact, the analysis in section
5 and the algorithm in section 6.1 can be generalized to the general case that the
derivatives of f n of order 1 through r (r # 2) vanish simultaneously at the singular
solution. Computing higher order derivatives, however, may be expensive. Those
two generalizations can also be combined, i.e., any rank deficiency and any order of
derivatives of f n that vanish. We will pursue these generalizations in the future.
Modification of the algorithm to verify complex roots that are not lying near the
real axis is possible.
Another future direction of this study is to apply the algorithms to bifurcation
problems and other physical models.
Finally, verification is possible
in a multidimensional analogue
of odd-multiplicity roots. We are presently writing up theoretical and experimental
results for this situation.



--R

Computation of topological degree using interval arithmetic
Introduction to Interval Computations
New York
Direkte Verfahren zur Berechnung der Nullstellen von Polynomen
Fixed Points and Topological Degree in Nonlinear Analysis
Fixed Points
Circular arithmetic and the determination of polynomial zeros
Global Optimization Using Interval Analysis
Applied and Computational Complex Analysis.


Computing the Degree of Maps and a Generalized Method of Bisection

A summary of recent experiments to compute the topological degree
A Fortran 90 environment for research and prototyping of enclosure algorithms for nonlinear equations and global optimization
Continuous Problems
Rigorous global optimization and the GlobSol package.
Existence Verification for Singular Zeros of Nonlinear Systems
The Poincar-e-Miranda theorem
Degree Theory
Methods and Applications of Interval Analysis
Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems
Interval Methods for Systems of Equations
New Computer Methods for Global Optimization
Verification methods for dense and sparse systems of equations
Optimal Solution of Nonlinear Equations
An algorithm for the topological degree of a mapping in R n
An Algorithm for the Numerical Calculation of the Degree of a Mapping
An algorithm for numerical calculation of the topological degree
A short proof and a generalization of Miranda's existence theorem
The topological degree theory for the location and computation of complex zeros of Bessel functions
--TR

--CTR
Jianwei Dian , R. Baker Kearfott, Existence verification for singular and nonsmooth zeros of real nonlinear systems, Mathematics of Computation, v.72 n.242, p.757-766, 1 April
R. Baker Kearfott , Jianwei Dian, Verifying topological indices for higher-order rank deficiencies, Journal of Complexity, v.18 n.2, p.589-611, June 2002
B. Mourrain , M. N. Vrahatis , J. C. Yakoubsohn, On the complexity of isolating real roots and computing with certainty the topological degree, Journal of Complexity, v.18 n.2, p.612-640, June 2002
interval computations;complex nonlinear systems;topological degree;verified computations;singularities