Principle

Let $x_1,\ldots,x_m$ be the set of unknowns and let ${\cal I}_1=\{(\underline{x_1^1},\overline{x_1^1}), \ldots (\underline{x_m^1},\overline{x_m^1})\}$ be the set of $m$ intervals in which you are searching the solutions of the $n$ equations ${\cal F}_1(x_1,\ldots,x_m)=0, \ldots {\cal F}_n(x_1,\ldots,x_m)=0$ (for the sake of simplicity we don't consider inequalities but the extension to inequalities is straightforward).

We will denote by $F_i$ the interval value of ${\cal F}_i$ when this function is evaluated for the box $(\underline{x_1},\overline{x_1}),\ldots,(\underline{x_m},\overline{x_m})$ of the unknowns while $F({\cal I}_j)$ will denote the $n$-dimensional interval vector constituted of the $F_i$ when the unknowns have the interval value defined by the set ${\cal I}_j$.

The algorithm will use a list of boxes ${\cal I}$ whose maximal size $M$ is an input of the program. This list is initialized with ${\cal I}_1$. The number of ${\cal I}$ currently in the list is $N$ and therefore at the start of the program $N=1$. The algorithm will also use an accuracy on the variable $\epsilon$ and on the functions $\epsilon_F$. The norm of a ${\cal I}$ is defined as:

$$\Vert{{\cal I}^j}\Vert={\rm Max}(\overline{x_k^j}-\underline{x_k^j})~~~{\rm for}~k \in [1,m]$$

The norm of the interval vector $F({\cal I}_j)$ is defined as:

$$\uparrow{F({\cal I}_j)}\uparrow ={\rm Max}(\overline{F_k({\cal I}_j)}- \underline{F_k({\cal I}_j)})~~~{\rm for}~k \in [1,n]$$

The algorithm uses an index $i$ and the result is a set ${\cal S}$ of interval vector ${\cal S}_k$ for the unknowns whose size is $S$. We assume that there is no $F_k$ with $k$ in $[1,n]$ such that $\overline{F({\cal I}_1)}<0$ or $\underline{F({\cal I}_1)}>0$, otherwise the equations have no solution in ${\cal I}_1$. Two lists of interval vectors ${\cal V}, {\cal W}$ whose size is $2^m$ will also be used. The algorithm is initialized with $i = 1, S=1$ and proceed along the following steps:

1. if $i = N+1$ return $S-1$ and ${\cal S}$ and exit
2. bisect ${\cal I}_i$ which produce $2^m$ new interval vectors ${\cal V}_l$ and set $j=1$
3. for $l=1,\ldots,2^m$
   1. evaluate $F({\cal V}_l)$
   2. if it exist $F_k$ with $k$ in $[1,n]$ such that $\overline{F_k({\cal V}_l)}<0$ or $\underline{F_k({\cal V}_l)}>0$, then $l = l+1$ and go to step 3
   3. if $\Vert{{\cal V}_l}\Vert < \epsilon$ or $\uparrow{F({\cal V}_l)}\uparrow <\epsilon_f$, then store ${\cal V}_l$ in ${\cal S}_S$, increment $S$ and go to step 3
   4. store ${\cal V}_l$ in ${\cal W}_j$, increment $j$ and go to step 3
4. if $j=1$ increment $i$ and go to step 1
5. if $N+j-2 >M$ return a failure code as there is no space available to store the new intervals
6. if $j>1$ store one of the ${\cal W}$ in ${\cal I}_i$, the other $j-2$ at the end of ${\cal I}$, starting at position $N+1$. Add $j-2$ to $N$ and go to step 1

Basically the algorithm just bisect the box until either their width is lower than $\epsilon$ or the width of the interval function is lower than $\epsilon_f$ (provided that there is enough space in the list to store the intervals). Then if all the intervals functions contain 0 we get a new solution, if one of them does not contain 0 there is no solution of the equations within the current box. A special case occurs when all the components of the box are reduced to a point, in which case a solution is obtained if the absolute value of the interval evaluation of the function is lower than $\epsilon_f$.

Now three problems have to be dealt with:

1. how to choose the ${\cal W}$ which will be put in place of the ${\cal I}_i$ and in which order to store the other ${\cal W}$ at the end of the list?
2. can we improve the management of the bisection process in order to conclude the algorithm with a limited number $M$?
3. how do we distinguish distinct solutions ?

The first two problems will be addressed in the next section.