Verifying Topological Indices for Higher-Order Rank Deficiencies

R. Baker Kearfott and Jianwei Dian Department of Mathematics University of Louisiana atLafayette Lafayette, LA 70504

               rbk@louisiana.edu

October 27, 2000

Abstract

It has been known how to use computational fixed point theorems to verify existence and uniqueness of a true solution to a nonlinear system of equations within a small region about an approximate solution. This can be done in O (n^{3}) operations, where n is the number of equations and unknowns. However, these standard techniques are only valid if the Jacobi matrix for the system is nonsingular at the solution. In previous work and a dissertation (of Dian), we have shown, both theoretically and practically, that existence and multiplicity can be verified in a complex setting, and in the real setting for odd multiplicity, when the rank defect of the Jacobi matrix at an isolated solution is 1. Here, after a brief introduction, we discuss the case of higher rank defect.

1 Background

Given a system of nonlinear equations, numerical methods can typically produce an approximation x̌ to a solution x^{∗}. It is then sometimes desirable to compute bounds

x = (x_{1}, x_{2}, …, x_{n})
= ([\underline{x}_{1}, \overline{x}_{1}], [\underline{x}_{2}, \overline{x}_{2}], …, [\underline{x}_{n}, \overline{x}_{n}],

such that x̌ ∈ x, and such that a computational fixed point theorem can verify that there is a true solution of the nonlinear system within x. Specifically, we examine the problem

Given F : x → ℝ^{n} and x ∈ 𝕀ℝ^{n}, rigorously verify:
• there exists a unique x^{∗} ∈ x such that F (x^{∗}) = 0.

(1)

If the Jacobi matrix F′(x^{∗}) is nonsingular, if x̌ is a sufficiently accurate approximation to x^{∗}, and if the dimensions of x are chosen appropriately, then interval Newton methods will do (1); there are numerous explanations of how this is done, such as those found in [3], [7], [8], and [5, §1.5]. These interval Newton methods are based on interval versions of linear algebra algorithms, such as the Gauss–Seidel method; such interval Newton methods achieve the verification in O (n^{3}) operations.

Common thinking has been that (1) cannot be done when F′(x^{∗}) is singular or excessively ill-conditioned. However, as explained in [6] and [2], if ℝ^{n} in (1) is replaced by ℂ^{n}, then, in principle, existence and uniqueness can still be verified. The steps of the algorithm for this singular-case verification are outwardly similar to the non-singular case, except that there is an extra low-dimensional search. In [6], we exhibited algorithms for the rank-defect-1 case, i.e. when the null space of F (x^{∗}) has dimension 1; we showed theoretically that these algorithms can verify existence and uniqueness in O (n^{3}) operations, and we illustrated this dependency on dimension with actual computations on a discretization of a model nonlinear eigenvalue problem, with dimension ranging from 2 to 320.

Our algorithms are based on rigorous computation the topological index of the map F over a box x of appropriate size centered at the approximate solution x̌; see [6] and [2] for a review and references. The algorithm in [6] was specific to the case where this topological index was 2, although we subsequently discovered that the algorithm easily generalizes to arbitrary index. This generalization is explained in [2]; the algorithm is also an O (n^{3}) algorithm.

The dissertation [2] also contains a theoretical study and an algorithm dealing directly with F : x → ℝ^{n}, where x ⊂ ℝ^{n}, rather than dealing with a complex extension. In particular, a heuristic is effective at guessing the topological index of the complex extension at a solution with a rank-1 singularity; if this topological index happens to be odd, then the topological index in ℝ^{n} can be verified to be either 1 or -1, much more efficiently (but still with O (n^{3}) operations) than the corresponding verification in the complex extension. This real-space verification has the additional theoretical advantage that an actual solution to F (x) = 0 has been verified to exist within x ⊂ ℝ^{n}, while the complex computations only verify that a solution, possibly with imaginary components, exists within a small region in complex space containing x ⊂ ℝ^{n}.

At present, we are completing implementations of the techniques for arbitrary topological index and for verifying the topological index in real space; we expect this to be straightforward and successful. However, interesting development remains for the case where the null space of F (x^{∗}) is greater than one. Here, we present this case, pointing out opportunities and difficulties. In §2, we explain our general framework for computing the topological index, while we discuss the higher-dimensional null-space case in §3.

2 Topological Index Computations: The General Setting and the Rank 1 Defect Case

Our computations are based on

1. preconditioning the system F (x) = 0 by multiplying by a constant matrix Y so the Jacobi matrix for Y F (x^{∗}) is approximately diagonal, except in p rows, where the dimension of the null space of F′(x^{∗}) is p;

2. constructing a box x, with astutely chosen coordinate widths, centered at the approximate solution x̌;

3. computing the Brouwer degree of Y F, and hence of F, over x by searching the (n − 1)-dimensional sides of x to verify solutions of a certain system of equations derived from the components of Y F.

Details of these ideas, as well as a review of properties of the Brouwer degree and references to comprehensive introductions, appear in [6] and [2]. Of interest here is the fact that, because of the form of the preconditioned system, the search on the boundary can be greatly streamlined. In particular, previous general algorithms for the topological degree, such as the heuristic algorithms in [10] or [4], as well as the rigorous algorithm in [1], have running times that depend exponentially on n. In contrast, in the preconditioned system Y F, one can, in effect, express (n − p) of the variables in terms of p variables. When p = 1, if the box dimensions are chosen appropriately, all but four of the 4n sides of the box in ℂ^{n} (treated as a box in ℝ^{2n}) may be eliminated with simple interval evaluations, and the remaining four (n − 1)-dimensional sides may be handled with onedimensional searches. Easily-obtainable approximations to the solutions of the system derived from Y F further facilitate these one-dimensional searches. Here, we present those details of that process relevant to studying generalization to p > 1.

Following [6] and [2], we observe that, if the rank defect of F′(x^{∗}) is p, then the preconditioner Y can be formed as one would compute an inverse of F′(x^{∗}), except with an incomplete LU-factorization based on full pivoting. The resulting preconditioned Jacobi matrix, to within a column permutation, has the form in Figure 1. Hence, if we assume F (z) = (f_{1}(z_{1}, …, z_{n}), …, f_{n}(z_{1}, …, z_{n}) has already been so preconditioned, then, for the rank-1 defect case p = 1, the components of f have the form

f_{k}(z) = (z_{k} − x_{k}^{∗}) + \frac{∂f_{k}}{∂z_{n}} (x^{∗})(z_{n} − x_{n}^{∗}) + O (∥z − x^{∗}∥^{2})  (2)

for 1 ≤ k ≤ n − 1,

f_{n}(z) = \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x^{∗})(z_{k_{1}} − x_{k_{1}}^{∗})(z_{k_{2}} − x_{k_{2}}^{∗}) + … + (3)

\frac{1}{d!} ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{d}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x^{∗})(z_{k_{1}} − x_{k_{1}}^{∗}) … (z_{k_{d}} − x_{k_{d}}^{∗})

Y F (x^{∗}) = ( \begin{matrix}1 & 0 & … & 0 & \overbrace{∗ … ∗}^{p} \\ 0 & 1 & 0 … & 0 & ∗ … ∗ \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & … & 0 & 1 & ∗ … ∗ \\ 0 & … & 0 & 0 & 0 … 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & … & 0 & 0 & 0 … 0 \end{matrix} ) .

Figure 1: A singular system of rank n − p preconditioned with an incomplete LU factorization, where “∗” represents a non-zero element.

+O (∥z − x^{∗}∥^{d+1}) .

(See [2].) If the actual solution x^{∗} is sufficiently close to the approximate solution x̌, then x^{∗} may be replaced by x̌ in the above equations. Furthermore, if we use the notation z_{k} = x_{k} + iy_{k}, x = (x_{1}, …, x_{n}), y = (y_{1}, …, y_{n}), and f_{k}(z) = u_{k}(x, y) + iv_{k}(x, y) for 1 ≤ i ≤ n, then the above forms for f_{k} and f_{n} may be expressed as

u_{k}(x, y) = (x_{k} − x_{k}^{∗}) + \frac{∂f_{k}}{∂x_{n}} (x^{∗})(x_{n} − x_{n}^{∗}) + O (∥(x − x^{∗}, y)∥^{2}) ,(4)
v_{k}(x, y) = y_{k} + \frac{∂f_{k}}{∂x_{n}} (x^{∗})y_{n} + O (∥(x − x^{∗}, y)∥^{2}) ,  (5)

We further define

\begin{matrix}x\underline{k} ≡ (x_{1}, y_{1}, …, x_{k−1}, y_{k−1}, \underline{x}_{k}, y_{k}, x_{k+1}, y_{k+1}, …, x_{n}, y_{n}), and \\ x\overline{k} ≡ (x_{1}, y_{1}, …, x_{k−1}, y_{k−1}, \overline{x}_{k}, y_{k}, x_{k+1}, y_{k+1}, …, x_{n}, y_{n}), \end{matrix}

and define y\underline{k} and y\overline{k} similarly. With this, define F̃ : D̃ ⊂ ℝ^{2n} → ℝ^{2n} by F̃ = (u_{1}, v_{1}, …, u_{n}, v_{n}), and define

F̃_{¬u_{n}}(x, y) ≡ (u_{1}(x, y), v_{1}(x, y), …, u_{n−1}(x, y), v_{n−1}(x, y), v_{n}(x, y)).

To compute the degree d(F̃, z, 0), we will consider F̃_{¬u_{n}} on the boundary of z. The boundary of z consists of the 4n faces x\underline{1}, x\overline{1}, y\underline{1}, y\overline{1}, …, x\underline{n}, x\overline{n}, y\underline{n}, y\overline{n}. The box coordinates x_{k}, y_{k}, 1 ≤ k ≤ n are chosen so x_{k} is centered on x̌_{k} and y_{k} is centered on 0, 1 ≤ k ≤ n, and so the widths w(x_{n}) and w(y_{n}) obey

w(x_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{ \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} \} ,  (6)
w(y_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{ \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|} \} .  (7)

Provided F̃_{¬u_{n}}(x, y) ≠ 0 for (x, y) ∈ x\underline{k}, x\overline{k}y\underline{k}, y\overline{k}, 1 ≤ k ≤ n − 1, the degree d(F̃, z, 0) may be computed with the formula

d(F̃, z, 0) = − ∑_{xn=\underline{x}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn | \frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}y_{n}} (x, y)| (8)

+ ∑_{xn=\overline{x}n ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn | \frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}y_{n}} (x, y)|

+ ∑_{yn=\underline{y}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn | \frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}x_{n}} (x, y)|

− ∑_{yn=\overline{y}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn | \frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}x_{n}} (x, y)| .

If the coordinate extents are chosen as in (6), then (4) makes it unlikely that u_{k} = 0 on x\underline{k} and x\overline{k}, 1 ≤ k ≤ n − 1; similarly, if the coordinate extents are chosen as in (7), then (5) makes it unlikely that v_{k} = 0 on y\underline{k} and y\overline{k}, 1 ≤ k ≤ n − 1. The verification algorithm verifies these facts by computing the interval values u(x\underline{k}), u(x\overline{k}), v(y\underline{k}), and v(y\overline{k}). Formula (8) is then used by systematic search (made rigorous with interval computations) of the four faces x\underline{n}, x\overline{n}, y\underline{n}, and y\overline{n} for solutions to F̃_{¬u_{n}} = 0. The search is reduced to a one-dimensional search over the y_{n} coordinate on x\underline{n} and x\overline{n} and a one-dimensional search over the x_{n} coordinate on y\underline{n} and y\overline{n} by formally solving the rigorous interval enclosure for v_{k} corresponding to (7) for y_{k} in terms of y_{n} and formally solving the rigorous interval enclosure for u_{k} corresponding to (6) for x_{k} in terms of x_{n}. For details, see [6] and [2].

The one-dimensional search is facilitated with accurate a priori approximations to the solutions of F̃_{¬u_{n}}(x, y) = 0 on x\underline{n}, x\overline{n}, y\underline{n}, and y\overline{n}. Denote

α_{k} ≡ \frac{∂f_{k}}{∂x_{n}} (x̌),  1 ≤ k ≤ n − 1,
α_{n} ≡ −1,
Δ_{1} ≡ | \frac{∂F}{∂x_{1} … ∂x_{n}} (x̌)| .
Δ_{d} ≡ ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{d}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x̌)α_{k_{1}} … α_{k_{d}},  2 ≤ d. (9)

Then, consider f_{n}(z) = (u_{n}(x, y), v_{n}(x, y)), assume the formulas (2) and (3) to be exact without the Landau symbols, and assume the actual topological index is d. In that case, all terms in (3) except terms of order d and higher vanish. Further assuming f_{k}(z) = 0, 1 ≤ k ≤ n − 1, solving

[figure: v_{n} = 0; x_{n}; y_{n}; v_{n} = 0; v_{n} = 0; v_{n} = 0; v_{n} = 0; v_{n} = 0; v_{n} = 0; v_{n} = 0; x\underline{n}; x\overline{n}; y\underline{n}; y\overline{n}]

Figure 2: When d = 4. v_{n} = 0 on solid lines and u_{n} = 0 on dashed lines. The thick dots represent the solutions of F̃_{¬u_{n}}(x, y) = 0 on the boundary.

for z_{k} in terms of z_{n} in (2), and plugging into (3) finally gives

f_{n}(z) = \frac{(−1)^{d}t^{d}Δ_{d}}{d!} (z_{n} − x_{n}^{∗})^{d} + O (∥z − x^{∗}∥^{d+1}) .  (10)

It follows from (10) that solutions of v_{n} = ℑf_{n} = 0 at those complex values z_{n} = (x_{n}, y_{n}) whose argument θ obeys θ^{d} = π/2 + ℓπ for some integer ℓ. This is illustrated for d = 4 in Figure 2. For details, see [2].

Thus, in the one-dimensional rank defect case, one-dimensional subintervals of four one-dimensional intervals can be constructed about approximate solutions of F̃_{¬u_{n}}(x, y) = 0, these intervals can be rapidly verified to contain unique solutions of F̃_{¬u_{n}}(x, y) = 0, and the remainder of the four original one-dimensional intervals can be rapidly but rigorously eliminated with interval evaluations of F̃_{¬u_{n}}0. Actual algorithms appear in [2], numerical results for d = 2 appear in [6], and we are presently completing numerical experiments for d > 2.

3 The Higher Rank-Defect Case

When the dimension p of the null space of F′(x^{∗}) is greater than 1, the forms corresponding to (2) and (3) are

f_{k}(z) = (z_{k} − x_{k}^{∗}) + \frac{∂f_{k}}{∂z_{n}} (x^{∗})(z_{n} − x_{n}^{∗}) + O (∥z − x^{∗}∥^{2})  (11)

for 1 ≤ k ≤ n − p,

f_{q}(z) = \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{q}}{∂x_{k_{1}}∂x_{k_{2}}} (x^{∗})(z_{k_{1}} − x_{k_{1}}^{∗})(z_{k_{2}} − x_{k_{2}}^{∗}) + … +

\frac{1}{d!} ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{d}f_{q}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x^{∗})(z_{k_{1}} − x_{k_{1}}^{∗}) … (z_{k_{d}} − x_{k_{d}}^{∗})

+O (∥z − x^{∗}∥^{d+1}) ,  for n − p + 1 ≤ q ≤ n. (12)

In this more general setting, (11) can be used as before to eliminate variables from (12). However, p variables remain, and there are p equations left. In general, this system is an arbitrary system of p homogeneous degree-d equations in p variables; to see this, let n = p and specify the complete original system by

f_{q}(z) = \frac{1}{d!} ∑_{k_{1}=1}^{p} … ∑_{k_{d}=1}^{p} \frac{∂^{d}f_{q}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x^{∗})(z_{k_{1}} − x_{k_{1}}^{∗}) … (z_{k_{d}} − x_{k_{d}}^{∗})
for 1 ≤ q ≤ p,  (13)

where the partial derivatives are set arbitrarily, subject only to the condition that corresponding mixed partial derivatives are equal. This implies that, in the analogue of the case when p = 1, a p-dimensional space must be searched. Furthermore, for approximate starting solutions, instead of a simple formula as for (10) and Figure 2, all solutions to a general d-homogeneous system of p equations in p unknowns would need to be found. For higher p and d, that could be expensive for a verification step that may be a small part of another overall algorithm.

The general formula from which (8) was derived is [6, Theorem 2.5]). A straightforward change of notation to our complex setting, selection of s = 1, and selection of p = 2n − 1 (corresponding to the component u_{n}) in that theorem gives the general formula

d(F, z, 0) = d(F̃, (x, y), 0)  (14)
= − ∑_{k∈\underline{K_{0}}} ∑_{{z∈x\underline{k} \atop F¬u_{n} (z)=0}} sgn(D_{1})

+ ∑_{k∈\overline{K_{0}}} ∑_{{z∈x\overline{k} \atop F¬u_{n} (z)=0}} sgn(D_{1})

+ ∑_{k∈\underline{K_{1}}} ∑_{{z∈y\underline{k} \atop F¬u_{n} (z)=0}} sgn(D_{2})

− ∑_{k∈\overline{K_{1}}} ∑_{{z∈y\overline{k} \atop F¬u_{n} (z)=0}} sgn(D_{2}),

where

D_{1} = | \frac{∂F_{¬u_{n}}}{∂x_{1}y_{1}x_{2}y_{2} … x_{k−1}y_{k−1}y_{k}x_{k+1}y_{k+1} … x_{n}y_{n}} (z)| and

D_{2} = | \frac{∂F_{¬u_{n}}}{∂x_{1}y_{1}x_{2}y_{2} … x_{k−1}y_{k−1}x_{k}x_{k+1}y_{k+1} … x_{n}y_{n}} (z)| ,

and where, following [6, p. 366], \underline{K_{0}} is that subset of the integers k ∈ \{1, …, n\} such that F_{¬u_{n}} = 0 has solutions on x\underline{k} and sgn(u_{n}) = 1 at these solutions, and \overline{K_{0}} is that subset of the integers k ∈ \{1, …, n\} such that F_{¬u_{n}} = 0 has solutions on x\overline{k} and sgn(u_{n}) = 1 at these solutions, \underline{K_{1}} is that subset of the integers k ∈ \{1, …, n\} such that F_{¬u_{n}} = 0 has solutions on y\underline{k} and sgn(u_{n}) = 1 at these solutions, and \overline{K_{1}} is that subset of the integers k ∈ \{1, …, n\} such that F_{¬u_{n}} = 0 has solutions on y\overline{k} and sgn(u_{n}) = 1 at these solutions (Note that now, in contrast to (8), any u_{q}, 1 ≤ n − p + 1 ≤ n may be chosen for F_{¬u_{q}}. For simplicity, we continue to choose q = n without loss of generality.)

Now, as in the rank-1 defect case p = 1, the box z can be constructed so u_{k} is probably nonzero on x\underline{k} and x\overline{k}, and v_{k} is probably nonzero on y\underline{k} and y\overline{k} for 1 ≤ k ≤ n − p. In particular, the forms corresponding to (4) and (5) are

u_{k}(x, y) = (x_{k} − x_{k}^{∗}) + ∑_{q=n−p+1}^{n} \frac{∂f_{k}}{∂x_{q}} (x^{∗})(x_{q} − x_{q}^{∗})  (15)
+O (∥(x − x^{∗}, y)∥^{2}) ,
v_{k}(x, y) = y_{k} + ∑_{q=n−p+1}^{n} \frac{∂f_{k}}{∂x_{q}} (x^{∗})y_{q} + O (∥(x − x^{∗}, y)∥^{2}) , (16)

from which it follows that conditions corresponding to (6) and (7) are

∑_{q=n−p+1}^{n} | \frac{∂f_{k}}{∂x_{q}} (x^{∗})| w(x_{q}) ≤ \frac{1}{2} w(x_{k}), 1 ≤ k ≤ n − p, (17)

∑_{q=n−p+1}^{n} | \frac{∂f_{k}}{∂x_{q}} (x^{∗})| w(y_{q}) ≤ \frac{1}{2} w(y_{k}), 1 ≤ k ≤ n − p. (18)

Therefore, the 4p faces x\underline{q} and x\overline{q}, y\underline{q}, and y\overline{q} for n−p+1 ≤ q ≤ n cannot be eliminated from consideration in the sums in (14). In fact, instead of being reduced to (8), (14) can only be reduced to

d(F, z, 0) = ∑_{q=1}^{p}  (19)

\{− ∑_{xn=\underline{x}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn(D_{1}) + ∑_{xn=\overline{x}n ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn(D_{1})

+ ∑_{yn=\underline{y}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn(D_{2}) − ∑_{yn=\overline{y}_{n} ; {F̃¬u_{n} (x,y)=0 \atop un(x,y)>0}} sgn(D_{2})\} .

Thus, instead of four one-dimensional intervals to be searched, 4p 2pdimensional boxes (corresponding to real and imaginary coordinates of the p variables that could not be eliminated) would need to be searched. All in all, it thus appears that the task increases exponentially with the rank defect p, and the complication of implementation jumps substantially from p = 1 to p = 2. However, it is conceivable that some kind of preconditioning based on higher-order, in addition to first derivatives, could simplify the process.

4 Uses and Limitations

Verification of topological indices is potentially useful in automatic theorem proving associated with bifurcation theory and practical bifurcation problems. It also could be useful in branch and bound optimization algorithms as described in [5, Ch. 5] or [3], to verify feasibility of a set of constraints that happen to be linearly dependent on isolated parts of the feasible set. Although such linear dependencies appear unlikely (or impossible) from a probabilistic point of view, they do occur in practice. Higher-order rank deficiencies also occur in practice. However, the difficulty of verification appears to increase rapidly with the dimension of the null space, and the implementation becomes significantly more complicated between a one-dimensional and a two-dimensional null space.

In any case, the techniques treated here are in general only applicable to isolated solutions. To see this, note that a condition that the topological degree d(F, x, 0) be defined is that there be no solutions to F (x) = 0 on the boundary of x; if a solution x^{∗} to F (x) = 0 is not isolated, then x in general cannot be so chosen.

An algorithm for exhaustively analyzing the solution sets of polynomial systems that have higher-dimensional solution sets is described in [9]. Although that algorithm, in its present form, does not claim to rigorously verify, its capabilities are nonetheless impressive.

References

[1] O. Aberth. Computation of topological degree using interval arithmetic, and applications. Math. Comp., 62(205):171–178, January 1994.

[2] J. Dian. Existence Verification of Higher Degree Singular Zeros of Nonlinear Systems. PhD thesis, University of Louisiana at Lafayette, 2000.

[3] E. R. Hansen. Global Optimization Using Interval Analysis. Marcel Dekker, Inc., New York, 1992.

[4] R. B. Kearfott. An efficient degree-computation method for a generalized method of bisection. Numer. Math., 32:109–127, 1979.

[5] R. B. Kearfott. Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht, Netherlands, 1996.

[6] R. B. Kearfott, J. Dian, and A. Neumaier. Existence verification for singular zeros of nonlinear systems. SIAM J. Numer. Anal., 38(2):360–379, 2000.

[7] A. Neumaier. Interval Methods for Systems of Equations. Cambridge University Press, Cambridge, England, 1990.

[8] H. Ratschek and J. Rokne. New Computer Methods for Global Optimization. Wiley, New York, 1988.

[9] A. J. Sommese, J. Verschelde, and C. W. Wampler. Numerical decomposition of the solution sets of polynomial systems into irreducible components, 2000. Accepted for publication in SIAM J. Numer. Anal.

[10] F. Stenger. An algorithm for the topological degree of a mapping in ℝ^{n}. Numer. Math., 25:23–38, 1976.