Root finding with interval arithmetic

Many books and articles explain how to solve the nonlinear equation ${f}\!\left(x\right)=0$ for $x\in\mathds{R}{}$, and some of them consider the verified solution of such equation, that is, finding solutions with rigorous bounds on them. Here we discuss the computation of verified solutions using interval arithmetic. This is an old subject, but the interval arithmetic literature is mostly concerned with the solution of the general equation ${f}\!\left(x\right)=0$ for $x\in{\mathds{R}}^{n}$ using variations of Newton’s method, since its very beginning [2, 10, 11]. For instance, the third chapter of [13] gives a nice description of interval Newton’s method and [2, 3] present interesting ways to improve it. Here we focus on the simplest case $n=1$, because it is very important in practice and there are techniques which are applicable only for $n=1$. As far as we know, the most detailed discussion of this particular case is presented in chapter 9 of [4], but we believe that there is more to be said about this subject than what is presented there.

In this article, $f$ is a function from $\mathds{R}{}$ to $\mathds{R}{}$. We denote $f$’s $\ell$th derivative by $f^{\left({\ell}\right)}$ and when we mention such derivatives we assume that they exist. We consider the class

of all closed intervals, including the empty and unbounded ones. For ${\cal{S}}\subset\mathds{R}{}$, we write

We assume that we have extensions of $f$ and its derivatives, in the sense that:

We usually identify the point $t\in\mathds{R}{}$ with the interval $[t,t]$, and for functions $\mathbb{f}:{\mathds{I}}\rightarrow{\mathds{I}}$ we write ${\mathbb{f}}\!\left(t\right):={\mathbb{f}}\!\left([t,t]\right)$. The set of roots is denoted by ${\cal{R}}$, and $r$ is a generic root.

When $f$ is differentiable, the main tool for root finding is Newton’s method

This method has a natural extension to interval arithmetic:

where $t_{k}$ is a point in $\mathbb{x}_{k}$ and $f$ and $f^{\prime}$ are extensions of $\mathbb{f}$ and $\mathbb{f}^{\prime}$. Traditionally, we compute $t_{k}$ by rounding the midpoint of $\mathbb{x}_{k}$. The first question an alert numerical analyst would ask about Equation $\left({\ref{newton_i}}\right)$ is:

However, we see few people asking the following question:

Both questions are quite relevant, but before we answer them we must say that the interval version of Newton’s method may be implemented without Equation $\left({\ref{newton_i}}\right)$. Instead, we could write $f$ in the centered form

and replace the extended derivative ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$ by an extension ${\mathbb{s}_{t_{k}}}\!\left(\mathbb{x}_{k}\right)$ of the slope $s_{t_{k}}$ of $f$ at $t_{k}$. This leads to an improved version of Newton’s method given by

The slope is defined as

and centered forms are already mentioned in Moore’s book [10]. They have been discussed in detail in the interval arithmetic literature [7, 12] and have generalizations called Taylor forms or Taylor models [9]. Our algorithm uses only the plain centered form $\left({\ref{center}}\right)$, and in situations in which practical implementations of such generalizations are available we would use them only as a tool to obtain more accurate centered forms.

The Mean Value Theorem shows that

and this implies that any extension of $\mathbb{f}^{\prime}$ is an extension of $s_{t_{k}}$, but there may be better ones, specially when the interval $\mathbb{x}_{k}$ is not small. For instance, if ${f}\!\left(t\right):=t^{2}$ then

and any extension of $f^{\prime}$ yields intervals twice as large than necessary for an slope.

In practice the gain by replacing derivatives by slopes may not be great, because usually ${f}^{\prime}\!\left(r\right)\neq 0$ at the roots $r$ and in this case, for $t_{k}$ close to $r$, the difference between $s_{t_{k}}$ and ${f}^{\prime}\!\left(t\right)$ for $t$ in a short interval $\mathbb{x}_{k}$ is not large. Moreover, extensions of derivatives have an important feature that extensions of slopes do not have: derivative extensions can detect monotonicity. In other words, if ${\mathbb{f}^{\prime}}\!\left(\mathbb{x}\right)\subset(0,+\infty)$ then $f$ is increasing in $\mathbb{x}$, but we cannot reach the same conclusion from the fact that ${s_{t_{k}}}\!\left(\mathbb{x}\right)\subset(0,+\infty)$. The ease with which we can use information about monotonicity when we have only one variable is what makes the case $n=1$ so special, and this is the reason why this article was written to start with.

We can ask questions similar to $\left({\ref{qd0}}\right)$ and $\left({\ref{qd1}}\right)$ about the modified Newton’s step $\left({\ref{newton_t1}}\right)$:

and there are at least two more questions we must ask:

After much thought about the questions above, we devised a root finding algorithm which combines the interval versions of Newton’s method in Equations $\left({\ref{newton_i}}\right)$ and $\left({\ref{newton_t1}}\right)$, and modifies them in order to exploit montonicity. Our algorithm tries to strike a balance between theory and practice, and takes into account the following practical issues:

1. 1.

   Usually, the evaluation of interval extensions over intervals is much less accurate than the evaluation of interval extensions at points, that is, the width of the interval ${\mathbb{f}}\!\left(t_{k}\right)$ is much smaller than the width of the interval ${\mathbb{f}}\!\left(\mathbb{x}_{k}\right)$. The same applies to the derivative $\mathbb{f}^{\prime}$.

2. 2.

   Usually, the floating point evaluation of $d_{k}={f}^{\prime}\!\left(t_{k}\right)$ yields a reasonable estimate this derivative, at a much lower cost than the evaluation of the extensions ${\mathbb{f}}^{\prime}\!\left(t_{k}\right)$ or ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$. The only defect of $d_{k}$ is that it does not come with a guarantee of its accuracy. As a result we can, and should, use floating point computations in order to obtain reasonable estimates (which may turn out to be inaccurate sometimes) and resort to interval computation only when we absolutely need guarantees about our results. In particular, the computed ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$ may be very wide even when the floating point $d_{k}={f}^{\prime}\!\left(t_{k}\right)$ would lead to a good Newton step $t_{k+1}=t_{k}-\hat{w}_{k}/d_{k}$, where $\hat{w}_{k}$ is the mid point of $\mathbb{w}_{k}={\mathbb{f}}\!\left(t_{k}\right)$.

3. 3.

   In interval arithmetic, our goal is to find short intervals $\mathbb{r}$ which may contain roots and to discard intervals guaranteed not to contain roots.

4. 4.

   The simplest way to ensure that an interval $\mathbb{r}=[\underline{r},\overline{r}]$ contains a root is to prove that ${f}\!\left(\underline{r}\right)$ and ${f}\!\left(\overline{r}\right)$ have opposite signs. We can do that by evaluating $\mathbb{f}$ at the points $\underline{r}$ and $\overline{r}$, and as we noted above this evaluation tends do yield sharp results.

5. 5.

   The simplest way to ensure that an interval $\mathbb{r}$ does not contain a root is to prove that ${f}\!\left(\underline{r}\right)$ and ${f}\!\left(\overline{r}\right)$ are different from zero and have the same sign and ${f}^{\prime}\!\left(\mathbb{r}\right)$ does not contain $0$. In practice, this is not as easy as in the previous item because the computed ${\mathbb{f}}^{\prime}\!\left(\mathbb{r}\right)$ tends to be inflated by the usual weaknesses of interval arithmetic. However, when we know that $f$ is monotonic we can dispense with ${\mathbb{f}}^{\prime}\!\left(\mathbb{r}\right)$ and check only the values of $f$ at the end points. This simplifies things immensely for monotonic functions. Actually, it makes little practical sense to compute ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$ once we already know that it does not contain 0, and we have this information for all nodes below a first node $\mathbb{x}_{k}$ in the branch and bound tree which is such that ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$ does not contain 0.

6. 6.

   Multiple roots, i.e., roots $r$ such that ${f}^{\prime}\left(r\right)=0$, are a major problem and should be handled differently from single roots. Usually, it is hopeless to try to get them with the same accuracy as single roots, and the cases in which we have some hope are too specific to deserve attention in a generic software.

Combining the items above, we propose an algorithm which can be outlined as follows. For each candidate interval we keep seven additional bits of information:

• •

  The sign $\sigma_{i}\in\left\{{-1,0,1}\right\}$ of $f$ at its infimum,

• •

  The sign $\sigma_{s}\in\left\{{-1,0,1}\right\}$ of $f$ at its supremum

• •

  The sign $\sigma_{d}\in\left\{{-1,0,1}\right\}$ of $f^{\prime}$ in $\mathbb{x}$

• •

  A point $t$ in its interior, indicating where it should be split in the Newton step

• •

  A point $\tilde{t}$, which lies near $\mathbb{x}$.

• •

  The value $\tilde{d}={f}^{\prime}\!\left(\tilde{t}\right)$.

• •

  The expected sign $\sigma_{t}$ of ${f}\!\left(t\right)$.

Regarding the signs above, $-1$ means definitely negative, $+1$ means definitely positive and $0$ means that we do not know the corresponding sign.

We then procedure in the usual branch and bound way:

1. 1.

   If the stack is empty then we are done. Otherwise we pop a tuple

   $$\left({\mathbb{x}_{k},\sigma_{i},\sigma_{s},\sigma_{d},t_{k},\tilde{t},\tilde{d},\sigma_{t}}\right)$$

   from the stack.

2. 2.

   If the width of $\mathbb{x}$ is below a given tolerance $\tau_{x}$ then we insert $\left({\mathbb{x},\sigma_{i},\sigma_{s},\sigma_{d}}\right)$ in a list of possible solutions and go to item 1.

3. 3.

   We compute ${\mathbb{f}}\!\left(\mathbb{x}\right)$. If it does not contain $0$, then drop $\mathbb{x}$ and go to item 1. (Tolerances are discussed in Section 3.)

4. 4.

   We compute the slope $\mathbb{s}_{k}:={\mathbb{s}_{t_{k}}}\!\left(\mathbb{x}_{k}\right)$. If it does not contain $0$ then we compute the derivative $\mathbb{d}_{k}:={\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$. If $\mathbb{d}_{k}$ does not contain $0$ then we change to the specific algorithm to handle strictly monotonic functions described in Section 5, and once we are done we go to item 1. Otherwise, we replace $\mathbb{s}_{k}$ by $\mathbb{s}_{k}\cap\mathbb{d}_{k}$ and continue.

5. 5.

   We compute ${\mathbb{f}}\!\left(t_{k}\right)$ and check whether it contains zero or is too close to zero. If it does then use the algorithm described in Section 3 to handle $\mathbb{x}_{k}$ and go to item 1.

6. 6.

   If ${\mathbb{f}}\!\left(t_{k}\right)$ is not too close to zero then and apply the version of the interval Newton step described in Section 2, obtaining at most two intervals $\mathbb{x}_{k+1}$ and $\mathbb{x}_{k+1}^{\prime}$. If there are no $\mathbb{x}_{k+1}$ then we drop $\mathbb{x}_{k}$ and go to item 1.

7. 7.

   If there is just one $\mathbb{x}_{k+1}$ then we check whether the sign $\sigma_{t}$ matches the sign of ${\mathbb{f}}\!\left(\mathbb{x}_{k}\right)$. If it does not then we set $\sigma_{t}$ to zero, take $t_{k+1}$ as the mid point of $\mathbb{x}_{k+1}$, obtain the extra information for $\mathbb{x}_{x+1}$ and push it on the stack and go to item 1. Otherwise, we compute $d={f}^{\prime}\!\left(t_{k}\right)$ using floating point arithmetic and use $\tilde{d}$ and $\tilde{t}$ to check whether the corrected Newton step $t_{k+1}$ described in Section 4 lies in $\mathbb{x}_{k+1}$. If it does then we use this $t_{k+1}$, otherwise we take $t_{k+1}$ as the midpoint of $\mathbb{x}_{k+1}$. We then obtain the additional information for $\mathbb{x}_{k+1}$, push it on the stack and go to item 1.

8. 8.

   If there are two $\mathbb{x}_{k+1}$s, then we check whether the width of $\mathbb{x}_{k}$ lies below the cluster threshold $\tau_{c}$. If it does then we deem $\mathbb{x}_{k}$ to contain a cluster of zeros, and proceed as in item 2. Otherwise we apply the same procedure as in item 7 to $\mathbb{x}_{k+1}$ and $\mathbb{x}_{k+1}^{\prime}$ and go to item 1.

We implemented the algorithm outlined above in C++, using our Moore interval arithmetic library [8]. This code is available with the arxiv version of this article. It is distributed under the Mozilla Public License 2.0. Unfortunately, the code is much more involved than the summary above, because it must deal with many technical details which were not mentioned in this article in order not to make it longer and more complex than it already is.

In the rest of the article we discuss in more detail several aspects of the algorithm outline above. We start with Section 2, in which we present a version of Newton’s method for interval arithmetic. This version is similar to the ones found in the literature, but it is slightly different because it ensures that $f$ is different from zero at some of the extreme points of the new intervals, and computes the signs of $f$ at these extremes at a low extra cost. As a result, our algorithm yields not only intervals containing the roots but also the sign of $f$ at the extreme of such intervals, and these signs certify the existence of a root in the interval when they are different from each other. In Section 3 we make some comments about interval arithmetic in general and discuss the thorny subject of tolerances, which are unavoidable for defining stopping criterion for interval root solvers. Section 4 presents yet another version of the classic Newton’s method for exact, “point arithmetic.” This “point version” is the motivation for the interval version of Newton’s method for monotone functions presented in Section 5. Finally, in Section 6 we discuss the important subject of testing.

This section is about the interval version of Newton’s step

where $\mathbb{d}_{k}=[\underline{d}_{k},\overline{d}_{k}]$ can be either the derivative ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$ or the slope ${\mathbb{s}_{t_{k}}}\!\left(\mathbb{x_{k}}\right)$. Here we make the simplifying assumption that the interval

does not contain $0$, answering to the questions $\left({\ref{qd0}}\right)$, $\left({\ref{qd1}}\right)$, $\left({\ref{qs0}}\right)$ and $\left({\ref{qs1}}\right)$ in the introduction in this case. By replacing $f$ by $-f$ if necessary, we can assume that $\overline{w}_{k}<0$, and we make this simplifying assumption from now on.

The answer to questions $\left({\ref{qd0}}\right)$ and $\left({\ref{qs0}}\right)$, in which case $\mathbb{d}_{k}$ contains $0$, is illustrated in Figure 2. In this Figure, the inclined lines have equations

and by intersecting these lines with the axis $w=0$ we obtain a better bound on ${\cal{R}}\cap\mathbb{x}_{k}$ (recall that ${\cal{R}}$ is the set of roots.) There are three possibilities for our bounds on ${\cal{R}}\cap\mathbb{x}_{k}$ after we take in to account the intersections above. We may find that:

• •

  ${\cal{R}}\cap\mathbb{x}_{k}=\emptyset$. In this case we drop $\mathbb{x}_{k}$.

• •

  ${\cal{R}}\cap\mathbb{x}_{k}$ is contained in a single interval $\mathbb{r}=[\underline{r},\overline{r}]$, as in the second and third cases in Figure 1. In this case we take $\mathbb{x}_{k+1}=\mathbb{r}$.

• •

  ${\cal{R}}\cap\mathbb{x}_{k}$ is contained in the union of two intervals $\mathbb{r}_{1}$ and $\mathbb{r}_{2}$ such that $\overline{r}_{1}\leq\underline{r}_{2}$, as in the last case of Figure 1.

There is nothing new in our approach up to this point. We just explained with a picture what most people do algebraically. Our slight improvement comes in the way we compute the intersections in Equations $\left({\ref{inter1}}\right)$. These intersections are

and, as usual in the Moore Library [8], we propose that $\overline{r}_{1}$ and $\underline{r}_{2}$ be computed with the rounding mode upwards, using the following expressions:

This order of evaluation ensure that $\overline{r}_{1}$ and $\underline{r}_{2}$ are rounded in the correct direction, so that we do not risk loosing roots. The expressions in Equation $\left({\ref{computedr1}}\right)$ and $\left({\ref{computedr2}}\right)$ are so simple that we find whether $\overline{r}_{1}$ and $\underline{r}_{2}$ where computed exactly by checking whether

and

If we find that $\overline{r}_{1}$ was computed exactly then we increase it to the next floating point number. If $\underline{r}_{2}$ was computed exactly then we decrease it to the previous floating point number. By doing so, we ensure that ${f}\!\left(\overline{r}_{1}\right)<0$ and ${f}\!\left(\underline{r}_{2}\right)<0$, without computing ${\mathbb{f}}\!\left(\overline{r}_{1}\right)$ or ${\mathbb{f}}\!\left(\underline{r}_{2}\right)$. We should also mention that it is possible to prove that, even after the rounding, incrementing and decrementing steps above, $\overline{r}_{1}\leq t_{k}\leq\underline{r}_{2}$ and there is no risk of $\overline{r}_{1}$ crossing $\underline{r}_{2}$.

Regarding the cost of all this, note that in the Moore library the rounding mode is always upwards. Therefore, there is no cost associated with changing rounding modes. Moreover, the expensive operations in Equations $\left({\ref{computedr1}}\right)$ and $\left({\ref{computedr2}}\right)$ are the divisions, and the extra couple of sums and multiplications do not increase the overall cost of computing the intersections by much. In fact, if we take into account branch prediction and speculative execution and the fact that computations are usually inexact, the extra cost due to the verification in Equations $\left({\ref{veri1}}\right)$ and $\left({\ref{veri2}}\right)$ and the decrement of $\overline{r}_{1}$ and the increment of $\underline{r}_{2}$ is likely to be minimal. We would not be surprised if by performing the verification in Equations $\left({\ref{veri1}}\right)$ and $\left({\ref{veri2}}\right)$ above the code would be faster than blindly incrementing $\overline{r}_{1}$ and decrementing $\underline{r}_{2}$ (we did not have the time to check this.)

Finally, by using symmetry, we can reduce the analysis of questions $\left({\ref{qd1}}\right)$ and $\left({\ref{qs1}}\right)$ in the case in which ${\mathbb{f}}\!\left(t_{k}\right)$ does not contain $0$ to the cases described in Figure 2. In this Figure, the inclined lines have equations

and we may either have no intersection or one intersection. We can use the same technique described above for find intersections which are correctly rounded and such that ${f}\!\left(u\right)$ is different from zero in the new extreme points $u$, without evaluating ${\mathbb{f}}\!\left(u\right)$ (the old extreme points stay as they were.)

In order to appreciate the content of this article, one must realize that interval arithmetic is quite different from floating point arithmetic, and it is used for different purposes. Floating point methods are usually faster than interval methods, but do not provide rigorous bounds on their results. We use interval arithmetic when we want global solutions to our problems, with rigorous bounds on them, and we are willing to pay more for that.

Interval arithmetic is fundamentally limited. When we evaluate the extension $\mathbb{f}$ at a point $t$ obtain an interval $\mathbb{w}={\mathbb{f}}\!\left(t\right)$ containing ${f}\!\left(t\right)$. The width of $\mathbb{w}$ depends on the quality of the implementation of $\mathbb{f}$ and the underlying arithmetic. When $f$ is not too complex and ${f}\!\left(t\right)$ is ${O}\!\left(1\right)$, we expect that the width of $\mathbb{w}$ to be ${O}\!\left(\epsilon\right)$, where $\epsilon$ is the machine precision. It is also reasonable to expect that the functions $\underline{\mathbb{f}},\overline{\mathbb{f}}:\mathbb{x}\rightarrow\mathds{R}{}$ will be as in Figure 3 when $r\in\mathbb{x}$ is a root of $f$.

When the evaluation of $f$ using interval arithmetic leads to wide intervals, we can (and do) use centered forms and Taylor models to reduce the width of ${\mathbb{f}}\!\left(\mathbb{x}\right)$, but this will not free us from facing the reality that the evaluation of interval extensions is usually imperfect. This is a fact, and we need to be prepared to handle not only the monotone situation described in Figure 3, but also the noisy situation described in Figure 4, which often occurs near a double root, and gets worse near roots of higher multiplicities.

Another fundamental fact about interval arithmetic is that it is a pessimistic theory: it always considers the worst scenario, and it would be inconsistent if it were otherwise. As consequence, of this fundamental fact, the functions $f_{1}$ and $f_{2}$ are indistinguishable by their common extension $\mathbb{f}$, that is, if $\mathbb{f}$ is all the information that we have, then all our conclusions about $f_{1}$ must also apply to $f_{2}$. This forces us to be very conservative, and include all possible roots for $f_{1}$ as candidates for roots for $f_{2}$ too, and vice versa. We emphasize this point because we found that it is not always clear in the minds of people trying to implement roots finders like the one we propose here, and this makes them under estimate the effort required for this task in real life.

Taking into account all that was said since the beginning of this section, we reached the conclusion that we should use three tolerances to try to handle the problems caused by the very nature of interval arithmetic. First, we should have a tolerance $\tau_{x}$ to handle “normal” roots $r$ like the one in Figure 3, at which ${f}^{\prime}\!\left(r\right)\neq 0$, so that we consider intervals around such roots with width below $\tau_{x}$ to be good enough. Clusters of zeros of $\mathbb{f}$, as in Figure 4, are a different kind of beast, as the experience with other problems in interval arithmetic has show [5, 6]. By using the same tolerance $\tau_{x}$ for them we can end up with literally thousands of small intervals as candidates for containing roots, in the favorable case in which the program actually ends and returns something. We say this from our experience with our own interval root finders as well as with interval root finders developed by other people.

Our algorithm asks the user to also provide a “cluster tolerance” $\tau_{c}$, which could be something of the order of $\sqrt{\tau_{x}}$, and a tolerance $\tau_{w}$ so that function values $w={f}\!\left(t\right)$ with magnitude below $\tau_{w}$ are deemed to be negligible. We take the liberty of increasing $\tau_{w}$ if we find it to be smaller than $16$ times the maximum width of the intervals ${\mathbb{f}}\!\left(t\right)$ which compute along the execution of algorithm. Contrary to what is proposed in [4], we believe that all tolerances should be absolute, and not relative (and users are free to choose which approach suits their needs best.) Using these tolerances, we can implement Algorithm 1 which “expands” a zero $z\in\mathbb{x}$. By applying this algorithms to the point $z$ and the interval $\mathbb{x}=[\underline{x},\overline{x}]$ in Figure 4 we would identify the cluster $[\underline{c},\overline{c}]$ and submit the intervals $[\underline{x},\underline{c}]$ and $[\overline{c},\overline{x}]$ to further examination. Finally, we must say that the actual algorithms for zero expansion in the C++ code are more involved than Algorithm 1, but the details are too cumbersome to be presented here.

As mentioned in the introduction, the main tool for root finding is Newton’s method

This method was devised to find a “point approximation” to the root $r$, that is, we for $\tilde{r}$ such that ${\left|{r-\tilde{r}}\right|}$ is small. It is our opinion that the goal of interval arithmetic is a bit different: we look for a short interval $\mathbb{r}$ such that $r\in\mathbb{r}$. Of course, in the end both goals amount to the same, but they suggest slightly different perspectives. We can rephrase the interval arithmetic goal as finding a short interval $\mathbb{r}=[\underline{r},\overline{r}]$ such that ${f}\!\left(\underline{r}\right)$ and ${f}\!\left(\overline{r}\right)$ have opposite signs (modulo degenerate cases in which ${f}\!\left(\underline{r}\right)=0$ or ${f}\!\left(\overline{r}\right)=0$.) From this perspective, we believe that we should modify the classic Newton step $\left({\ref{newton2}}\right)$ so that it produces iterates $t_{k}$ such that the signs of $w_{k}:={f}\!\left(t_{k}\right)$ alternate, in order to ensure that the interval

always contain a root.

A simple modification with this purpose is described in Figure 5. This Figure illustrates only the case in which $w_{k}<0$ and $d_{k}>0$, but cases in which $w_{k}>0$ or $d_{k}<0$ are similar and can be reduced to the case $w_{k}<0$ and $d_{k}>0$ by replacing $f$ by $-f$ or $t$ by $-t$. A simple way to obtain a reasonable approximation $h_{k}$ for the second derivative mentioned in Figure 5 is to use the quotient

with data $d_{k-m}$ and $t_{k-m}$ from previous iterations.

Motivated by Figure 5, we propose Algorithm 2 below for finding a short interval $\mathbb{r}$ containing the single root of $f$ in the interval $\mathbb{x}$ using exact arithmetic, in the case in which ${f}\!\left(\underline{x}\right)<0<{f}\!\left(\overline{x}\right)$ and ${f}^{\prime}\!\left(t\right)>0$ for $t\in\mathbb{x}$. It is not difficult to prove that Algorithm 2 has the same properties as the many other versions of Newton’s method that we find in the literature:

• •

  When the model

  $${f}\!\left(t\right)\approx{f}\!\left(t_{k}\right)+{f}^{\prime}\!\left(t_{k}\right)\left({t-t_{k}}\right)+{f}^{\prime\prime}\!\left(t\right)\left({t-t_{k}}\right)^{2}/2$$

  (19)

  is inaccurate for $t\in\mathbb{x}$, we fallback into a bisection step. This guarantees that eventually the interval $\mathbb{x}$ will become very short, and the analysis in the next item will apply

• •

  When the quadratic model $\left({\ref{quadratic}}\right)$ is accurate in $\mathbb{x}$, and the interval $\mathbb{x}$ is very short, $h$ will be a good approximation of ${f}^{\prime\prime}\!\left(t_{k}\right)$ and our modification in the Newton step will ensure that bisection will stop being necessary, the signs of the $w_{k}$ will alternate and the iterates will converge at the same quadratic rate as the classic method.

This Section presents an interval Newton’s method for functions with ${f}^{\prime}\!\left(t\right)>\kappa>0$ for $t\in\mathbb{x}$ (in order to used for functions with ${f}^{\prime}\!\left(t\right)<-\kappa<0$ we can simply replace $f$ by $-f$. The method is a simple adaptation of the “point method” presented in Section 4, taking into account the fundamental differences between exact (point) arithmetic and interval arithmetic. Among such differences we have the following ones:

• •

  In interval arithmetic, for a strictly increasing function, we may have that $0\in{\mathbb{f}}\!\left(t\right)$ for many $t$’s. In exact arithmetic this can only happen for one value of $t$.

• •

  In interval arithmetic we have a natural measure for the error in the the values of $f$, given by the maximum width of the intervals ${\mathbb{f}}\!\left(t_{k}\right)$. This allows us to correct a poor choice of the tolerance $\tau_{w}$ by the user and dispense with the cluster tolerance $\tau_{c}$ mentioned in Section 3.

As we mentioned in the introduction, the increasing case is much simpler than the general one. In this case, we do not need to compute neither ${\mathbb{f}}\!\left(\mathbb{x}_{k}\right)$ nor ${\mathbb{f}}^{\prime}\!\left(\mathbb{x}_{k}\right)$. We only need the interval ${\mathbb{f}}\!\left(t_{k}\right)$ and the floating point number ${f}^{\prime}\!\left(t_{k}\right)$, and this entails a considerable reduction in the cost of the steps. The analysis of degenerate cases is also much simpler in the increasing case. These ideas are explored in Algorithm 3, which is presented after the bibliography. For brevity, we omit the expand_zero label in this code, but it is similar to the zero expansion algorithms presented in Section 3.

In practice, it is quite difficult to implement the algorithm presented here due to the need of much attention to detail. It is quite likely that the first attempts at such an implementation will fail, as ours did. Even worse, it is also likely that the failure will be undetected, specially if the details involved in the actual coding of the algorithm are considered to be of little relevance. As a result, in order to have some assurance of the reliability of our code in practice it is necessary to have a good suit of tests. Readers should not underestimate this point (and should not underestimate our remark that they should not underestimate this point, and, recursively…)

It is also important to realize that a good test suit may reveal problems not only in the code but also in the theory upon which it is based. A good example of this point is the choice of the stopping criterion presented in chapter 9 of [4]. When using that stopping criterion for the simple polynomial