EXISTENCE VERIFICATION FOR SINGULAR ZEROS OF NONLINEAR SYSTEMS^{∗}

R. BAKER KEARFOTT^{†}, JIANWEI DIAN^{‡}, AND A. NEUMAIER^{§}

Abstract. Computational fixed point theorems can be used to automatically verify existence and uniqueness of a solution to a nonlinear system of equations F (x) = 0, F : ℝ^{n} → ℝ^{n} within a given region x of n-space. But such computations succeed only when the Jacobi matrix F′(x) is nonsingular everywhere in x. However, in many practical problems, the Jacobi matrix is singular, or nearly so, at the solution x^{∗}, F (x^{∗}) = 0. In such cases, arbitrarily small perturbations of the problem result in problems F̃(x) = 0 either with no solutions in x or with more than one solution in x; thus no general computational technique can prove existence and uniqueness. This leads to a fundamental philosophical problem: “What is meant by existence and uniqueness in such cases?”

Here, an interpretation of verification is given in the singular context: proof that a given number of true solutions exist within a region in complex space containing x.

Proof that a given number of true solutions exist within a given region of complex space is possible by computation of 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.

Key words. nonlinear systems, interval computations, verified computations, singularities, topological degree

AMS subject classifications. 65G10, 65H10

1. Introduction. Given an approximate solution x̌ to a nonlinear system of equations F (x) = 0, F : ℝ^{n} → ℝ^{n}, it is useful in various contexts to construct bounds around x̌ in which it is proven that there exists a unique solution x^{∗}, F (x^{∗}) = 0. For continuously differentiable 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 [5, 9, 10]. However, if F′(x^{∗}) is ill-conditioned or singular, such computations necessarily must fail. In this singular situation, what is meant by “unique solution” can be defined in several ways. A particular interpretation may be most appropriate in a particular practical context. Also, given a particular interpretation, various computational procedures can be developed to verify existence and uniqueness. In this paper, one interpretation and its corresponding computational procedure are considered.

1.1. Notation. We assume familiarity with the fundamentals of interval arithmetic; see [9, 10] for an introduction in the present context. (The works [2, 5, 11] 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, x = (x_{1}, …, x_{n}) denotes

∗This work was supported by National Science Foundation grant DMS-9701540.

†Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504, USA

(rbk@usl.edu).

‡Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504, USA

(dian@usl.edu).

§Institut für Mathematik, Universität Wien, Strudhofgasse 4, A-1050 Wien, Austria

(neum@cma.univie.ac.at).

an interval vector, A = (a_{i,j}) denotes a point matrix, and A = (a_{i,j}) denotes an interval matrix. Real n-space will be denoted by ℝ^{n}, while the set of n-dimensional interval vectors, i.e. n-dimensional boxes, will be denoted by 𝕀ℝ^{n}. Similarly, complex n-space will be denoted by ℂ^{n}, while the space of n-dimensional complex interval vectors will be denoted by 𝕀ℂ^{n}. The midpoint of an interval or interval vector x will be denoted by m(x). The non-oriented boundary of a box x will be denoted by ∂x while its oriented boundary will be denoted by b(x). (See §2.)

1.2. 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 operator G(x); G(x) ⊆ x then implies existence and uniqueness of F (x) = 0 within x. To describe these, we review

Definition 1.1. ( [10, p. 174], etc.) Let F : ℝ^{n} → ℝ^{m}. The matrix A is said to be a Lipschitz matrix for F over x provided, for every x ∈ x and y ∈ x, F (x) − F (y) = A(x − y) for some A ∈ A.

Most interval Newton methods for F : ℝ^{n} → ℝ^{n}, abstractly, are of the general form

x̃ = N(F; x, x̌) = x̌ + v, (1.1)

where v is computed to contain the solution set to the interval linear system

Av = −F (x̌), (1.2)

and where, for initial uniqueness verification, A is generally a Lipschitz matrix^{1} 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

Lemma 1.2. ([9, 10]) Suppose x̃ = x̌ + v 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 equation (1.2) to be considered here is the interval Gauss–Seidel method, defined by:

Definition 1.3. The preconditioned interval Gauss–Seidel image GS(F; x, x̌) of a box x is defined as GS(F; x, x̌) ≡ (x̃_{1}, …, x̃_{n}), where x̃_{i} is defined sequentially for i = 1 to n by:

x̃_{i} ≡ x_{i} ∩ (x̌_{i} − N_{i}/(Y_{i}A_{i})),

where

N_{i} = Y_{i}F (x̌) + ∑_{j=1}^{i−1} Y_{i}A_{j}(x̃_{j} − x̌_{j}) + ∑_{j=i+1}^{n} Y_{i}A_{j}(x_{j} − x̌_{j}),

and where x̌ = (x̌_{1}, …, x̌_{n})^{T} is an initial guess point, Y A ∈ 𝕀ℝ^{n×n} and Y F (x̌) are the matrix and right-hand-side vector for the preconditioned interval system Y A(x−x̌) = −Y F (x̌), Y ∈ ℝ^{n×n} is a point preconditioning matrix, Y_{i} denotes the i-th row of Y, and A_{j} denotes the j-th column of A.

1However, see [9, 12] for techniques for using slope matrices.

Lemma 1.2 applies when N(F; x, x̌) = GS(F; x, x̌), provided we specify that GS(F; x, x̌) be in the interior int(x) of x. (We must specify the interior because of the intersection step in Definition 1.3.) In particular, we have

Theorem 1.4. ([9, 10]) Suppose F : x ⊂ ℝ^{n} → ℝ^{n} and 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 unique x^{∗} ∈ x with F (x^{∗}) = 0. Various authors have proven Theorem 1.4; see [9, 10]. In particular, Miranda’s theorem can be used to easily prove Theorem 1.4 for N(F; x, x̌) = GS(F; x, x̌); see [9, p. 60].

Example 1. Take

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − x_{2},
f_{2}(x_{1}, x_{2}) = x_{1} − x_{2}^{2},

and x = (x_{1}, x_{2})^{T} = ([−0.1, 0.1], [−0.1, 0.3])^{T} . There is a unique root x^{∗} = (0, 0)^{T} of F = (f_{1}, f_{2})^{T} within x. In Example 1, if x̌ = (0, .1)^{T}, then F (x̌) = (−0.1, −0.01)^{T}, and an interval extension of the Jacobi matrix is

F′(x) = ( \begin{matrix}2x_{1} & −1 \\ 1 & −2x_{2} \end{matrix} ) = ( \begin{matrix}[−.2, .2] & −1 \\ 1 & [−0.6, 0.2] \end{matrix} ) .

If the preconditioning matrix Y is taken to be the inverse of the midpoint matrix of F′(x), then

Y = \{m(F′(X))\}^{−1} = ( \begin{matrix}−0.2 & 1 \\ −1 & 0 \end{matrix} ) .

We obtain, rounded out to four digits,

GS(F; x, x̌) = ([−0.07292, 0.09375], [−0.01875, 0.01875])^{T} ⊂ int(x).

Therefore, this computation proves that there is a unique solution to F (x) = 0 within x.

Inclusion in the interval Gauss–Seidel method is made possible because the inverse midpoint preconditioner reduces the interval Jacobi matrix to approximately a diagonal matrix.

1.3. Singularities: Philosophical Considerations. Theorem 1.4 is applicable only when the matrix A is regular, i.e. when A does not contain singular matrices. But consider

Example 2. Take

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − x_{2},
f_{2}(x_{1}, x_{2}) = x_{1}^{2} + x_{2},

and x = (x_{1}, x_{2})^{T} = ([−0.001, 0.001], [−0.001, 0.001])^{T} . Even though there is a unique root x^{∗} = (0, 0)^{T} of F = (f_{1}, f_{2})^{T} within x when F is as in Example 2, the interval Gauss–Seidel method cannot prove this, since the Jacobi matrix F′(x^{∗}) is singular. In fact, the interval Jacobi matrix is computed to be

F′(x) = ( \begin{matrix}2x_{1} & −1 \\ 2x_{1} & 1 \end{matrix} ) = ( \begin{matrix}[-0.002,0.002] & −1 \\ [-0.002,0.002] & 1 \end{matrix} ) ,

and the midpoint

m(F′(x)) = ( \begin{matrix}0 & −1 \\ 0 & 1 \end{matrix} ) ,

matrix the inverse of which is often used as the preconditioner matrix Y, is not invertible^{2}.

Symbolic methods can be used to show that Example 2 has a unique solution in a small region containing x_{1} = 0, x_{2} = 0. However, arbitrarily small perturbations of the problem result in either no solutions or two solutions. Consider

Example 3. Take

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − x_{2},
f_{2}(x_{1}, x_{2}) = x_{1}^{2} + x_{2} + ϵ,

and x = (x_{1}, x_{2})^{T} = ([−0.001, 0.001], [−0.001, 0.001])^{T} . Here, |ϵ| is very small.

The system in Example 3 has two solutions for ϵ < 0 and no solutions for ϵ > 0. But roundout in computer arithmetic and, perhaps, uncertainties in the system itself due to modelling or measurement uncertainties, make it impossible to distinguish systems such as in Example 3 for different ϵ, especially when computer arithmetic is used as part of the verification process. In such instances, a new interpretation needs to be given to existence / uniqueness. One expects the appropriate interpretation to be application-dependent in general, but the following is possible: verify that the system has an exact number of solutions within a larger space containing the original space. This paper will be based on this interpretation.

To be specific, we can extend an n-dimensional box in ℝ^{n} to an n-dimensional box in ℂ^{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 §2 for an explanation of multiplicity.) When the imaginary parts are small, we can view the complex roots as approximating a multiple real root. For example, the degree of the system in Example 3 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 §8.) The topological degree corresponds roughly to algebraic degree in one dimension; for example, the degree of z^{n} in a small region in ℂ^{1} containing 0 is n.

1.4. Organization of This Paper. A review of properties of the topological degree, to be used later, appears in §2. The issue of preconditioning appears in §3. Construction of the box in the complex space appears in §4.

Several algorithms have previously been proposed for computing the topological degree [1, 7, 14], but these require computational effort equivalent to finding all solutions to 4n (2n-1)-dimensional nonlinear systems within a given box, or worse. In §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.

2Alternate preconditioners can nonetheless be computed; see [9]. However, it can be shown that uniqueness cannot be proven in this case; see [9, 10].

In §6, the actual algorithm is presented and its computational complexity is given. Test problems and the test environment are described in §7. Numerical results appear in §8. Future directions appear in §9.

2. Review of Some Elements of Degree Theory. The topological degree or Brouwer degree, well-known within algebraic topology and nonlinear functional analysis, 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, 4, 7]. We only present the material we need.

Here, we explain what we mean by “multiplicity”. Actually, there is a more general concept index (see [4, 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. So, we use multiplicity as an alternative term for index since it’s more suggestive.

Suppose that F : D ⊂ ℂ^{n} → ℂ^{n} is analytic. Then the real and imaginary components of F and its argument z ∈ ℂ^{n} may be viewed as real components in ℝ^{2n}. Let z = x + iy and F (z) = u(x, y) + iv(x, y), where x = (x_{1}, …, x_{n}), y = (y_{1}, …, y_{n}), u(x, y) = (u_{1}(x, y), …, u_{n}(x, y)) and v(x, y) = (v_{1}(x, y), …, v_{n}(x, y)). We can define D̃ by

D̃ ≡ {(x_{1}, y_{1}, …, x_{n}, y_{n})|(x_{1} + iy_{1}, …, x_{n} + iy_{n}) ∈ D}

and F̃ : D̃ ⊂ ℝ^{2n} → ℝ^{2n} by

F̃ = (u_{1}, v_{1}, …, u_{n}, v_{n}).

Then, we have the following property of topological degree d(F̃, D̃, 0), and relationships between d(F̃, D̃, 0) and the solutions of the system F (z) = 0 in D.

Theorem 2.1. Suppose F : D ⊂ ℂ^{n} → ℂ^{n} is analytic, with F (z) ≠ 0 for any z ∈ ∂D, and suppose D̃ and F̃ : D̃ → ℝ^{2n} are defined as above. Then

1. d(F̃, D̃, 0) ≥ 0.

2. d(F̃, D̃, 0) > 0 if and only if there is a solution z^{∗} ∈ D, F (z^{∗}) = 0.

3. d(F̃, D̃, 0) is equal to the number of solutions z^{∗} ∈ D, F (z^{∗}) = 0, counting multiplicities.

4. If the Jacobi matrix F′(z^{∗}) is non-singular at every z^{∗} ∈ D with F (z^{∗}) = 0, then d(F̃, D̃, 0) is equal to the number of solutions z^{∗} ∈ D, F (z^{∗}) = 0. To obtain Theorem 2.1, we need to notice that F (z^{∗}) = 0 for z^{∗} ∈ D is equivalent

to F̃(x^{∗}, y^{∗}) = 0 for (x^{∗}, y^{∗}) ∈ D̃, where z^{∗} = x^{∗} + iy^{∗}. In Theorem 2.1, the first conclusion is easily obtained by knowing |d(F_{R}, D_{R}, 0)| ≤ d(F̃, D̃, 0), where F_{R} is function F (z) with z in a corresponding n-dimensional real domain D_{R} ∈ ℝ^{n}. (See [4, p. 49].) The third conclusion is actually Theorem (5.2) in Chapter I of [4], considering the definition of topological degree. The fourth conclusion can be obtained by the third conclusion, Theorem (7.2) and Lemma (9.3-2) in Chapter I of [4]. As for the second conclusion, it’s easy to see that d(F̃, D̃, 0) > 0 implies there is a solution (x^{∗}, y^{∗}) ∈ D̃, F̃(x^{∗}, y^{∗}) = 0. (See Existence Theorem (6.6) in Chapter I of [4].) The reverse is obtained by noticing the third conclusion and the fact that the multiplicities are positive.

The following three theorems will lead to the degree computation formula in Theorem 5.1 in §5, the formula used in our computational scheme.

Theorem 2.2. (See [13, §4.2].) Let D be an n-dimensional connected, oriented region in ℝ^{n} and F = (f_{1}, …, f_{n}), where f_{k}, k = 1, …, n, are continuous 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 β_{n−1}^{k}, k = 1, …, r, and there is p, 1 ≤ p ≤ n, such that:

1. F_{¬p} ≡ (f_{1}, …, f_{p−1}, f_{p+1}, …, f_{n}) ≠ 0 on the oriented boundary b(β_{n−1}^{k}) of β_{n−1}^{k}, k = 1, …, r; and

2. f_{p} has the same sign at all solutions of F_{¬p} = 0, if there are any, on β_{n−1}^{k}, 1 ≤ k ≤ r.

Choose s ∈ {−1, +1} and let K_{0}(s) denote the subset of the integers k ∈ {1, …, r} such that F_{¬p} = 0 has solutions on β_{n−1}^{k} and sgn(f_{p}) = s at each of those solutions. Then

d(F, D, 0) = (−1)^{p−1}s ∑_{k∈K_{0}(s)} d(F_{¬p}, β_{n−1}^{k}, 0).

The formula in Theorem 2.2 is a combination of formulas (4.15) and (4.16) in [13]. The orientation of D is positive and the orientations of β_{n−1}^{k}, whether positive or negative, are induced by the orientation of D. If we assume that the Jacobi matrices of F_{¬p} are non-singular at all solutions of F_{¬p} = 0 on β_{n−1}^{k}, then

d(F_{¬p}, β_{n−1}^{k}, 0) = t(β_{n−1}^{k}) ∑_{{x∈β{k \atop n−1} \atop F¬p=0}} sgn(JF_{¬p}(x)),

where t(β_{n−1}^{k}) = +1 or −1 depending on whether β_{n−1}^{k} 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 (7.2) in Chapter I of [4].) So, under this assumption, we can simplify the formula in Theorem 2.2.

Theorem 2.3. Suppose that all the conditions of Theorem 2.2 are satisfied and the Jacobi matrices of F_{¬p} are non-singular at all solutions of F_{¬p} = 0 on β_{n−1}^{k} for each k ∈ K_{0}(s). Then

d(F, D, 0) = (−1)^{p−1}s ∑_{k∈K_{0}(s)} t(β_{n−1}^{k}) ∑_{{x∈β{k \atop n−1} \atop F¬p(x)=0}} sgn(JF_{¬p}(x)),

where t(β_{n−1}^{k}) = +1 or −1 depending on whether β_{n−1}^{k} 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

x = (x_{1}, …, x_{n}),

where n ≥ 2 and x_{k} = [\underline{x}_{k}, \overline{x}_{k}]. The boundary ∂x of x consists of 2n (n − 1)- dimensional boxes

x\underline{k} ≡ (x_{1}, …, x_{k−1}, \underline{x}_{k}, x_{k+1}, …, x_{n})

and

x\overline{k} ≡ (x_{1}, …, x_{k−1}, \overline{x}_{k}, x_{k+1}, …, x_{n}),

where k = 1, …, n.

[figure: x\underline{1}; x\overline{1}; x\underline{2}; x\overline{2}]

Fig. 2.1. Positive center decomposition of x when n = 2.

Before introducing the next theorem, for completeness in the exposition, we construct the so-called positive or negative center decomposition of x (although such decompositions are undoubtedly known among experts). This center decomposition shows us how to specify orientations of the faces of x. We follow the convention that an n-simplex < P_{0}, P_{1}, …, P_{n} > in ℝ^{n} has positive or negative orientation if the determinant Δ_{n}((1, P_{0}), (1, P_{1}), …, (1, P_{n})) is positive or negative. Here, P_{i} = (x_{i1}, …, x_{in}), i = 0, 1, …, n, are n + 1 points in ℝ^{n} represented by row vectors, and (1, P_{i}) = (1, x_{i1}, …, x_{in}), i = 0, 1, …, n, are (n+1)-dimensional vectors that form the n + 1 rows of the determinant. See [3, 4, 7] for further details.

We construct the center decompositions of x inductively. Without loss of generality, assume x_{k} = [−1, +1], k = 1, …, n. So, the center of x is (0, …, 0). When n = 2 the positive center decomposition is shown in Figure 2.1, in which P_{0} = (0, 0). It’s easy to check that each of the 2^{2} · 2! 2-simplexes has positive orientation. So, x has positive orientation. Sometimes we say the orientation is +1 or −1 if it’s positive or negative. It’s clear that the orientations of x\underline{1}, x\overline{1}, x\underline{2} and x\overline{2} are −1, +1, +1 and −1. The negative center decomposition of x can be obtained by taking all the simplexes with the negative orientation. Sometimes we say a center decomposition with sign +1 or −1 if it’s positive or negative, respectively.

Suppose we already have center decompositions of (n − 1)-dimensional (n > 2) boxes with each decomposition consisting of 2^{n−1} · (n − 1)! (n − 1)-simplexes . We construct the positive center decomposition of an n-dimensional box x as follows. We take the center decomposition with sign (−1)^{k} of the (n − 1)-dimensional box x\underline{k} and we take the center decomposition with sign (−1)^{k+1} of the (n − 1)-dimensional box x\overline{k}, where k = 1, …, n.

We use

< Q_{0}^{(k,−1)}, Q_{1}^{(k,−1)}, …, Q_{n−1}^{(k,−1)} >

or

< Q_{0}^{(k,+1)}, Q_{1}^{(k,+1)}, …, Q_{n−1}^{(k,+1)} >

to denote an (n − 1)-simplex in the center decomposition of x\underline{k} or x\overline{k}, where

Q_{i} = (x_{i1}, …, x_{i(k−1)}, x_{i(k+1)}, …, x_{in})

and

Q_{i}^{(k,±1)} = (x_{i1}, …, x_{i(k−1)}, ±1, x_{i(k+1)}, …, x_{in}).

So,

sgn(Δ_{n−1}((1, Q_{0}), (1, Q_{1}), …, (1, Q_{n−1})))
= (−1)^{k} on x\underline{k} and (−1)^{k+1} on x\overline{k}.

For each (n − 1)-simplex

< Q_{0}^{(k,−1)}, Q_{1}^{(k,−1)}, …, Q_{n−1}^{(k,−1)}) >

or

< Q_{0}^{(k,+1)}, Q_{1}^{(k,+1)}, …, Q_{n−1}^{(k,+1)}) >,

we add the center P_{0} = (0, 0, …, 0) of x as the first point to form an n-simplex

< P_{0}, Q_{0}^{(k,−1)}, Q_{1}^{(k,−1)}, …, Q_{n−1}^{(k,−1)} >

or

< P_{0}, Q_{0}^{(k,+1)}, Q_{1}^{(k,+1)}, …, Q_{n−1}^{(k,+1)} > .

It’s easy to see that there are altogether 2^{n} · n! n-simplexes. Next we verify that each of those n-simplexes has the positive orientation.

Δ_{n}((1, P_{0}), (1, Q_{0}^{(k,−1)}), (1, Q_{1}^{(k,−1)}), …, (1, Q_{n−1}^{(k,−1)}))

= Δ_{n−1}(Q_{0}^{(k,−1)}, Q_{1}^{(k,−1)}, …, Q_{n−1}^{(k,−1)})
= (−1)^{k}Δ_{n−1}((1, Q_{0}), (1, Q_{1}), …, (1, Q_{n−1})).
Δ_{n}((1, P_{0}), (1, Q_{0}^{(k,+1)}), (1, Q_{1}^{(k,+1)}), …, (1, Q_{n−1}^{(k,+1)}))
= Δ_{n−1}(Q_{0}^{(k,+1)}, Q_{1}^{(k,+1)}, …, Q_{n−1}^{(k,+1)})
= (−1)^{k−1}Δ_{n−1}((1, Q_{0}), (1, Q_{1}), …, (1, Q_{n−1})).

So,

sgn(Δ_{n}((1, P_{0}), (1, Q_{0}^{(k,−1)}), (1, Q_{1}^{(k,−1)}), …, (1, Q_{n−1}^{(k,−1)})))
= (−1)^{k}(−1)^{k} = +1, and
sgn(Δ_{n}((1, P_{0}), (1, Q_{0}^{(k,+1)}), (1, Q_{1}^{(k,+1)}), …, (1, Q_{n−1}^{(k,+1)})))
= (−1)^{k−1}(−1)^{k+1} = +1.

This means we actually do get a positive center decomposition of x. We can get a negative center decomposition by taking all the 2^{n} · n! n-simplexes with negative orientation.

From the above we know that if we decompose x by the positive center decomposition, then x has positive orientation, the induced orientation of x\underline{k} is (−1)^{k}, and the induced orientation of x\overline{k} is (−1)^{k+1}. b(x) can be divided into x\underline{k} and x\overline{k}, k = 1, …, n, with the associated orientations. So, F ≠ 0 on b(x) is the same as F ≠ 0 on ∂x. Similarly, F_{¬p}(x) = 0 on b(x\underline{k}) or b(x\overline{k}) is the same as F_{¬p}(x) = 0 on ∂x\underline{k} or ∂x\overline{k}. Let \underline{K_{0}(s)} denote the subset of the integers k ∈ \{1, …, n\} such that F_{¬p} = 0 has solutions on x\underline{k} and sgn(f_{p}) = s at all the solutions, and \overline{K_{0}(s)} denote the subset of the integers k ∈ \{1, …, n\} such that F_{¬p} = 0 has solutions on x\overline{k} and sgn(f_{p}) = s at all the solutions, where s ∈ {−1, +1}. Then, by Theorem 2.3, we have

Theorem 2.4. Suppose F ≠ 0 on ∂x, and there is p, 1 ≤ p ≤ n, such that:

1. F_{¬p} ≠ 0 on ∂x\underline{k} or ∂x\overline{k}, k = 1, …, n;

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} ) .

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

2. f_{p} has the same sign at all solutions of F_{¬p} = 0, if there are any, on x\underline{k} or x\overline{k}, 1 ≤ k ≤ n; and

3. the Jacobi matrices of F_{¬p} are non-singular at all solutions of F_{¬p} = 0 on ∂x.

Then

d(F, x, 0)

= (−1)^{p−1}s\{ ∑_{k∈\underline{K_{0}(s)}}(−1)^{k} ∑_{{x∈x\underline{k} \atop F¬p(x)=0}} sgn |\frac{∂F_{¬p}}{∂x_{1}x_{2} … x_{k−1}x_{k+1} … x_{n}} (x)|

+ ∑_{k∈\overline{K_{0}(s)}}(−1)^{k+1} ∑_{{x∈x\overline{k} \atop F¬p(x)=0}} sgn |\frac{∂F_{¬p}}{∂x_{1}x_{2} … x_{k−1}x_{k+1} … x_{n}} (x)| \}.

3. On Preconditioning. The inverse midpoint preconditioner approximately diagonalizes the interval Jacobi matrix, when F′(x^{∗}) is non-singular (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′(x^{∗}) is n − p for some p > 0, then 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.

In the analysis to follow, we assume that the system has already been preconditioned, so that it is, to within second-order terms with respect to w(x), of the form in Figure 3.1. In the rest of this paper, we concentrate on the case p=1, although the idea can be applied to the general case.

4. The Complex Setting and System Form. In the remainder of this paper, we assume

1. F : D ⊂ ℝ^{n} → ℝ^{n} can be extended to an analytic function in ℂ^{n}.

2. x = (x_{1}, …, x_{n}) = ([\underline{x}_{1}, \overline{x}_{1}], …, [\underline{x}_{n}, \overline{x}_{n}]) is a small box that will be constructed centered at the approximate solution x̌, i.e. m(x) = (x̌_{1}, …, x̌_{n}).

3. x̌ is near a point x^{∗} with F (x^{∗}) = 0, such that ∥x̌ − x^{∗}∥ is much smaller than the width of the box x, and width of the box x is small enough for a quadratic model to be accurate over the box x.

4. F has been preconditioned as in Figure 3.1, and F′(x^{∗}) has null space of dimension 1.

So,

f_{k}(x) = (x_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}} (x̌)(x_{n} − x̌_{n}) + O (∥x − x̌∥^{2})
for 1 ≤ k ≤ n − 1.
f_{n}(x) = \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}} (x̌)(x_{k} − x̌_{k})(x_{l} − x̌_{l}) + O (∥x − x̌∥^{3})

For F : ℝ^{n} → ℝ^{n}, extend F to complex space: x + iy, with y in a small box

y = (y_{1}, …, y_{n}) = ([\underline{y}_{1}, \overline{y}_{1}], …, [\underline{y}_{n}, \overline{y}_{n}]),

where y is centered at (0, …, 0). Define z ≡ (x_{1}, y_{1}, …, x_{n}, y_{n}) = ([\underline{x}_{1}, \overline{x}_{1}],[\underline{y}_{1}, \overline{y}_{1}],…, [\underline{x}_{n}, \overline{x}_{n}],[\underline{y}_{n}, \overline{y}_{n}]), u_{k}(x, y) ≡ ℜ(f_{k}(x + iy)) and v_{k}(x, y) ≡ ℑ(f_{k}(x + iy)). With this, define

F̃(x, y) ≡ (u_{1}(x, y), v_{1}(x, y), …, u_{n}(x, y), v_{n}(x, y)) : ℝ^{2n} → ℝ^{2n}.

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, for 1 ≤ k ≤ (n − 1),

\begin{matrix}u_{k}(x, y) & = & (x_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}} (x̌)(x_{n} − x̌_{n}) \\ & & +O (∥(x − x̌, y)∥^{2}) , \\ v_{k}(x, y) & = & y_{k} + \frac{∂f_{k}}{∂x_{n}} (x̌)y_{n} + O (∥(x − x̌, y)∥^{2}) , \end{matrix} \} (4.1)

and

\begin{matrix}u_{n}(x, y) & = & \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}} (x̌)(x_{k} − x̌_{k})(x_{l} − x̌_{l}) \\ & & −\frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}} (x̌)y_{k}y_{l} \\ & & +O (∥(x − x̌, y)∥^{3}) , \\ v_{n}(x, y) & = & ∑_{k=1}^{n} ∑_{l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}} (x̌)(x_{k} − x̌_{k})y_{l} \\ & & +O (∥(x − x̌, y)∥^{3}) . \end{matrix} \} (4.2)

5. Simplification of a Degree Computation Procedure. To use Theorem 2.4 to compute the topological degree d(F̃, z, 0) directly in a verification algorithm would require a global search of the 4n (n−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 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 §3 and §4, the verification can be reduced to 4n − 4 interval evaluations and four 1-dimensional searches.

To describe the simplification, define

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}).

Similarly define y\underline{k} and y\overline{k}. Also 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}.

Observe that, for 1 ≤ k ≤ (n − 1), F̃_{¬u_{n}}(x, y) = 0 on x\underline{k}

⟹  u_{k}(x, y) ≈ (\underline{x}_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}} (x̌)(x_{n} − x̌_{n}) = 0
⟹  |\underline{x}_{k} − x̌_{k}| = |∂f_{k}/∂x_{n}(x̌)||x_{n} − x̌_{n}|
⟹w(x_{k}) ≤ |∂f_{k}/∂x_{n}(x̌)|w(x_{n}), or \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(x_{n}).

Similarly, F̃_{¬u_{n}}(x, y) = 0 on x\overline{k}

⟹w(x_{k}) ≤ |∂f_{k}/∂x_{n}(x̌)|w(x_{n}), or \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(x_{n}).

Also observe that, for 1 ≤ k ≤ (n − 1), F̃_{¬u_{n}}(x, y) = 0 on y\underline{k},

⟹  u_{k}(x, y) ≈ \underline{y}_{k} + \frac{∂f_{k}}{∂x_{n}} (x̌)y_{n} = 0
⟹  |\underline{y}_{k}| = |∂f_{k}/∂x_{n}(x̌)||y_{n}|
⟹w(y_{k}) ≤ |∂f_{k}/∂x_{n}(x̌)|w(y_{n}), or \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(y_{n}).

Similarly, F̃_{¬u_{n}}(x, y) = 0 on y\overline{k}

⟹w(y_{k}) ≤ |∂f_{k}/∂x_{n}(x̌)|w(y_{n}), or \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(y_{n}).

Thus, if x_{n} is chosen so that

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

then it is unlikely that u_{k}(x, y) = 0 on either x\underline{k} or x\overline{k}.

Similarly, if y_{n} is chosen so that

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

then it is unlikely that v_{k}(x, y) = 0 on either y\underline{k} or y\overline{k}.

Here, the coefficient \frac{1}{2} is to take into consideration the fact that u_{k}(x, y) ≈ (x_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}} (x̌)(x_{n} − x̌_{n}) and v_{k}(x, y) ≈ y_{k} + \frac{∂f_{k}}{∂x_{n}} (x̌)y_{n} are only approximate equalities.

When ∂f_{k}/∂x_{n}(x̌) happens to be 0, we can think that \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} and \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|}

are +∞.

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̃_{¬u_{n}}(x, y) ≠ 0 on each of these faces. Elimination of these 4n − 4 faces only needs 4n − 4 interval evaluations. Then, we only need to search the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} for solutions of F̃_{¬u_{n}}(x, y) = 0, regardless of how large n is. This reduces total computational cost dramatically, since searching a face is expensive. The following theorem forms the theoretical basis of our algorithm in §6.1.

Theorem 5.1. Suppose

1. u_{k} ≠ 0 on x\underline{k} and x\overline{k}, and v_{k} ≠ 0 on y\underline{k} and y\overline{k}, k = 1, …, n − 1;

2. F̃_{¬u_{n}} = 0 has a unique solution on x\underline{n} and x\overline{n} with y_{n} in the interior of y_{n}, and F̃_{¬u_{n}} = 0 has a unique solution on y\underline{n} and y\overline{n} with x_{n} in the interior of x_{n};

3. u_{n} ≠ 0 at the four solutions of F̃_{¬u_{n}} = 0 in condition 2; and

4. the Jacobi matrices of F̃_{¬u_{n}} are non-singular at the four solutions of F̃_{¬u_{n}} = 0

in condition 2. Then

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

+ ∑_{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)| .

Proof. Condition 1 implies F̃ ≠ 0 on x\underline{k}, x\overline{k}, y\underline{k} and y\overline{k}, k = 1, …, n − 1, and 2 and 3 imply F̃ ≠ 0 on x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}. So, F̃ ≠ 0 on ∂z.

Condition 1 implies F̃_{¬u_{n}} ≠ 0 on ∂x\underline{k}, ∂x\overline{k}, ∂y\underline{k} and ∂y\overline{k}, k = 1, …, n − 1. ∂x\underline{n} consists of 2(n − 1) (n − 2)-dimensional boxes each of which is either embedded in some x\underline{k}, x\overline{k}, y\underline{k} or y\overline{k}, 1 ≤ k ≤ n − 1, or is embedded in y\underline{n} or y\overline{n}. So, by 2 and 3, F̃_{¬u_{n}} ≠ 0 on ∂x\underline{n}. Similarly, F̃_{¬u_{n}} ≠ 0 on ∂x\overline{n}, ∂y\underline{n} and ∂y\overline{n}. So, condition 1 in Theorem 2.4 is satisfied.

Condition 2 in Theorem 2.4 is automatically satisfied since F_{¬p} = 0 either has no solutions or a unique solution on x\underline{k}, x\overline{k}, y\underline{k}, or y\overline{k}, 1 ≤ k ≤ n.

Then, by noticing 4, the conditions of Theorem 2.4 are satisfied. The formula is thus obtained with s = +1.

The conditions of Theorem 5.1 will be satisfied when the system is that as described in §3 and §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 and then add them to get the degree.

In Theorem 5.1 we choose s = +1. We can also choose s = −1. That doesn’t make any difference in our context if we ignore higher order terms in the values of u_{n} at the solutions of F̃_{¬u_{n}} = 0 on the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}. To be specific, the four values of u_{n} are

u_{n} = \frac{1}{2}Δ(\underline{x}_{n} − x̌_{n})^{2} + O (∥(x − x̌, y)∥^{3}) ,
u_{n} = \frac{1}{2}Δ(\overline{x}_{n} − x̌_{n})^{2} + O (∥(x − x̌, y)∥^{3}) ,
u_{n} = −\frac{1}{2}Δ\underline{y}_{n}^{2} + O (∥(x − x̌, y)∥^{3}) ,
u_{n} = −\frac{1}{2}Δ\overline{y}_{n}^{2} + O (∥(x − x̌, y)∥^{3}) ,

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 \underline{y}_{n} (or \overline{y}_{n}) will be the same (or roughly the same) as the values of u_{n} as a function of \underline{x}_{n} − x̌_{n} (or \overline{x}_{n} − x̌_{n}). So, if we ignore higher order terms, the cost of verifying u_{n} < 0 and searching for solutions of F̃_{¬u_{n}} = 0 with u_{n} > 0 is approximately the same as the cost of verifying u_{n} > 0 and searching for solutions of F̃_{¬u_{n}} = 0 with u_{n} < 0.

Next, we will give a theorem which will further reduce the search cost by telling us how we should search. Before introducing the theorem, define

α_{k} ≡ \frac{∂f_{k}}{∂x_{n}} (x̌),  1 ≤ k ≤ n − 1,
α_{n} ≡ −1,
β_{kl} ≡ \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}} (x̌)  1 ≤ k ≤ n, 1 ≤ l ≤ n,
(5.3)  Δ ≡ ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}

Note β_{kl} = β_{lk}

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 d(F̃, z, 0) = 2.

Proof. Under the assumptions,

(5.4)  u_{k} = (x_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}),  1 ≤ k ≤ n − 1,

(5.5)  v_{k} = y_{k} + α_{k}y_{n},  1 ≤ k ≤ n − 1,
(5.6)  u_{n} = \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}(x_{k} − x̌_{k})(x_{l} − x̌_{l}) − \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}y_{k}y_{l}

(5.7)  v_{n} = ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}(x_{k} − x̌_{k})y_{l}

Due to the construction of the box z,

u_{k} = (\underline{x}_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}) ≠ 0

on x\underline{k} and x\overline{k}, and

v_{k} = \underline{y}_{k} + α_{k}y_{n} ≠ 0

on y\underline{k} and y\overline{k}, where k = 1, …, n − 1.

Next, we locate the solutions of F̃_{¬u_{n}} = 0 on x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}.

1) On x\underline{n};

(5.8)  u_{k} = 0 ⟹ x̃_{k} = x̌_{k} − α_{k}(\underline{x}_{n} − x̌_{n}),  1 ≤ k ≤ n − 1,
(5.9)  v_{k} = 0 ⟹ ỹ_{k} = −α_{k}y_{n},  1 ≤ k ≤ n − 1,

Plugging equations (5.8) and (5.9) into equations (5.6) and (5.7), we get

(5.10)  u_{n} = \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}(\underline{x}_{n} − x̌_{n})^{2} − \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}y_{n}^{2}
= \frac{1}{2}Δ(\underline{x}_{n} − x̌_{n})^{2} − \frac{1}{2}Δy_{n}^{2}
(5.11)  v_{n} = ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}(\underline{x}_{n} − x̌_{n})y_{n}
= Δ(\underline{x}_{n} − x̌_{n})y_{n}.

Then,

(5.12)  v_{n} = 0 ⟹ ỹ_{n} = 0,

since Δ ≠ 0. Thus, by equation (5.9)

(5.13)  ỹ_{n} = 0 ⟹ ỹ_{k} = 0,  1 ≤ k ≤ n − 1,

So, F̃_{¬u_{n}} = 0 has a unique solution

(x̃, ỹ) = (x̃_{1}, 0, …, x̃_{n−1}, 0, \underline{x}_{n}, 0)

on x\underline{n}.

Plugging (5.12) into (5.10), we get the u_{n} value at this solution, which is

(5.14)  u_{n} = \frac{1}{2}Δ(\underline{x}_{n} − x̌_{n})^{2}.

Next, we compute the determinant of the Jacobi matrix of F̃_{¬u_{n}} at this solution. Define

γ_{k} ≡ ∑_{l=1}^{n} β_{lk}α_{l}.

Noting equations (5.4), (5.5) and (5.7), we have

|\frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}y_{n}} (x̃, ỹ)| (5.15)

= | \begin{matrix}\frac{∂u_{1}}{∂x_{1}} & \frac{∂u_{1}}{∂y_{1}} & … & \frac{∂u_{1}}{∂x_{n−1}} & \frac{∂u_{1}}{∂y_{n−1}} & \frac{∂u_{1}}{∂y_{n}} \\ \frac{∂v_{1}}{∂x_{1}} & \frac{∂v_{1}}{∂y_{1}} & … & \frac{∂v_{1}}{∂x_{n−1}} & \frac{∂v_{1}}{∂y_{n−1}} & \frac{∂v_{1}}{∂y_{n}} \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ \frac{∂u_{n−1}}{∂x_{1}} & \frac{∂u_{n−1}}{∂y_{1}} & … & \frac{∂u_{n−1}}{∂x_{n−1}} & \frac{∂u_{n−1}}{∂y_{n−1}} & \frac{∂u_{n−1}}{∂y_{n}} \\ \frac{∂v_{n−1}}{∂x_{1}} & \frac{∂v_{n−1}}{∂y_{1}} & … & \frac{∂v_{n−1}}{∂x_{n−1}} & \frac{∂v_{n−1}}{∂y_{n−1}} & \frac{∂v_{n−1}}{∂y_{n}} \\ \frac{∂v_{n}}{∂x_{1}} & \frac{∂v_{n}}{∂y_{1}} & … & \frac{∂v_{n}}{∂x_{n−1}} & \frac{∂v_{n}}{∂y_{n−1}} & \frac{∂v_{n}}{∂y_{n}} \end{matrix} |

= | \begin{matrix}1 & 0 & … & 0 & 0 & 0 \\ 0 & 1 & … & 0 & 0 & α_{1} \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ 0 & 0 & … & 1 & 0 & 0 \\ 0 & 0 & … & 0 & 1 & α_{n−1} \\ 0 & −γ_{1}(\underline{x}_{n} − x̌_{n}) & … & 0 & −γ_{n−1}(\underline{x}_{n} − x̌_{n}) & −γ_{n}(\underline{x}_{n} − x̌_{n}) \end{matrix} |

= −(\underline{x}_{n} − x̌_{n}) | \begin{matrix}1 & 0 & … & 0 & 0 & 0 \\ 0 & 1 & … & 0 & 0 & α_{1} \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ 0 & 0 & … & 1 & 0 & 0 \\ 0 & 0 & … & 0 & 1 & α_{n−1} \\ 0 & γ_{1} & … & 0 & γ_{n−1} & γ_{n} \end{matrix} |

= −(\underline{x}_{n} − x̌_{n})(− ∑_{k=1}^{n} α_{k}γ_{k})

= (\underline{x}_{n} − x̌_{n}) ∑_{k=1}^{n} α_{k} ∑_{l=1}^{n} β_{lk}α_{l}

= (\underline{x}_{n} − x̌_{n}) ∑_{k=1}^{n} ∑_{l=1}^{n} β_{lk}α_{k}α_{l}

= (\underline{x}_{n} − x̌_{n}) ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}
= (\underline{x}_{n} − x̌_{n})Δ.

2) On x\overline{n}; Similarly, F̃_{¬u_{n}} = 0 has a unique solution (x̃, ỹ) on x\overline{n}. The u_{n} value at this solution is

(5.16)  u_{n} = \frac{1}{2}Δ(\overline{x}_{n} − x̌_{n})^{2}.

The determinant of the Jacobi matrix of F̃_{¬u_{n}} at this solution is

(5.17)  |\frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}y_{n}} (x̃, ỹ)| = (\overline{x}_{n} − x̌_{n})Δ

3) On y\underline{n};

(5.18)  u_{k} = 0 ⟹ x̃_{k} = x̌_{k} − α_{k}(x_{n} − x̌_{n}),  1 ≤ k ≤ n − 1,
(5.19)  v_{k} = 0 ⟹ ỹ_{k} = −α_{k}\underline{y}_{n},  1 ≤ k ≤ n − 1,

Plugging equations (5.18) and (5.19) into equations (5.6) and (5.7), we get

(5.20)  u_{n} = \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}(x_{n} − x̌_{n})^{2} − \frac{1}{2} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}\underline{y}_{n}^{2}

= \frac{1}{2}Δ(x_{n} − x̌_{n})^{2} − \frac{1}{2}Δ\underline{y}_{n}^{2}
(5.21)  v_{n} = ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}(x_{n} − x̌_{n})\underline{y}_{n}
= Δ(x_{n} − x̌_{n})\underline{y}_{n}.

Then,

(5.22)  v_{n} = 0 ⟹ x̃_{n} = x̌_{n},

since Δ ≠ 0. Thus, by equation (5.18)

(5.23)  x̃_{n} = x̌_{n} ⟹ x̃_{k} = x̌_{k},  1 ≤ k ≤ n − 1,

So, F̃_{¬u_{n}} = 0 has a unique solution

(x̃, ỹ) = (x̌_{1}, ỹ_{1}, …, x̌_{n−1}, ỹ_{n−1}, x̌_{n}, \underline{y}_{n})

on y\underline{n}.

Plugging (5.22) into (5.20), we get the u_{n} value at this solution, which is

(5.24)  u_{n} = −\frac{1}{2}Δ\underline{y}_{n}^{2}.

Next, we compute the determinant of the Jacobi matrix of F̃_{¬u_{n}} at this solution.

Noting the equations (5.4), (5.5) and (5.7), we have

|\frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}x_{n}} (x̃, ỹ)| (5.25)

= | \begin{matrix}\frac{∂u_{1}}{∂x_{1}} & \frac{∂u_{1}}{∂y_{1}} & … & \frac{∂u_{1}}{∂x_{n−1}} & \frac{∂u_{1}}{∂y_{n−1}} & \frac{∂u_{1}}{∂x_{n}} \\ \frac{∂v_{1}}{∂x_{1}} & \frac{∂v_{1}}{∂y_{1}} & … & \frac{∂v_{1}}{∂x_{n−1}} & \frac{∂v_{1}}{∂y_{n−1}} & \frac{∂v_{1}}{∂x_{n}} \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ \frac{∂u_{n−1}}{∂x_{1}} & \frac{∂u_{n−1}}{∂y_{1}} & … & \frac{∂u_{n−1}}{∂x_{n−1}} & \frac{∂u_{n−1}}{∂y_{n−1}} & \frac{∂u_{n−1}}{∂x_{n}} \\ \frac{∂v_{n−1}}{∂x_{1}} & \frac{∂v_{n−1}}{∂y_{1}} & … & \frac{∂v_{n−1}}{∂x_{n−1}} & \frac{∂v_{n−1}}{∂y_{n−1}} & \frac{∂v_{n−1}}{∂x_{n}} \\ \frac{∂v_{n}}{∂x_{1}} & \frac{∂v_{n}}{∂y_{1}} & … & \frac{∂v_{n}}{∂x_{n−1}} & \frac{∂v_{n}}{∂y_{n−1}} & \frac{∂v_{n}}{∂x_{n}} \end{matrix} |

= | \begin{matrix}1 & 0 & … & 0 & 0 & α_{1} \\ 0 & 1 & … & 0 & 0 & 0 \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ 0 & 0 & … & 1 & 0 & α_{n−1} \\ 0 & 0 & … & 0 & 1 & 0 \\ −γ_{1}\underline{y}_{n} & 0 & … & −γ_{n−1}\underline{y}_{n} & 0 & −γ_{n}\underline{y}_{n} \end{matrix} |

= −\underline{y}_{n} | \begin{matrix}1 & 0 & … & 0 & 0 & α_{1} \\ 0 & 1 & … & 0 & 0 & 0 \\ ⋮ & ⋮ & & ⋮ & ⋮ & ⋮ \\ 0 & 0 & … & 1 & 0 & α_{n−1} \\ 0 & 0 & … & 0 & 1 & 0 \\ γ_{1} & 0 & … & γ_{n−1} & 0 & γ_{n} \end{matrix} |

= −\underline{y}_{n}(− ∑_{k=1}^{n} α_{k}γ_{k})

= \underline{y}_{n} ∑_{k=1}^{n} α_{k} ∑_{l=1}^{n} β_{lk}α_{l}

= \underline{y}_{n} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{lk}α_{k}α_{l}

= \underline{y}_{n} ∑_{k=1}^{n} ∑_{l=1}^{n} β_{kl}α_{k}α_{l}
= \underline{y}_{n}Δ.

4) On y\overline{n}; Similarly, F̃_{¬u_{n}} = 0 has a unique solution (x̃, ỹ) on y\overline{n}. The u_{n} value at this solution is

(5.26)  u_{n} = −\frac{1}{2}Δ\overline{y}_{n}^{2}.

The determinant of the Jacobi matrix of F̃_{¬u_{n}} at this solution is

(5.27)  |\frac{∂F̃_{¬u_{n}}}{∂x_{1}y_{1} … x_{n−1}y_{n−1}x_{n}} (x̃, ỹ)| = \overline{y}_{n}Δ

Finally, we can use the formula in Theorem 5.1 to compute the topological degree d(F̃, z, 0).

If Δ > 0, then we know from the equations (5.14), (5.16), (5.24) and (5.26) that u_{n} is only positive at the solutions of F̃_{¬u_{n}} = 0 on x\underline{n} and x\overline{n}. We also know the signs of the determinants of the Jacobi matrices at the two solutions from equations (5.15) and (5.17). So,

d(F̃, z, 0) = −(−1) + (+1) = 2

If Δ < 0, then we know from the equations (5.14), (5.16), (5.24) and (5.26) that u_{n} is only positive at the solutions of F̃_{¬u_{n}} = 0 on y\underline{n} and y\overline{n}. We also know the signs of the determinants of the Jacobi matrices at the two solutions from equations (5.25) and (5.27). So,

d(F̃, z, 0) = +(+1) − (−1) = 2

The proof of Theorem 5.2 tells us approximately where we can expect to find the solutions of F̃_{¬u_{n}} = 0 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, y) = 0 for x_{k} gives sharper bounds x̃_{k} with w(x̃_{k}) = O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1. Similarly, if y_{n} is known precisely, formally solving v_{k}(x, y) = 0 for y_{k} gives sharper bounds ỹ_{k} with w(ỹ_{k}) = O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1. So, when we

search x\underline{n} (or x\overline{n}) for solutions of F̃_{¬u_{n}} = 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_{n}^{0} of y_{n}, we can solve v_{k}(x, y) = 0 for y_{k} to get sharper bounds ỹ_{k} with w(ỹ_{k}) = O (max(∥(x − x̌, y)∥^{2}, ∥y_{n}^{0}∥)), 1 ≤ k ≤ n − 1. Thus we get a small subface of x\underline{n} (or x\overline{n}) over which we can either use an interval Newton method to verify the existence and uniqueness of the zero of F̃_{¬u_{n}} or use mean-value extensions to verify that F̃_{¬u_{n}} has no zeros, depending on whether y_{n}^{0} is in the middle of y_{n} or not. So, we end up with searching over a 1-dimensional interval y_{n}. This further reduces the search cost. We can search y\underline{n} or y\overline{n} in a similar way.

The analysis in this section leads to a practical algorithm in the next section.

6. The Algorithm and Its Computational Complexity. In this section, we will present the actual algorithm and give its 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\underline{k} and x\overline{k}, and v_{k} ≠ 0 on y\underline{k} and y\overline{k}, where 1 ≤ k ≤ n − 1. In the search phase, we verify the unique solution of F̃_{¬u_{n}} = 0 on x\underline{n} and x\overline{n} with y_{n} in the interior of y_{n}, and on y\underline{n} and y\overline{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̃_{¬u_{n}} at the four solutions of F̃_{¬u_{n}} = 0, compute the degree contributions of the 4 faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{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 1 ≤ k ≤ n − 1

(a) Do for x\underline{k} and x\overline{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\underline{k} and y\overline{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\underline{n} and x\overline{n}

(a) i. Use mean-value extensions for u_{k}(x, y) = 0 to solve for x_{k} to get

sharper bounds x̃_{k} with width O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1. ii. If x̃_{k} ∩ x_{k} = ∅, then return the degree contribution of that face as

0. 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_{n}^{0} of y_{n} which is centered at 0.

(d) i. Use mean-value extensions for v_{k}(x, y) = 0 to solve for y_{k} to get sharper bounds ỹ_{k} with width O (max(∥(x − x̌, y)∥^{2}, ∥y_{n}^{0}∥)), 1 ≤ k ≤ n − 1, and thus to get a

subface x0 ; \underline{n} (or x0 ; \overline{n}) of x\underline{n} (or x\overline{n}.) ii. If ỹ_{k} ∩ y_{k} = ∅, then stop and signal failure.

(e) i. Set up an interval Newton method for F̃_{¬u_{n}} to verify existence and

uniqueness of a zero in the subface x0 ; \underline{n} (or x0 ; \overline{n}.) ii. If the zero can not be verified, then stop and signal failure.

(f) Inflate y_{n}^{0} as much as possible provided the existence and uniqueness of the zero of F̃_{¬u_{n}} can be verified over the corresponding subface, and thus get a subinterval y_{n}^{1} of y_{n}.

(g) In this step, we verify F̃_{¬u_{n}} = 0 has no solutions when y_{n} ∈ y_{n} \ y_{n}^{1}. y_{n} \ y_{n}^{1} has two separate parts; we denote the lower part by y_{n}^{l} and the upper part by y_{n}^{u}. We only present the processing of the lower part. The upper part can be processed in a similar way.

i. Do

A. Use mean-value extensions for v_{k}(x, y) = 0 to solve for y_{k} to get sharper bounds for y_{k}, 1 ≤ k ≤ n − 1, and thus to get a subface of x\underline{n} (or x\overline{n}.)

B. Compute the mean-value extensions F̃_{¬u_{n}} over the subface ob-

tained in the last step. C. If 0 ∈ F̃_{¬u_{n}}, then bisect y_{n}^{l}, update the lower part as a new y_{n}^{l} and cycle. If 0 ∉ F̃_{¬u_{n}}, then exit the loop.

ii. Do A. If \underline{y}_{n}^{1} ≤ \overline{y}_{n}^{l}, then exit the loop. B. y_{n}^{l} ⟵ [\overline{y}_{n}^{l}, \overline{y}_{n}^{l} + w(y_{n}^{l})]. C. Use mean-value extensions for v_{k}(x, y) = 0 to solve for y_{k} to

get sharper bounds for y_{k}, 1 ≤ k ≤ n − 1, and thus to get a subface of x\underline{n} (or x\overline{n}.)

D. Compute the mean-value extensions F̃_{¬u_{n}} over the subface ob-

tained in the last step. E. If 0 ∉ F̃_{¬u_{n}}, then cycle.

If 0 ∈ F̃_{¬u_{n}}, then y_{n}^{l} ⟵ [\underline{y}_{n}^{l}, mid(y_{n}^{l})] and cycle.

(h) i. Compute the mean-value extension of u_{n} over x0 ; \underline{n} (or x0 ; \overline{n}.) ii. If u_{n} < 0, then return the degree contribution of that face as 0.

(i) i. Compute |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \underline{n})| (or |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \overline{n})|).

ii. If 0 ∈ |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \underline{n})| (or 0 ∈ |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \overline{n})|), then

stop and signal failure.

(j) Use the formula in Theorem 5.1 to compute the degree contribution of

that face.

2. Do for y\underline{n} and y\overline{n}

(a) Same as Step 1a except change x_{k} to y_{k}, x̃_{k} to ỹ_{k}, x_{k} to y_{k} and u_{k} to v_{k}.

(b) Same as Step 1b.

(c) Same as Step 1c except change y_{n}^{0} to x_{n}^{0}, y_{n} to x_{n} and 0 to x̌_{n}.

(d) Same as Step 1d except change y_{k} to x_{k}, ỹ_{k} to x̃_{k}, y_{k} to x_{k}, x0 ; \underline{n} to y0 ; \underline{n},

x0 ; \overline{n} to y0 ; \overline{n}, x\underline{n} to y\underline{n} and x\overline{n} to y\overline{n}.

(e) Same as Step 1e except change x0 ; \underline{n} to y0 ; \underline{n} and x0 ; \overline{n} to y0 ; \overline{n}.

(f) Same as Step 1f except change y_{n}^{0} to x_{n}^{0}, y_{n}^{1} to x_{n}^{1} and y_{n} to x_{n}.

(g) Same as Step 1g except change y_{n} \ y_{n}^{1} to x_{n} \ x_{n}^{1}.

(h) Same as Step 1h except change x0 ; \underline{n} to y0 ; \underline{n} and x0 ; \overline{n} to y0 ; \overline{n}.

(i) Same as Step 1i except change

|\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \underline{n})| to |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}x_{n}} (y0 ; \underline{n})| and

|\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}y_{n}} (x0 ; \overline{n})| to |\frac{∂F̃_{¬un}}{∂x_{1}y_{1}...x_{n−1}y_{n−1}x_{n}} (y0 ; \overline{n})|.

(j) Same as Step 1j.

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 (x) = 0. w(x_{k}) should be much larger than |x̌_{k} − x_{k}^{∗}|. But, at the same time, it should not be too large since the quadratic model needs to be accurate over the box.

2. In the search phase, in Step 1b (or 2b), 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 u_{n} < 0 on that face. It did happen for all the test problems. (See §8 for the test results.)

3. In the search phase, in Step 1e ( or 2e), we precondition the system F̃_{¬u_{n}} before we use an interval Newton method, so that the method will succeed (see §1.2 and §3). The system F̃_{¬u_{n}} is nonsingular over the subfaces under consideration. In the actual implementation, LINPACK routines DGECO and DGESL were used to compute the preconditioners.

4. In the search phase, in Step 1f (or 2f), we first expand the subinterval y_{n}^{0} (or x_{n}^{0}) by ϵ = \frac{1}{2}w(y_{n}^{0}) at both ends. If the existence and uniqueness of the zero of F̃_{¬u_{n}} can be verified over the corresponding subface, then we expand the subinterval by 2ϵ at both ends, and then 4ϵ and so on until the existence and uniqueness of the zero can not be established.

5. In the search phase, in Step 1g (or 2g), the underlying idea is that the farther away the interval y_{n}^{l} is from the interval y_{n}^{0} whose corresponding subface of x\underline{n} (or x\overline{n}) contains a unique solution of F̃_{¬u_{n}} = 0 and/or the narrower the interval y_{n}^{l} is, the more probable it is that we can verify that F̃_{¬u_{n}} ≠ 0 over the subface of x\underline{n} (or x\overline{n}) which corresponds to y_{n}^{l}.

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 1ai and 1bi are of order O (n^{2}). Step 1aii and 1bii are of order O (1). Thus, the order of this phase is O (n^{3}).

Search Phase: Step 1a and 2a are of order O (n^{3}). Step 1b and 2b are of order O (n^{2}). Step 1c and 2c are of order O (1). Step 1d and 2d are of order O (n^{3}). Step 1e and 2e are of order O (n^{3}). Step 1f and 2f are of order N_{infl}*O (n^{3}). (See explanation below.) Step 1g and 2g are of order N_{proc}*O (n^{3}). (See explanation below.) Step 1h and 2h are of order O (n^{2}). Step 1i and 2i are of order O (n^{3}). Step 1j and 2j 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}).

The order of the overall algorithm is thus O (n^{3}).

\underline{Remarks}

1. There are two performance measures in the algorithm, N_{infl} and N_{proc}. N_{infl} is the number of inflations the algorithm did in Step 1f or 2f. N_{proc} is the number of subintervals of y_{n} \ y_{n}^{1} the algorithm processed in Step 1g or subintervals of x_{n} \ x_{n}^{1} the algorithm processed in Step 2g, i.e. number of y_{n}^{l}’s plus number of y_{n}^{u}’s in Step 1g or number of x_{n}^{l}’s plus number of x_{n}^{u}’s in Step 2g . (See the algorithm in §6.1).

2. The order of the algorithm cannot be improved, since computing preconditioners of the original system and the system F̃_{¬u_{n}} 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 [6, Lemma 2.4], we considered systems of the following form.

Example 4. Set

f(x) = h(x, t) = (1 − t)(Ax − x^{2}) − tx,

where A ∈ ℝ^{n×n} is the matrix corresponding to central difference discretization of the boundary value problem −u′′ = 0, u(0) = u(1) = 0 and x^{2} = (x_{1}^{2}, …, x_{n}^{2})^{T}. t was chosen to be equal to

t_{1} = λ_{1}/(1 + λ_{1}),

where λ_{1} is the largest eigenvalue of A.

The homotopy h in Example 4 has a simple bifurcation point at t = t_{1}, where the two paths cross obliquely. That is, there are two solutions to f(x) = 0 near x = 0, for 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 t = t_{1}.

The test set consists of Example 2, Example 3 with ϵ = +10^{−6} and −10^{−6}, and Example 4 with n = 5, 10, 20, 40, 80, 160, 320. For all the test problems, we used (0, 0, …, 0) as a good approximate solution to the problem F (x) = 0. Actually, it’s the exact solution in Example 2 and Example 4. w(x_{k}) and w(y_{k}) were set to 10^{−3} for 1 ≤ k ≤ n − 1. w(x_{n}) and w(y_{n}) were computed automatically by the algorithm. In fact, w(x_{k}) and w(y_{k}), 1 ≤ k ≤ n−1, 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 2 and Example 4 and the known approximate solution to Example 3 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 §6.1 was programmed in the Fortran 90 environment developed and described in [8, 9]. 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 1e, 2e, 1f and 2f of the search phase to compute the preconditioners. (See the algorithm and its explanation in §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 §7.1. n: number of independent variables. Success: whether the algorithm was successful. Degree: yopological degree returned by the algorithm. CPU Time: CPU time in seconds of the algorithm.

Table 8.1 Numerical Results

Problem  n Success Degree CPU Time Time Ratio
Example 2 2 Yes 2 0.0761

Example 3

(ϵ = +10^{−6}) 2 Yes 2 0.0511

Example 3

(ϵ = −10^{−6})  2  Yes  2  0.0513
Example 4  5  Yes  2  0.6806
Example 4 10  Yes  2  3.3403  4.91
Example 4 20  Yes  2  19.440  5.82
Example 4 40  Yes  2  140.34  7.22
Example 4 80  Yes  2  1123.6  8.01
Example 4 160  Yes  2  8891.3  7.91
Example 4 320  Yes  2  65395.5  7.36

Time Ratio: This only applies to Example 4. It’s the ratio of two successive CPU

times. The column labels of Table 8.2 are as follows.

Problem: names of the problems identified in §7.1. n: number of independent variables. N_{infl}: number of inflations the algorithm did in Step 1f or 2f for the indicated face

x\underline{n}, x\overline{n}, y\underline{n} or x\overline{y} (see explanation of the computational complexity of the algorithm in §6.2).

N_{proc}: number of intervals the algorithm processed in Step 1g or 2g for the indicated face x\underline{n}, x\overline{n}, y\underline{n} or x\overline{y} (see explanation of the computational complexity of the algorithm in §6.2).

We can see from Table 8.1 that the algorithm was successful on all the problems in the test set. The algorithm is of order 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. So, when n was small, those lower order structures had significant influence on the CPU time. But, when n was large, the O (n^{3}) terms dominated. We can see this from the time ratios of Example 4 in Table 8.1.

In Table 8.2, for all the problems, there were two faces of x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} for which N_{infl} = 0. This is because the algorithm verified that u_{n} < 0 on each of those two faces in Step 1b or 2b, and returned a degree contribution of each of those two faces as 0. So, the algorithm didn’t proceed to Step 1f or 2f. For the same reason, N_{proc} = 0 for those two faces. For the remaining two faces for which the algorithm did proceed to Step 1f or 2f, N_{infl} is small.

In Step 1g or 2g which immediately follows the inflations, N_{proc} = 0 for Example 2 and Example 3. This is because the inflations had covered the whole interval y_{n}. (See the algorithm in §6.1 and the explanation of its computational complexity in §6.2.) More significant is that N_{proc} = 2 in Example 4 regardless of small n or large n. This is because only one interval was processed to verify that F̃_{¬u_{n}} = 0 has no solutions when x_{n} ∈ x_{n}^{l} and only one interval was processed to verify that F̃_{¬u_{n}} = 0 has no solutions when x_{n} ∈ x_{n}^{u}. (See the algorithm in §6.1 and the explanation of its computational complexity in §6.2.) This means that the algorithm was quite efficient.

Table 8.2 Statistical Data

N_{infl}  N_{proc}
Problem  n x\underline{n} x\overline{n} y\underline{n} y\overline{n} x\underline{n} x\overline{n} y\underline{n} y\overline{n}
Example 2  2 6 6 0 0 0 0 0 0

Example 3

(ϵ = +10^{−6})  2 2 2 0 0 0 0 0 0

Example 3

(ϵ = −10^{−6})  2 2 2 0 0 0 0 0 0
Example 4  5 0 0 5 5 0 0 2 2
Example 4 10 0 0 5 5 0 0 2 2
Example 4 20 0 0 4 4 0 0 2 2
Example 4 40 0 0 4 4 0 0 2 2
Example 4 80 0 0 4 4 0 0 2 2
Example 4 160 0 0 4 4 0 0 2 2
Example 4 320 0 0 3 3 0 0 2 2

9. Conclusions and Future Work. When we tested the algorithm, we took advantage of knowing the true solutions. (See §7.1.) For this reason, we set w(x_{k}) and w(y_{k}), 1 ≤ k ≤ n − 1 somewhat 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 §5 and the algorithm in §6.1 can be generalized to general rank deficiency. Also, at present, it is assumed that the second derivatives \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}}, 1 ≤ k ≤ n, 1 ≤ l ≤ n don’t vanish simultaneously at the singular solution. In fact, the analysis in §5 and the algorithm in §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. But computing higher order derivatives 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.

Another future direction of this study is to apply the algorithms to bifurcation problems and other physical models.

REFERENCES

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

[2] G. Alefeld and J. Herzberger, Introduction to Interval Computations, Academic Press, New York, 1983.

[3] P. Alexandroff and H. Hopf, Topologie, Chelsea, 1935.

[4] J. Cronin, Fixed Points and Topological Degree in Nonlinear Analysis, American Mathematical Society, Providence, RI, 1964.

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

[6] H. Jürgens, H.-O. Peitgen, and D. Saupe, Topological perturbations in the numerical nonlinear eigenvalue and bifurcation problems, in Analysis and Computation of Fixed Points, S. M. Robinson, ed., New York, 1980, Academic Press, pp. 139–181.

[7] R. B. Kearfott, Computing the Degree of Maps and a Generalized Method of Bisection, PhD thesis, University of Utah, Department of Mathematics, 1977.

[8] , A Fortran 90 environment for research and prototyping of enclosure algorithms for nonlinear equations and global optimization, ACM Trans. Math. Software, 21 (1995), pp. 63–78.

[9] , Rigorous Global Search: Continuous Problems, Kluwer, Dordrecht, Netherlands, 1996.

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

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

[12] S. M. Rump, Verification methods for dense and sparse systems of equations, in Topics in Validated Computations, J. Herzberger, ed., Amsterdam, 1994, Elsevier Science Publishers, pp. 63–135.

[13] F. Stenger, Computing the topological degree of a mapping in ℝ^{n}, Numer. Math., 25 (1975), pp. 23–38.

[14] M. Stynes, An Algorithm for the Numerical Calculation of the Degree of a Mapping, PhD thesis, Oregon State University, Department of Mathematics, Corvallis, Oregon, 1977.