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

R. BAKER KEARFOTT^{†} AND JIANWEI DIAN^{‡}

Abstract. It is known that, in general, no computational techniques can verify the existence of a singular solution of the nonlinear system of n equations in n variables within a given region x of n-space. However, computational verification that a given number of true solutions exist within a region in complex space containing x is possible. That can be done by computation of the topological degree. In a previous paper, we presented theory and algorithms for the simplest case, when the rank-defect of the Jacobi matrix at the solution is one and the topological index is 2. Here, we will generalize that result to arbitrary topological index d ≥ 2: We present theory, algorithms, and experimental results. We also present a heuristic for determining the degree, obtaining a value that we can subsequently verify with our algorithms.

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

AMS subject classifications. 65G10, 65H10

1. Introduction. Our fundamental problem is

Given F : x → ℝ^{n} and x ∈ 𝕀ℝ^{n}, rigorously verify: • there exists a unique x^{∗} ∈ x such that F (x^{∗}) = 0, where • x = \{(x_{1}, x_{2}, …, x_{n})^{T} ∈ ℝ^{n} | \underline{x}_{i} ≤ x_{i} ≤ \overline{x}_{i}, 1 ≤ i ≤ n\}, (1.1)

where the \underline{x}_{i} and \overline{x}_{i} represent constant bounds on the problem variables x_{i}. When the Jacobi matrix F′(x^{∗}) well-conditioned and not too quickly varying, interval computations have no trouble proving that there is a unique solution within small boxes with x^{∗} reasonably near the center; see [4, 9, 11]. When F′(x^{∗}) is ill-conditioned or singular, in general, no computational techniques can verify the existence of a solution within a given region x of ℝ^{n}. However, in the singular case, computational but rigorous verification that a given number of true solutions exist within a region in complex space containing x is possible, as we indicated in [10]. In [10], we studied the simplest case, when the rank-defect of the Jacobi matrix at the solution is one, and we developed and experimentally validated algorithms for the case when the topological index is 2. There, we proved the special case of Theorem 3.1 when d = 2 under the same assumptions as those in §2, we developed specialized versions of the algorithms in §4, and we presented varying-dimensional experimental results.

We were surprised and pleased that the results in [10] could be generalized so easily. In particular, we developed an alternate simple, general proof for Theorem 3.1. Furthermore, the algorithms in §4, although not taking advantage of special efficiencies in the degree-2 case, are similar in structure and have the same computational complexity as the algorithms in [10].

1.1. Notation. We assume familiarity with the fundamentals of interval arithmetic; see [1, 4, 9, 11, 13] for introductory material.

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

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

(rbk@louisiana.edu).

‡Software Design Engineer, Hewlett–Packard Company, 3000 Waterview Parkway, Richardson,

TX 75080, USA (jianwei@rsn.hp.com).

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 an interval vector, A = (a_{i,j}) denotes a point matrix, and A = (a_{i,j}) denotes an interval matrix. The midpoint of an interval or interval vector x will be denoted by m(x). Real n-space will be denoted by ℝ^{n}, while complex n-space will be denoted by ℂ^{n}.

Suppose x = (x_{1}, …, x_{n}) is an n-dimensional real box, where x_{k} = [\underline{x}_{k}, \overline{x}_{k}]. The non-oriented boundary of x, denoted by ∂x, consists of 2n (n − 1)-dimensional real 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. If x is positively oriented, then the derived orientation of x\underline{k} is (−1)^{k} and the derived orientation of x\overline{k} is (−1)^{k+1}; see [10].

1.2. Formulas from Degree Theory. In [10], we reviewed the topological degree in the context of this paper. Also see [2, 3, 6, 7, 12, 14]. Here, we repeat several properties used in the proofs in subsequent sections.

Theorem 1.1. ([12, p. 150]) Suppose that the Jacobian matrix F′(x) is nonsingular at each zero of F. Then, the degree d(F, D, 0) is equal to the number of zeros of F at which the determinant of the Jacobian matrix F′(x) is positive minus the number of zeros of F at which the determinant of the Jacobian matrix F′(x) is negative.

Theorem 1.1 gives some intuition of what the degree is, under the conditions in the theorem. It’s a kind of counting of zeros of F in D.

Theorem 1.2. ([12, p. 150]) Let F, G : \overline{D} ⊂ ℝ^{n} → ℝ^{n} be two continuous functions. If F (x) = G(x) for x ∈ ∂D, then d(F, D, 0) = d(G, D, 0).

Theorem 1.2 states one of the most important properties of degree: the degree depends only on the function values on the boundary.

Theorem 1.3. ([12, p. 152]) Let α = min{∥F (x)∥_{2}|x ∈ ∂D}. If

sup\{∥F (x) − G(x)∥_{2}|x ∈ \overline{D}\} < \frac{1}{7}α,

then

d(F, D, 0) = d(G, D, 0).

Theorem 1.3 tells us that small perturbations of F don’t change the degree. Theorem 1.4. ([12, p. 157]) Let F, G : \overline{D} ⊂ ℝ^{n} → ℝ^{n} be two continuous

functions. If

0 ∉ {tF (x) + (1 − t)G(x)|x ∈ ∂D and t ∈ [0, 1]},

then

d(F, D, 0) = d(G, D, 0).

Theorem 1.4 is the famous Poincaré-Bohl Theorem. It’s a particular case of the homotopy invariant property of the topological degree.

Suppose F : \overline{D} ⊂ ℂ^{n} → ℂ^{n} is analytic, and view the real and imaginary components of F and its argument z ∈ ℂ^{n} 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 define D̃ by

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

and F̃ : \overline{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 1.5. ([10]) Suppose F : \overline{D} ⊂ ℂ^{n} → ℂ^{n} is analytic, with F (z) ≠ 0 for any z ∈ ∂D, and suppose D̃ and F̃ : \overline{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.

1.3. A Basic Degree Computation Formula. If we let

F_{¬k}(x) ≡ (f_{1}(x), …, f_{k−1}(x), f_{k+1}(x), …, f_{n}(x))

and select s ∈ {−1, 1}, then d(F, x, 0) is equal to the number of zeros of F_{¬k} on ∂x with positive orientation at which sgn(f_{k}) = s, minus the number of zeros of F_{¬k} on ∂x with negative orientation at which sgn(f_{k}) = s. The orientation of each zero can be computed by computing the sign of the determinant of the Jacobian of F_{¬k} and by taking into account the orientation of the face of x on which the zero lies.

Next, we present a degree computation formula that is similar to a formula used in [10]; see Theorem 2.5 of [10]. We can get the formula in the following theorem by noticing formulas (4.12) and (4.14) in [14] and by taking the orientations of the faces of x into account. We will use this formula to derive the computational procedures in §4.

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

1. F_{¬p} ≡ (f_{1}, …, f_{p−1}, f_{p+1}, …, f_{n}) ≠ 0 on ∂x\underline{k} or ∂x\overline{k}, k = 1, …, n; and 2. 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=1}^{n}(−1)^{k} ∑_{x∈x\underline{k} ; {F¬p(x)=0 \atop sgn(fp(x))=s}} sgn |\frac{∂F_{¬p}}{∂x_{1}x_{2} … x_{k−1}x_{k+1} … x_{n}} (x)|

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

where s = +1 or −1.

2. Assumptions and Choice of Box. In this section, we present the basic assumptions. We also introduce how we choose the coordinate bounds x_{i} = [\underline{x}_{i}, \overline{x}_{i}] to satisfy the assumptions and enable more efficient algorithms. When the rank of F′(x^{∗}) is n − p for some p > 0, an appropriate preconditioner can be used to reduce F′(x) to approximately the pattern shown in Figure 2.1. (See [9] and [10] for details on preconditioning.)

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. 2.1. A preconditioned singular system of rank n − p, where “∗” represents a non-zero element.

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 2.1. Here as in [10], we concentrate on the case p=1.

2.1. The Basic Assumptions. As in the special case d = 2 of [10], we assume 1. F : \overline{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 constructed to be centered at an 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 so that mean value interval extensions lead, after preconditioning, to a system like Figure 2.1, with small intervals replacing the zeros.

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

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.

The following representation of f(x) near x̌ is appropriate under these assumptions.

(2.1) f_{k}(x) = (x_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}) + O (∥x − x̌∥^{2})
for 1 ≤ k ≤ n − 1.
(2.2) f_{n}(x) = \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)(x_{k_{1}} − x̌_{k_{1}})(x_{k_{2}} − x̌_{k_{2}}) + …

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

+O (∥x − x̌∥^{d+1}) ,

or

(2.3) f_{k}(x) ≈ (x_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}) for 1 ≤ k ≤ n − 1.
(2.4) f_{n}(x) ≈ \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)(x_{k_{1}} − x̌_{k_{1}})(x_{k_{2}} − x̌_{k_{2}}) + …

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

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 x̃ ≡ (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}.

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

Then, based on (2.1) and (2.2), for 1 ≤ k ≤ (n − 1),

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

or

\begin{matrix}u_{k}(x, y) & ≈ & (x_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}), \\ v_{k}(x, y) & ≈ & y_{k} + α_{k}y_{n}. \end{matrix} \} (2.6)

2.2. Choosing the Coordinate Bounds. We use a similar scheme to that of §5 of [10]. In particular, having defined x\underline{k} and x\overline{k} in §1.1, we define y\underline{k} and y\overline{k} similarly:

y\underline{k} ≡ (x_{1}, y_{1}, …, x_{k−1}, y_{k−1}, x_{k}, \underline{y}_{k}, x_{k+1}, y_{k+1}, …, x_{n}, y_{n}) and
y\overline{k} ≡ (x_{1}, y_{1}, …, x_{k−1}, y_{k−1}, x_{k}, \overline{y}_{k}, x_{k+1}, y_{k+1}, …, x_{n}, y_{n}).

To compute the degree d(F̃, x̃, 0), we consider F̃_{¬u_{n}} on the boundary of x̃. This boundary 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}. We set x_{n} and y_{n} in such a way that

(2.7) w(x_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{\frac{w(x_{k})}{|α_{k}|}\} and w(y_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{\frac{w(y_{k})}{|α_{k}|}\} .

Constructing the box widths this way makes it is unlikely that u_{k}(x, y) = 0 on either x\underline{k} or x\overline{k} and unlikely that v_{k}(x, y) = 0 on either y\underline{k} or y\overline{k}, where k = 1, …, n − 1.

This, in turn, allows us to replace searches on 4n − 4 of the 4n faces of ∂x̃ by simple interval evaluations, reducing the total computational cost dramatically. See [10] for details.

A difference between the scheme used here and that of [10] is the way the ratio w(y_{n})/w(x_{n}) is chosen. In [10], w(y_{n}) was chosen large relative to x_{n}, to arrange no solutions of u_{n} = 0 on y\underline{n} and y\overline{n}. When the degree is odd, that is not possible, and we have found the strategy represented by formula (4.3) below, implying w(y_{n}) small relative to w(x_{n}) as in Figure 4.1 below, to be more convenient.

3. A Theorem. In [10] we proved that, under the assumptions in §2, if (2.3) and (2.4) are exact for d = 2, and if Δ_{1} = 0 but Δ_{2} ≠ 0, then d(F̃, x̃, 0) = 2. Here, we generalize that result to Δ_{1} = … = Δ_{d−1} = 0, Δ_{d} ≠ 0.

Theorem 3.1. Suppose

1. all the assumptions in §2 are true;

2. (2.3) and (2.4) are exact; and

3. Δ_{1} = … = Δ_{d−1} = 0, Δ_{d} ≠ 0, where 2 ≤ d.

Then d(F̃, x̃, 0) = d.

In contrast to the proof in [10], we use a homotopy argument to prove Theorem 3.1.

Proof. Let z = (z_{1}, …, z_{n}) = (x_{1} + iy_{1}, …, x_{n} + iy_{n}). Then

F (z) = (f_{1}(z), …, f_{n−1}(z), f_{n}(z)),

where

f_{k}(z) = (z_{k} − ž_{k}) + \frac{∂f_{k}}{∂x_{n}} (x̌)(z_{n} − ž_{n})

for 1 ≤ k ≤ n − 1,

(3.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}} − ž_{k_{1}})(z_{k_{2}} − ž_{k_{2}}) + …

+\frac{1}{d!} ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x̌)(z_{k_{1}} − ž_{k_{1}}) … (z_{k_{d}} − ž_{k_{d}}).

We construct G : ℂ^{n} → ℂ^{n} by

G(z) = (g_{1}(z), …, g_{n−1}(z), g_{n}(z)),

where

g_{k}(z) = (z_{k} − ž_{k}) for 1 ≤ k ≤ n − 1,
(3.2)  g_{n}(z) = \frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}.

Let p_{k}(x, y) ≡ ℜ(g_{k}(x + iy)) and q_{k}(x, y) ≡ ℑ(g_{k}(x + iy)). With this, define G̃ : ℝ^{2n} → ℝ^{2n} by

G̃(x, y) ≡ (p_{1}(x, y), q_{1}(x, y), …, p_{n}(x, y), q_{n}(x, y)).

We will first prove d(F̃, x̃, 0) = d(G̃, x̃, 0). Define

H̃((x, y), t) ≡ tF̃(x, y) + (1 − t)G̃(x, y)
and  H(z, t) ≡ tF (z) + (1 − t)G(z).

We will prove that H̃((x, y), t) ≠ 0 when (x, y) ∈ ∂x̃ and t ∈ [0, 1]. It’s clear that H̃((x, y), t) = 0 is equivalent to H(z, t) = 0, so we consider H(z, t).

H(z, t) = (H_{1}(z, t), …, H_{n}(z, t))
= (tf_{1}(z) + (1 − t)g_{1}(z), …, tf_{n−1}(z) + (1 − t)g_{n−1}(z),

tf_{n}(z) + (1 − t)g_{n}(z))

= (t((z_{1} − ž_{1}) + α_{1}(z_{n} − ž_{n})) + (1 − t)(z_{1} − ž_{1}), …,

t((z_{n−1} − ž_{n−1}) + α_{n−1}(z_{n} − ž_{n})) + (1 − t)(z_{n−1} − ž_{n−1}),

tf_{n}(z) + (1 − t)\frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d})

= ((z_{1} − ž_{1}) + tα_{1}(z_{n} − ž_{n}), …, (z_{n−1} − ž_{n−1}) + tα_{n−1}(z_{n} − ž_{n}),

tf_{n}(z) + (1 − t)\frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d})

Thus, H(z, t) = 0 implies z_{k} = ž_{k} − tα_{k}(z_{n} − ž_{n}) for k = 1, …, n − 1. Plugging z_{k} − ž_{k} = −tα_{k}(z_{n} − ž_{n}) for each such k into (3.1), we get

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̌)(−1)^{2}t^{2}α_{k_{1}}α_{k_{2}}(z_{n} − ž_{n})^{2} + … +

\frac{1}{d!} ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x̌)(−1)^{d}t^{d}α_{k_{1}} … α_{k_{d}}(z_{n} − ž_{n})^{d}

= \frac{(−1)^{2}t^{2}}{2!} (z_{n} − ž_{n})^{2} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)α_{k_{1}}α_{k_{2}} + …

+\frac{(−1)^{d}t^{d}}{d!} (z_{n} − ž_{n})^{d} ∑_{k_{1}=1}^{n} … ∑_{k_{d}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x̌)α_{k_{1}} … α_{k_{d}}

= \frac{(−1)^{2}t^{2}Δ_{2}}{2!} (z_{n} − ž_{n})^{2} + … + \frac{(−1)^{d−1}t^{d−1}Δ_{d−1}}{(d − 1)!} (z_{n} − ž_{n})^{d−1}

+\frac{(−1)^{d}t^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}

= \frac{(−1)^{d}t^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}.

Thus, the last component of H(z, t) is:

tf_{n}(z) + (1 − t)\frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}

= \frac{(−1)^{d}t^{d+1}Δ_{d}}{d!} (z_{n} − ž_{n})^{d} + (1 − t)\frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}
= \frac{(−1)^{d}(1 − t + t^{d+1})Δ_{d}}{d!} (z_{n} − ž_{n})^{d}.

When t = 0 and t = 1, 1 − t + t^{d+1} = 1. When t ∈ (0, 1), 1 − t + t^{d+1} > t^{d+1} > 0. Thus, 1 − t + t^{d+1} ≠ 0 for t ∈ [0, 1]. Then, H(z, t) = 0 implies (z_{n} − ž_{n})^{d} = 0, and consequently, z_{n} − ž_{n} = 0 or z_{n} = ž_{n}.

Now we know H(z, t) has a unique zero at (z̃_{1}, …, z̃_{n−1}, z̃_{n}) = (ž_{1} − tα_{1}(z_{n} − ž_{n}), …, ž_{n−1} − tα_{n−1}(z_{n} − ž_{n}), ž_{n}). Accordingly, H̃((x, y), t) has a unique zero at (x̃_{1}, ỹ_{1}, …, x̃_{n−1}, ỹ_{n−1}, x̃_{n}, ỹ_{n}) = (x̌_{1} − tα_{1}(x_{n} − x̌_{n}), (y̌_{1} − tα_{1}(y_{n} − y̌_{n}), …, x̌_{n−1} − tα_{n−1}(x_{n} −x̌_{n}), y̌_{n−1} −tα_{n−1}(y_{n} −y̌_{n}), x̌_{n}, y̌_{n}). Based on the way we have constructed x_{n} and y_{n}, we have

|x̃_{k} − x̌_{k}| = |tα_{k}(x_{n} − x̌_{n})| ≤ |α_{k}(x_{n} − x̌_{n})| ≤ |α_{k}|\frac{w(x_{n})}{2} < \frac{w(x_{k})}{2}
and |ỹ_{k} − y̌_{k}| = |tα_{k}(y_{n} − y̌_{n})| ≤ |α_{k}(y_{n} − y̌_{n})| ≤ |α_{k}|\frac{w(y_{n})}{2} < \frac{w(y_{k})}{2},

where k = 1, …, n − 1. Thus, x̃_{k} ∉ ∂x_{k} and ỹ_{k} ∉ ∂y_{k} for k = 1, …, n − 1. Obviously, x̃_{n} = x̌_{n} ∉ ∂x_{n} and ỹ_{n} = y̌_{n} ∉ ∂y_{n}. This implies H̃((x, y), t) ≠ 0 for (x, y) ∈ ∂x̃ and t ∈ [0, 1]. Then, by Theorem 1.4,

d(F̃, x̃, 0) = d(G̃, x̃, 0).

Next, we prove d(G̃, x̃, 0) = d. Perturb G(z) by an arbitrary small ϵ to define

G_{ϵ}(z) = (g_{1ϵ}(z), …, g_{(n−1)ϵ}(z), g_{nϵ}(z)),

where

g_{kϵ}(z) = g_{k}(z) = (z_{k} − ž_{k}) for 1 ≤ k ≤ n − 1,
(3.3)  g_{nϵ}(z) = g_{n}(z) + ϵ = \frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d} + ϵ.

Let p_{kϵ}(x, y) ≡ ℜ(g_{kϵ}(x + iy)) and q_{kϵ}(x, y) ≡ ℑ(g_{kϵ}(x + iy)). With this, define

G̃_{ϵ}(x, y) ≡ (p_{1ϵ}(x, y), q_{1ϵ}(x, y), …, p_{nϵ}(x, y), q_{nϵ}(x, y)).

It is obvious that p_{kϵ}(x, y) = x_{k}−x̌_{k} and q_{kϵ}(x, y) = y_{k}−y̌_{k} for k = 1, …, n−1. Assume ϵ is small enough. Then, G_{ϵ}(z), and thus G̃_{ϵ}(x, y), have d zeros z̃ = (z̃_{1}, …, z̃_{n−1}, z̃_{n}), or x̃ = (x̃_{1}, ỹ_{1}, …, x̃_{n−1}, ỹ_{n−1}, x̃_{n}, x̃_{n}) in x̃, with z̃_{k} − ž_{k} = 0, or x̃_{k} − x̌_{k} = 0 and ỹ_{k} − y̌_{k} = 0 for k = 1, …, n − 1, and (z̃_{n} − ž_{n})^{d} = \frac{d!ϵ}{(−1)^{d+1}Δ_{d}} ≠ 0. \frac{∂g_{nϵ}}{∂z_{n}} (z̃) =

\frac{(−1)^{d}Δ_{d}}{(d−1)!} (z̃_{n} − ž_{n})^{d−1} ≠ 0.

|\frac{∂G̃_{ϵ}}{∂x_{1}∂y_{1} … ∂x_{n}∂y_{n}} (x̃)| = | \begin{matrix}1 & 0 & … & 0 & 0 & 0 \\ 0 & 1 & … & 0 & 0 & 0 \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & … & 1 & 0 & 0 \\ 0 & 0 & … & 0 & \frac{∂p_{nϵ}}{∂x_{n}} & \frac{∂p_{nϵ}}{∂y_{n}} \\ 0 & 0 & … & 0 & \frac{∂q_{nϵ}}{∂x_{n}} & \frac{∂q_{nϵ}}{∂y_{n}} \end{matrix} |

(3.4)  = | \begin{matrix}\frac{∂p_{nϵ}}{∂x_{n}} & \frac{∂p_{nϵ}}{∂y_{n}} \\ \frac{∂q_{nϵ}}{∂x_{n}} & \frac{∂q_{nϵ}}{∂y_{n}} \end{matrix} | = | −\frac{∂p_{nϵ}}{∂y_{n}}  \frac{∂p_{nϵ}}{∂x_{n}}^{\frac{∂p_{nϵ}}{∂x_{n}}  \frac{∂p_{nϵ}}{∂y_{n}}} |

= (\frac{∂p_{nϵ}}{∂x_{n}} )^{2} + (\frac{∂p_{nϵ}}{∂y_{n}} )^{2} = |\frac{∂g_{nϵ}}{∂z_{n}} (z̃)|^{2} > 0.

Thus, by Theorem 1.1, d(G̃_{ϵ}, x̃, 0) = d, and then d(G̃, x̃, 0) = d by Theorem 1.3. Finally,

d(F̃, x̃, 0) = d(G̃, x̃, 0) = d.

Unless the components of F are exactly linear and degree d polynomials, the equalities (2.3) and (2.4) in §2 are not exact, but are only approximately true to second order. However, if second-order approximations are accurate, we can expect the degree to be equal to d. The disadvantage of the proof of Theorem 3.1 is that it does not lead to a practical computational technique. If we try to verify H(z, t) ≠ 0 or H̃((x, y), t) ≠ 0 when (x, y) ∈ ∂x̃ and t ∈ [0, 1], then it would require 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. First, we would need to compute Δ_{d}, which involves all partial derivatives of order 1 and order d. This is expensive when both n and d are large. Second, we would need to know where the solutions of u_{n}(x) = 0 and v_{n}(x) = 0 are on x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} when z_{k} = ž_{k} − tα_{k}(z_{n} − ž_{n}), and the search process for such solutions is expensive.

We could try to verify H(z, t) ≠ 0 when (x, y) ∈ ∂x̃ and t ∈ [0, 1] in another way: verify H(z, t) = 0 has a unique solution in the interior of x̃ when t ∈ [0, 1]. However, we will run into the singular situation again if we do that.

In fact, there is an alternative algorithm to compute the degree. That will be the subject of next section.

4. Algorithm. The algorithm we present here is similar to the algorithm in [10]. Based on Theorem 1.6 in §1.2, the following theorem underlies our algorithm.

Theorem 4.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 solutions, if there are any, on x\underline{n} and x\overline{n} with y_{n} in the interior of y_{n}, and F̃_{¬u_{n}} = 0 has solutions, if there are any, on y\underline{n} and y\overline{n} with x_{n} in the interior of x_{n};

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

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

condition 2. Then, for a fixed s ∈ {−1, 1},

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

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

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

−s ∑_{yn=\overline{y}_{n} ; {F̃¬u_{n} (x,y)=0 \atop sgn(un(x,y))=s}} 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 conditions 2 and 3 imply F̃ ≠ 0 on x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}. Thus, F̃ ≠ 0 on ∂x̃. Now,

9

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) (2n − 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}. Thus, by 1 and 2, F̃_{¬u_{n}} ≠ 0 on ∂x\underline{n}. Similarly, F̃_{¬u_{n}} ≠ 0 on ∂x\overline{n}, ∂y\underline{n} and ∂y\overline{n}. Thus, condition 1 in Theorem 1.6 is satisfied. Finally, with condition 4, all the conditions of Theorem 1.6 are satisfied. The formula is thus obtained.

By constructing the box x̃ according to (2.7), we can verify 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, since u_{k}(x, y) ≈ (x_{k} −x̌_{k})+α_{k}(x_{n}−x̌_{n}) ≠ 0 on x\underline{k} and x\overline{k}, and v_{k}(x, y) ≈ y_{k} + α_{k}y_{n} ≠ 0 on y\underline{k} and y\overline{k}. This 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, thus dramatically reducing the total computational cost. The four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} remaining to be searched are (2n − 1)-dimensional boxes. However, exploitation of (2.5) will reduce the search of solutions of F̃_{¬u_{n}}(x, y) = 0 on the (2n−1)-dimensional boxes to actually a one dimensional search. We use x\underline{n} as an example to explain this.

On x\underline{n}, x_{n} = \underline{x}_{n}. We know from (2.5) 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. Then, we can divide y_{n} into smaller subintervals. For a small subinterval y_{n}^{0} of y_{n}, we can formally 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 have reduced the search to searching the one dimensional interval y_{n}, much less costly than searching a (2n − 1)-dimensional box when n is large. Furthermore, if we know approximately where the solutions of F̃_{¬u_{n}}(x, y) = 0 are, we can even reduce the cost of the one dimensional search. To this end, we will next analyze the solutions of F̃_{¬u_{n}}(x, y) = 0 on the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}.

For convenience of analysis, assume (2.3) and (2.4) are exact. In other words, for 1 ≤ k ≤ (n − 1),

\begin{matrix}u_{k}(x, y) & = & (x_{k} − x̌_{k}) + α_{k}(x_{n} − x̌_{n}), \\ v_{k}(x, y) & = & y_{k} + α_{k}y_{n}. \end{matrix} \} (4.1)

and

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}} − ž_{k_{1}})(z_{k_{2}} − ž_{k_{2}}) + … +

\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}} − ž_{k_{1}}) … (z_{k_{d}} − ž_{k_{d}}).

From this, F̃_{¬u_{n}}(x, y) = 0 implies u_{k}(x, y) = 0 and v_{k}(x, y) = 0, i.e. f_{k}(z) = 0, whence z_{k} = ž_{k} − α_{k}(z_{n} − ž_{n}) for k = 1, …, n − 1. In terms of the real and imaginary coordinates,

\begin{matrix}x_{k} & = & x̌_{k} − α_{k}(x_{n} − x̌_{n}), \\ y_{k} & = & −α_{k}y_{n}, \end{matrix} \} (4.2)

for k = 1, …, n − 1. Plugging z_{k} − ž_{k} = −α_{k}(z_{n} − ž_{n}), k = 1, …, n − 1, into f_{n}(z) as in the proof of Theorem 3.1 then gives

f_{n}(z) = \frac{(−1)^{d}Δ_{d}}{d!} (z_{n} − ž_{n})^{d}.

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

Fig. 4.1. When d is odd. Here, d = 3. v_{n} = 0 on solid lines and u_{n} = 0 on dashed lines. The thick dots are the solutions of F̃_{¬u_{n}}(x, y) = 0 on ∂x̃.

Thus, u_{n}(x, y) = ℜ(f_{n}(z)) = \frac{(−1)^{d}Δ_{d}}{d!} ℜ((z_{n} − ž_{n})^{d}) and v_{n}(x, y) = ℑ(f_{n}(z))

= \frac{(−1)^{d}Δ_{d}}{d!} ℑ((z_{n} − ž_{n})^{d}). Setting z_{n} − ž_{n} = r(cos(θ) + i sin(θ)), we obtain u_{n}(x, y) =

\frac{(−1)^{d}Δ_{d}}{d!} r cos(dθ) and v_{n}(x, y) = \frac{(−1)^{d}Δ_{d}}{d!} r sin(dθ), so u_{n}(x, y) = 0 is equivalent to cos(dθ) = 0 and v_{n}(x, y) = 0 is equivalent to sin(dθ) = 0. If we choose x_{n} and y_{n} such that

(4.3)  \begin{matrix}w & (y_{n}) \\ w & (x_{n}) \end{matrix} = tan (\frac{π}{4d}) ,  that is,  w(y_{n}) = tan (\frac{π}{4d}) w(x_{n}),

then all solutions of v_{n}(x, y) = 0, and consequently all solutions of F̃_{¬u_{n}}(x, y) = 0 are arranged in a known pattern on x\underline{n}, x\overline{n}, y\underline{n}, and y\overline{n}. In particular, on x\underline{n}, x̃_{n} = \underline{x}_{n}. v_{n}(x, y) = 0 has a unique solution ỹ_{n} = 0. Plugging these into (4.2), we get the unique solution of F̃_{¬u_{n}}(x, y) = 0 with

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

Similarly, F̃_{¬u_{n}}(x, y) = 0 has a unique solution on x\overline{n} with

(x̃, ỹ) = (x̌_{1} − α_{1}(\overline{x}_{n} − x̌_{n}), 0, …, x̌_{n−1} − α_{n−1}(\overline{x}_{n} − x̌_{n}), 0, \overline{x}_{n}, 0).
On y\underline{n}, ỹ_{n} = \underline{y}_{n}. v_{n}(x, y) = 0 has d − 1 solutions with

x̃_{n} = \begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix}, m = d − 1, d − 2, …, 1. (4.4)

Plugging these into (4.2) gives the d − 1 solutions (x̃, ỹ) of F̃_{¬u_{n}}(x, y) = 0 with

(x̃, ỹ) = (x̌_{1} − α_{1}(\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}), α_{1}\underline{y}_{n}, …, x̌_{n−1} − α_{n−1}(\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}) ,

−α_{n−1}\underline{y}_{n}, x̌_{n} − (\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}), \underline{y}_{n}) .

Similarly, F̃_{¬u_{n}}(x, y) = 0 has d − 1 solutions on y\overline{n} with

(x̃, ỹ) = (x̌_{1} − α_{1}(\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}), …, α_{n−1}(\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}) ,

−α_{n−1}\overline{y}_{n}, x̌_{n} − (\begin{matrix}w & (y_{n}) \\ tan & (\frac{mπ}{d} ) \end{matrix} − x̌_{n}), \overline{y}_{n}) .

For example, Figure 4.1 gives the solutions of v_{n}(x, y) = 0 on the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} when d = 3.

To use the above analysis to find approximations to the solutions of F̃_{¬u_{n}} = 0 on the faces we search, we need to know d; we present a heuristic for d in the next section.

Now, we present our algorithm. The algorithm consists of three phases:

1. the box-construction phase where we set x̃,

2. the elimination phase where we use interval evaluations to 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, and thus eliminate those 4n − 4 faces, and

3. the search phase, where we (a) search x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} to locate the solutions of F̃_{¬u_{n}}(x, y) = 0,

(b) compute the signs of u_{n} and determinants of the Jacobi matrices of F̃_{¬u_{n}} at those solutions,

(c) compute the degree contributions of each of the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} according to Theorem 4.1, and

(d) finally sum up to get the degree.

ALGORITHM 1 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 width of x_{n} as in (2.7).

4. Set the width of y_{n} to be the minimum of that obtained from conditions (2.7) and (4.3).

Elimination Phase Do for 1 ≤ k ≤ n − 1

1. Do for x\underline{k} and x\overline{k} (a) Compute the mean-value extension of u_{k} over that face.

(b) If 0 ∈ u_{k}, then stop and signal failure.

2. Do for y\underline{k} and y\overline{k} (a) Compute the mean-value extension of v_{k} over that face.

(b) If 0 ∈ v_{k}, then stop and signal failure. Search Phase

1. Set the value of s ∈ {+1, −1}.

(a) Initialize s to be +1. Initialize search lower and search upper to be false. (See the second note below.)

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

and thus to get a subface x0 ; \underline{n} (or x0 ; \overline{n}) of x\underline{n} (or x\overline{n}.). ii. If x̃_{k} ∩ x_{k} = ∅, then cycle. iii. Compute the mean-value extension u_{n} over x0 ; \underline{n} (or x0 ; \overline{n}). iv. If u_{n} contains 0, then set search lower (or search upper) to be true

and cycle. v. If u_{n} does not contain 0, then set s = −sgn(u_{n}).

(c) If u_{n} does not contain 0 on both x\underline{n} and x\overline{n}, then set s to be the opposite sign to the sign of u_{n} on x\overline{n}, and if u_{n} has different signs on x\underline{n} and x\overline{n}, then set search lower to be true.

2. For x\underline{n} (or x\overline{n}), if search lower (or search upper) is true, then take x\underline{n} and 0 as inputs and apply Algorithm 2 to compute the degree contribution of x\underline{n} (or x\overline{n}).

3. For y\underline{n} (or y\overline{n})

(a) Use (4.4) to compute the x̃_{n}^{m}, m = d − 1, d − 2, …, 1, x̃_{n}^{d−1} < x̃_{n}^{d−2} < … < x̃_{n}^{1}, corresponding to the d − 1 approximate solutions of F̃_{¬u_{n}} = 0 on y\underline{n}.

(b) Divide x_{n} into d − 1 parts x_{n}^{m}, m = 1, …, d − 1 as follows:

x_{n}^{1} = [\underline{x}_{n}, (x̃_{n}^{1} + x̃_{n}^{2})/2], x_{n}^{m} = [(x̃_{n}^{m−1} + x̃_{n}^{m})/2, (x̃_{n}^{m} + x̃_{n}^{m+1})/2]

for m = 2, …, d − 2, and x_{n}^{d−1} = [(x̃_{n}^{d−2} + x̃_{n}^{d−1})/2, \overline{x}_{n}].

(c) Do for m = 1, …, d − 1.

i. Set a subface ym ; \underline{n} of y\underline{n} (or ym ; \overline{n} of y\overline{n}) by replacing x_{n} by x_{n}^{m}. ii. Apply Algorithm 3 with ym ; \underline{n} and x̃_{n}^{m} as inputs, to compute the degree

contribution of ym ; \underline{n} (or ym ; \overline{n}.)

(d) Add the degree contributions in the last step to get the degree contri-

bution of y\underline{n} (or y\overline{n}.)

4. Add the degree contributions of x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} to get the overall degree. END OF ALGORITHM 1

Notes for Algorithm 1

1. In Step 3 of the box-setting phase, the width w(x_{n}) of x_{n} depends on the accuracy of the approximate solution x̌ of the system F (x) = 0: w(x_{n}) should be much larger than |x̌_{k} − x_{k}^{∗}|, but also should be sufficiently small for a quadratic model to be accurate over the box.

2. We may set s to minimize the amount of work required to evaluate the sum in Theorem 4.1. In particular, if we know sgnu_{n} = σ on a large number of faces, then setting s = −σ will eliminate the need to search those faces.

ALGORITHM 2 Inputs: x\underline{n} and y̌_{n} (or x\overline{n} and y̌_{n})

1. (a) 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.

(b) If x̃_{k} ∩ x_{k} = ∅, then return the degree contribution of that face as 0.

(c) Update x_{k}.

2. (a) Compute the mean-value extension u_{n} over that face.

(b) If s × sgn(u_{n}) < 0, then return the degree contribution of that face as 0.

3. Construct a small subinterval y_{n}^{0} of y_{n} centered at y̌_{n}.

4. Step 4 to step 9 are identical to step 1(d) to step 1(i), respectively, of the search phase in the algorithm in [10].

10. Apply Theorem 4.1 to compute the degree contribution of x\underline{n} or x\overline{n}. END OF ALGORITHM 2

ALGORITHM 3 Inputs: y\underline{n} and x̌_{n} (or y\overline{n} and x̌_{n}.)

1. (a) Use mean-value extensions for v_{k}(x, y) = 0 to solve for y_{k} to get sharper

bounds ỹ_{k} with width O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1.

(b) If ỹ_{k} ∩ y_{k} = ∅, then return the degree contribution of that face as 0.

(c) Update y_{k}.

2. (a) Compute the mean-value extension u_{n} over that face.

(b) If s × sgn(u_{n}) < 0, then return the degree contribution of that face as 0.

3. Construct a small subinterval x_{n}^{0} of x_{n} which is centered at x̌_{n}.

4. Step 4 to step 9 are identical to step 2(d) to step 2(i), respectively, of the search phase in the algorithm in [10].

10. Same as step 10 of Algorithm 2. END OF ALGORITHM 3

Notes for Algorithm 2 and Algorithm 3

1. Algorithms 2 and 3 are identical to steps 1 and 2 of the search phase of the algorithm in [10], except, in Algorithm 2, y̌_{n} can be any interior point of y_{n}, while y̌_{n} is assumed to equal zero in step 1 of the search phase in the algorithm in [10]. Similarly, in Algorithm 3, x̌_{n} can be any interior point of x_{n}, whereas x̌_{n} is assumed to equal the center of x_{n} in step 2 of the search phase in the algorithm in [10].

2. In the overall algorithm, Algorithm 1, the actual inputs are ym ; \underline{n} and x̃_{n}^{m} when Algorithm 3 is applied. However, for notational simplicity, we use y\underline{n} and x̌_{n} as inputs in the presentation of Algorithm 3.

The computational complexity of Algorithms 1, 2, and 3 is O (n^{3}). (See [10] for detailed analysis.) Thus, the computational complexity of the overall algorithm, Algorithm 1, is O (n^{3}). This is the best possible order, 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}).

5. A Heuristic for the Degree. The algorithms in §4 require a value for d to locate the approximate positions of solutions of F̃_{¬u_{n}} = 0 on the faces we search. Here, we present a practical heuristic for the value of d.

Suppose (2.3) and (2.4) are exact. Then, if we set x_{k} − x̌_{k} = −α_{k}(x_{n} − x̌_{n}), k = 1, …, n and plug those equalities into f_{n}, we obtain a univariate function

g(x_{n} − x̌_{n}) = f_{n}(x_{1}, …, x_{n})
= \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)(x_{k_{1}} − x̌_{k_{1}})(x_{k_{2}} − x̌_{k_{2}}) + … +

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

= \frac{1}{2!} ∑_{k_{1}=1}^{n} ∑_{k_{2}=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)(−1)^{2}α_{k_{1}}α_{k_{2}}(x_{n} − x̌_{n})^{2} + … +

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

= \frac{(−1)^{2}}{2!} (x_{n} − x̌_{n})^{2} \frac{∂^{2}f_{n}}{∂x_{k_{1}}∂x_{k_{2}}} (x̌)α_{k_{1}}α_{k_{2}} + …

+\frac{(−1)^{d}}{d!} (x_{n} − x̌_{n})^{d} \frac{∂^{2}f_{n}}{∂x_{k_{1}} … ∂x_{k_{d}}} (x̌)α_{k_{1}} … α_{k_{d}}

= \frac{(−1)^{2}Δ_{2}}{2!} (x_{n} − x̌_{n})^{2} + … + \frac{(−1)^{d−1}Δ_{d−1}}{(d − 1)!} (x_{n} − x̌_{n})^{d−1}

+\frac{(−1)^{d}Δ_{d}}{d!} (x_{n} − x̌_{n})^{d}

= \frac{(−1)^{d}Δ_{d}}{d!} (x_{n} − x̌_{n})^{d} = \frac{Δ_{d}}{d!} (x̌_{n} − x_{n})^{d}.

Setting

K(r, x_{n} − x̌_{n}) ≡ \frac{g(x_{n} − x̌_{n})}{(x_{n} − x̌_{n})^{r}} = \frac{Δ_{d}}{d!} (x̌_{n} − x_{n})^{d−r},

it is clear that K(d, x_{n} − x̌_{n}) = Δ_{d}/d! is independent of x_{n}, while K(r, x_{n} − x̌_{n}) depends on x_{n} for any other r value. We have the following ratios.

\frac{K(d, δ(x_{n} − x̌_{n}))}{K(d, x_{n} − x̌_{n})} = \frac{\frac{Δ_{d}}{d!}}{\frac{Δ_{d}}{d!}} = 1, while

R(r) = \frac{K(r, δ(x_{n} − x̌_{n}))}{K(r, x_{n} − x̌_{n})} = \frac{\frac{Δ_{d}}{d!} (δ(x̌_{n} − x_{n}))^{d−r}}{\frac{Δ_{d}}{d!} (x̌_{n} − x_{n})^{d−r}} = δ^{d−r}

for any other r value. The first ratio R(d) always equals 1, but R(r), r ≠ d, depends on the δ value. We can choose δ to distinguish d from other r values. For example, if we choose δ = 100, then R(r) is not smaller than 100 when r is smaller than d and is not larger than 0.01 when r is larger than d. Both values are sufficiently different from 1. We can also vary the δ value to check our detection of d. Thus, R(r) is a good heuristic to determine the value of d.

The above discussion is based on the assumptions in §2. However, unless the first n − 1 components of F are exactly linear and the last component is a degree-d polynomial of n variables, those assumptions are only approximately true. There are some finer issues to consider. The above analysis is valid only when the equality g(x_{n} − x̌_{n}) = \frac{Δ_{d}}{d!} (x̌_{n} − x_{n})^{d} is accurately approximated. That implies (x̌_{n} − x_{n})^{d} should dominate the value of g(x_{n} − x̌_{n}). Actually,

g(x_{n} − x̌_{n}) = ∑_{k=1}^{d−1} c_{k}Δ_{k}(x_{n} − x̌_{n})^{k} + c_{d}Δ_{d}(x_{n} − x̌_{n})^{d} + ∑_{k=d+1}^{∞} c_{k}Δ_{k}(x_{n} − x̌_{n})^{k},

where, approximately, Δ_{1} = … = Δ_{d−1} = 0, Δ_{d} ≠ 0. Thus, x_{n} − x̌_{n} and δ(x_{n} − x̌_{n}) should not be too small, since ∑_{k=1}^{d−1} c_{k}Δ_{k}(x_{n} − x̌_{n})^{k} could dominate otherwise. They should not be too big either, since ∑_{k=d+1}^{∞} c_{k}Δ_{k}(x_{n} −x̌_{n})^{k} could dominate otherwise. If Δ_{k} ≈ 0, k = 1, …, d − 1 are quite accurate, then we can choose x_{n} − x̌_{n} very small, so both ∑_{k=1}^{d−1} c_{k}Δ_{k}(x_{n} − x̌_{n})^{k} and ∑_{k=d+1}^{∞} c_{k}Δ_{k}(x_{n} − x̌_{n})^{k} can be ignored in the detection of d.

The choice of x_{n} − x̌_{n} is independent of the settings of x_{k}, k = 1, …, n, since we only want to know what d is at that stage.

An alternative choice for detecting d is to compute the values of Δ_{k}, k = 1, 2, … by interval evaluations until we get some Δ_{k_{0}} that is sufficiently different from 0. Then, we can decide d = k_{0}. The obvious disadvantage of this method is that it’s too expensive for just detecting the value of d, since computation of Δ_{k} involves computations of all k-th order derivatives. Furthermore, even if we actually evaluate Δ_{k}, k = 1, 2, …, spending much time in the process, we still can not detect the value of d if the magnitudes of Δ_{k}, k = 1, …, d − 1, d, are not sufficiently different either due to the problem itself or due to the range overestimation in interval computations.

6. Numerical Results. In this section, we present numerical results for the algorithm in §4.

6.1. Test Problems. Example 1. (The same as Example 3 from [10], motivated from considerations in [5] 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.

In Example 1, if we change the exponent of x from 2 to 3, then we get another problem.

Example 2. This example is identical to Example 1, except that we set f(x) = h(x, t) = (1 − t)(Ax − x^{3}) − tx.

The test set consists of Example 1 and Example 2 with n = 5, 10, 20, 40, 80, 160. For each test problem, we used (0, 0, …, 0), the exact solution to F (x) = 0, as the approximate solution to the problem F (x) = 0. We set the widths w(x_{k}) and w(y_{k}) to 10^{−2} for 1 ≤ k ≤ n − 1, then the algorithm automatically computed w(x_{n}) and w(y_{n}). For all the problems, the algorithm succeeded.

6.2. Test Environment. The algorithm in §4 was programmed in the Fortran 90 environment developed and described in [8, 9]. Similarly, all the test functions were programmed using the same Fortran 90 system, which generated internal symbolic representations of the functions to execution of the numerical tests. In the actual tests, generic routines then interpreted the internal representations to obtain both floating point and interval values.

The Sun Fortran 95 compiler version 6.0 was used on a Sparc Ultra-1 model 140 with optimization level 0. Execution times were measured with the Port library routine ETIME. All times are given in CPU seconds.

6.3. Test Results. We present the numerical results in Table 6.1.

The column labels of Table 6.1 are as follows. Problem: names of the problems identified in §6.1. n: number of independent variables. Heuristic Degree: the heuristic value of the degree computed by the heuristic de-

scribed in §5. Success: whether the algorithm was successful. Verified Degree: topological degree verified by the algorithm. CPU Time: CPU time in seconds of the algorithm. Time Ratio: The ratio of two successive CPU times for a particular problem.

We can see from Table 6.1 that the algorithm was successful on all the problems in the test set. We also see from the CPU time ratios that the algorithm is approximately of order O (n^{3}) in practice. However, when the degree is higher, the system F (x) is

Table 6.1 Numerical Results

Heuristic  Verified
Problem  n Degree Success Degree CPU Time Time Ratio
Example 1  5  2  Yes  2  1.13
Example 1 10  2  Yes  2  5.99  5.30
Example 1 20  2  Yes  2  38.40  6.41
Example 1 40  2  Yes  2  273.61  7.13
Example 1 80  2  Yes  2  2198.14  8.03
Example 1 160  2  Yes  2  13033.22  5.93
Example 2  5  3  Yes  3  39.27
Example 2 10  3  Yes  3  10.31  0.26
Example 2 20  3  Yes  3  74.32  7.21
Example 2 40  3  Yes  3  481.23  6.48
Example 2 80  3  Yes  3  3805.06  7.91
Example 2 160  3  Yes  3  33944.20  8.92

flatter at the singular solution. We can expect that the condition number of the Jacobi matrix of the system F̃_{¬u_{n}} will be larger, and thus, make the interval Newton method method to verify the unique solutions of F̃_{¬u_{n}} in Step 5 of Algorithm 2 and Step 5 of Algorithm 3 less efficient: More iterations should be expected when the condition number is larger. The experimental results are consistent with our expectations.

REFERENCES

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

[2] P. S. Aleksandrov and H. Hopf. Topologie. Springer, Berlin, 1935.

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

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

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

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

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

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

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

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

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

[12] J. M. Ortega and W. C. Rheinboldt. Iterative Solution of Nonlinear Equations in Several Variables. Academic Press, New York, 1970.

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

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