The multivariate bisection algorithm

TL;DR

This paper introduces a multivariate bisection algorithm based on the Poincaré-Miranda theorem for solving nonlinear systems in ^n, demonstrating local convergence and providing a numerical implementation in two dimensions.

Contribution

It proposes a novel multivariate bisection method supported by the Poincaré-Miranda theorem, with proven local convergence and practical numerical implementation.

Findings

Algorithm converges locally for systems satisfying Poincare9-Miranda conditions.

Numerical implementation successfully applied in two-dimensional cases.

Provides a new root-finding approach for nonlinear systems in higher dimensions.

Abstract

The aim of this paper is the study of the bisection method in $\mathbb{R}^n$. In this work we propose a multivariate bisection method supported by the Poincar\'e-Miranda theorem in order to solve non-linear system of equations. Given an initial cube verifying the hypothesis of Poincar\'e-Miranda theorem the algorithm performs congruent refinements throughout its center by generating a root approximation. Throughout preconditioning we will prove the local convergence of this new root finder methodology and moreover we will perform a numerical implementation for the two dimensional case.

Tables

Table 1. Table 2: Root’s approximation, evaluation, iteration performed and initial guess for testing maps.

F	$c$	$\text{F}\left(c\right)$	iter	$K_0$
${\text{F}}_1$	0.618033988749895	0.004965068306495 1e-14	51	$[0,1]\times [0,1]$
	0.786151377757422	-0.123942463016433 1e-14
$\text{F}2$	0.567143290409784	0.111022302462516 1e-15	50	$[0,1]\times [0,1]$
	0.567143290409784	0.111022302462516 1e-15
${\text{F}}_3$	0.378316940137480	0.139577647543639 1e-15	51	$[0,1]\times [0,1]$
	0.507403383528753	-0.072495394968176 1e-15
${\text{F}}_4$	0.926174872358938	0.129347223584252 1e-15	49	$[0,1]\times [-1,0]$
	-0.582851662173280	-0.115653908517277 1e-15
${\text{F}}_5$	0.353246619596717	-0.244439451327881 1e-15	52	$[0,1.1]\times [0,2]$
	0.606081736641465	0.047257391058546 1e-15
${\text{F}}_6$	0.510030862987151	-0.045236309398304 1e-13	42	$[0.4,1]\times [0,0.4]$
	0.048996913701194	-0.904901681894059 1e-13

Equations

$$K=\{x\in {\mathbb{R}}^n:|x_j-\widehat{x_j}|\le \rho ,j=1(1)n\}$$

$$F_i^{+}=\{x\in K:x_i=\widehat{x_i}+\rho \},F_i^{-}=\{x\in K:x_i=\widehat{x_i}-\rho \}$$

$$f_i(x)f_i(y)\le 0,x\in F_i^{+},y\in F_i^{-}$$

$$\Vert M{\Vert }_p=\underset{\Vert x{\Vert }_p=1}{\max }\Vert Mx{\Vert }_p\text{mbox}wherep=\infty ,2\text{mbox}or1$$

$$\text{mbox}\mathbf{G}(X)\approx \text{mbox}\mathbf{G}(X_0)+D\text{mbox}\mathbf{G}(X_0)(X-X_0)=\text{mbox}\mathbf{G}(X_0)+X-X_0$$

$$\text{mbox}\mathbf{G}(X)\approx X-X_0$$

$$c_{1}=\bigg{(}\frac{a_{1_{0}}+b_{1_{0}}}{2},\frac{a_{2_{0}}+b_{2_{0}}}{2},...,\frac{a_{n_{0}}+b_{n_{0}}}{2}\bigg{)}$$

$$c_{2}=\bigg{(}\frac{a_{1_{1}}+b_{1_{1}}}{2},\frac{a_{2_{1}}+b_{2_{1}}}{2},...,\frac{a_{n_{1}}+b_{n_{1}}}{2}\bigg{)}$$

$$\text{mbox}{\mathbf{G}}_1(X):=D\text{mbox}\mathbf{F}(c_1{)}^{-1}\text{mbox}\mathbf{F}(X)$$

$$\mbox{\bf{G}}_{k}(X):=\left\{\begin{array}[]{lcl}\mbox{\bf{G}}_{k-1}(X)\ \mbox{ if}\ \mbox{\bf{G}}_{k-1}\ \mbox{verifies P.M. in}\ Q_{k}\\ &&\\ D\mbox{\bf{G}}_{k-1}(c_{k})^{-1}\mbox{\bf{G}}_{k-1}(X)\ \mbox{if}\ \mbox{\bf{G}}_{k-1}\ \mbox{does not verify P.M. in}\ Q_{k}\end{array}\right.$$

$$K_{k+1}\subset K_k=[a_{1_k},b_{1_k}]\times ...\times [a_{n_k},b_{n_k}]\text{mbox}wheretheverticesverify$$

$$a_{j_0}\le a_{j_1}\le a_{j_2}\le \dots \le a_{j_k}\le \dots \le b_{j_0}$$

$$b_{j_0}\ge b_{j_1}\ge b_{j_2}\ge \dots \ge b_{j_k}\ge \dots \ge a_{j_0}$$

$$a_{j_k}-b_{j_k}=\frac{a_{j_{k-1}}-b_{j_{k-1}}}{2}=\dots =\frac{a_{j_0}-b_{j_0}}{2^k}$$

$$c_{k}=\bigg{(}\frac{a_{1_{k}}+b_{1_{k}}}{2},\frac{a_{2_{k}}+b_{2_{k}}}{2},...,\frac{a_{n_{k}}+b_{n_{k}}}{2}\bigg{)}$$

$$\text{mbox}{\mathbf{G}}_k(X)=D\text{mbox}{\mathbf{G}}_{k-1}(c_k{)}^{-1}\text{mbox}{\mathbf{G}}_{k-1}(X)$$

$$=D\text{mbox}{\mathbf{G}}_{k-1}(c_k{)}^{-1}D\text{mbox}\mathbf{F}(c_{k-1}{)}^{-1}\text{mbox}\mathbf{F}(X)=D\text{mbox}\mathbf{F}(c_k{)}^{-1}D\text{mbox}\mathbf{F}(c_{k-1})D\text{mbox}\mathbf{F}(c_{k-1}{)}^{-1}\text{mbox}\mathbf{F}(X)$$

$$=D\text{mbox}\mathbf{F}(c_k{)}^{-1}\text{mbox}\mathbf{F}(X).$$

$$\underset{k\to \infty }{\lim }{a}_{j_k}=\underset{k\to \infty }{\lim }{b}_{j_k}=r_j$$

$$g_{1_k}(a_{1_k},x_2,..,x_n)g_{1_k}(b_{1_k},x_2,..,x_n)\le 0\forall x\in {\widehat{K_k}}^1$$

$$........$$

$$g_{j_k}(x_1,..,x_{j-1},a_{j_k},x_{j+1},..,x_n)g_{j_k}(x_1,..,x_{j-1},b_{j_k},x_{j+1},..,x_n)\le 0\forall x\in {\widehat{K_k}}^j$$

$$........$$

$$g_{n_k}(x_1,..,x_{n-1},a_{n_k})g_{n_k}(x_1,..,x_{n-1},b_{n_k})\le 0\forall x\in {\widehat{K_k}}^n$$

$$\{p\}=\bigcap _{n=0}^{\infty }{K}_k$$

$$g_{1_k}(a_{1_k},r_2,..,r_n)g_{1_k}(b_{1_k},r_2,..,r_n)\le 0\forall k\in \mathbb{N}$$

$$........$$

$$g_{n_k}(r_1,..,r_{n-1},a_{n_k})g_{n_k}(r_1,..,r_{n-1},b_{n_k})\le 0\forall k\in \mathbb{N}$$

$$c_{k}=\bigg{(}\frac{a_{1_{k}}+b_{1_{k}}}{2},\frac{a_{2_{k}}+b_{2_{k}}}{2},...,\frac{a_{n_{k}}+b_{n_{k}}}{2}\bigg{)}\underset{k\rightarrow\infty}{\longrightarrow}(r_{1},..,r_{n})=r$$

$$D\text{mbox}\mathbf{F}(c_k{)}^{-1}\to D\text{mbox}\mathbf{F}(r{)}^{-1}$$

$$\Vert D\text{mbox}\mathbf{F}(c_k{)}^{-1}F(X)-D\text{mbox}\mathbf{F}(r{)}^{-1}F(X){\Vert }_{\infty }\le \Vert D\text{mbox}\mathbf{F}(c_k{)}^{-1}-D\text{mbox}\mathbf{F}(r{)}^{-1}{\Vert }_{\infty }\Vert F(X){\Vert }_{\infty }\to 0$$

$$g_{j_k}(X)\underset{k\to \infty }{\to }{g}_j(X).$$

$$g_{j_k}(r_1,..,r_{j-1},a_{j_k},r_{j+1},..,r_n)\to g_j(r_1,..,r_j,..,r_n)\text{mbox}foreachj=1(1)n$$

$$g_{j_k}(r_1,..,r_{j-1},b_{j_k},r_{j+1},..,r_n)\to g_j(r_1,..,r_j,..,r_n)\text{mbox}foreachj=1(1)n,$$

Full text

1 Introduction

The problem of finding numerical approximations to the roots of a non-linear system of equations was subject of various studies, different methodologies have been proposed between optimization and Newton’s procedures. In [2] D. H. Lehmer proposed a method for solving polynomial equations in the complex plane testing increasingly smaller disks for the presence or absence of roots. In other work, Herbert S. Wilf developed a global root finder of polynomials of one complex variable inside any rectangular region using Sturm sequences[11].

The classical Bolzano’s theorem or Intermediate Value theorem ensure that a continuous function that changes sign in an interval has a root, that is, if $f:[a,b]\rightarrow\mathbb{R}$ is continuous and $f(a)f(b)<0$ then there exist $c\in(a,b)$ that $f(c)=0$. In the multidimensional case the generalization of this result is the known Poincaré-Miranda theorem that ensures that if we have $f_{1},...,f_{n}$ $n$-continuous functions of $n$ variables $x_{1},...,x_{n}$ and the variables are subjected to vary between $a_{i}$ and $-a_{i}$ then if $f_{i}(x_{1},..,a_{i},..,x_{n})f_{i}(x_{1},..,-a_{i},..,x_{n})<0$ for all $x_{i}$ then there exist $c\in[-a_{i},a_{i}]^{n}$ such that $f(c)=0$. This result was announced the first time by Poincaré in 1883 [7] and published in 1884 [8] with reference to a proof using homotopy invariance of the index. The result obtained by Poincaré has come to be known as the theorem of Miranda, who in 1940 showed that it is equivalent to the Brouwer fixed point [6]. For different proofs of the Poincaré-Miranda theorem in the $n$-dimensional case, see [1], [10].

Theorem 1.1**.**

(Poincaré-Miranda theorem). Let $K$ be the cube

[TABLE]

where $\rho\geq 0$ and $\mbox{\bf{F}}=(f_{1},f_{2},..,f_{n}):K\rightarrow\mathbb{R}^{n}$ a continuous map on $K$. Also let,

[TABLE]

be the pairs of parallel opposite faces of the cube $K$.

If for $i=1(1)n$ the $i$-th component $f_{i}$ of F has opposite sign or vanishes on the corresponding opposite faces $F_{i}^{+}$ and $F_{i}^{-}$ of the cube $K$, i.e.

[TABLE]

then the mapping F has at least one zero point $r=(r_{1},r_{2},..,r_{n})$ in K.

Throughout this paper we will recall the opposite signs condition 1.1 as the Poincaré-Miranda property P.M.. The aim of this work is to develop a bisection method that allows us solve non-linear system of equations $\mbox{\bf{F}}(X)=0,\ X=(x_{1},x_{2},..,x_{n})$ using the above Poincaré-Miranda theorem. The idea of the algorithm will be similar as in the classical one dimensional algorithm, we perform refinements of the cube domain in order to check the sign conditions on the parallel faces. In one dimension it is clear that an initial sign change in the border of an interval produces other sign change in a half partition of it but in several dimension we cannot guarantee that the Poincaré-Miranda conditions maintain after a refinement. Even if $r$ is an exact solution, there may not be any such $K$ (for which 1.1 holds). However, J. B. Kioustelidis [3] has pointed out that, for $\hat{x}$ close to a simple solution (where the Jacobian is nonsingular) of $\mbox{\bf{F}}(X)$, Miranda’s theorem will be applicable to some equivalent system for suitable $K$. Therefore in case of a fail in the sign conditions with the original system, we should try to transform it. The idea will be find an equivalent system throughout non-linear preconditioning where the equations are better balanced in the sense that the new system could be close to some hyperplane in order to improve the chances to check the sign conditions in some member of the refinement.

We will denote the infinite norm by $\|x\|_{\infty}=\max\{|x_{1}|,..,|x_{n}|\}$, the Euclidean norm by $\|x\|_{2}=\sqrt{x_{1}^{2}+..+x_{n}^{2}}$ and the $1$-norm by $\|x\|_{1}=|x_{1}|+..+|x_{n}|$. Given a vector norm on $\mathbb{R}^{n}$, the associated matrix norm for a matrix $M\in\mathbb{R}^{n\times n}$ is defined by

[TABLE]

It is know that in the case of the $\infty$-matrix norm it can be expressed as a maximum sum of its row, that is if $M=(m_{ij})$ then $\|M\|_{\infty}=\max_{i=1,..,n}\sum_{j=1}^{n}|m_{ij}|$, therefore it is easy to see that a sequence of matrices $(M_{k})_{k}$ converge if and only if their coordinates converge. Since the domains involved are multidimensional cubes, the most proper norm to handle the distance will be the $\infty$-norm.

We will accept $r$ as a root with a small tolerance level $\delta$ if $\|\mbox{\bf{F}}(r)\|_{p}\leq\delta$.

2 The algorithm and its description

This section gives a step-by-step description of our algorithm, the core of it meets in the classical bisection algorithm in one dimension.

Definition 2.1**.**

A $2^{n}$-refinement of a cube $K\subset\mathbb{R}^{n}$ is a refinement into $2^{n}$ congruent cubes $Q=\{K^{1},K^{2},..,K^{2^{n}}\}$.

We say that a $2^{n}$-refinement of $Q$ verifies the Poincaré-Miranda condition if there exist $K^{l}\in Q$ such that $\mbox{\bf{F}}:K^{l}\rightarrow\mathbb{R}^{n}$ verifies the condition of Theorem 1.1.

Given a system $\mbox{\bf{F}}(X)=0$, the preconditioned system is $\mbox{\bf{G}}(X)=M\mbox{\bf{F}}(X)=0$ for some matrix $M$ such that the jacobian at $X_{0}$ verifies $D\mbox{\bf{G}}(X_{0})=\mbox{Id}$. Since $D\mbox{\bf{G}}(X_{0})=MD\mbox{\bf{F}}(X_{0})$ it turns out that $M=D\mbox{\bf{F}}(X_{0})^{-1}$ and it is clear that the preconditioned system is an equivalent system of F and both have the same roots. After preconditioning, the equations in $\mbox{\bf{G}}(X)=0$ are close to a hyperplane having equation $x=k_{1}$ and $y=k_{2}$, where $k_{i}$ is some constant. This fact comes from the Taylor expansion of G around $X_{0}$, indeed if $X$ is close to $X_{0}$ then

[TABLE]

and therefore it is clear that the equations are close to some hyperplane. Moreover if $X_{0}$ is nearly a zero point of F then,

[TABLE]

and therefore it will behave like the components of $X-X_{0}$, and take nearly opposite values on the corresponding opposite faces of the cube.

2.1 Algorithm procedure

The multivariate bisection algorithm proceeds as follows:

We start choosing an initial guess $K_{0}=[a_{1_{0}},b_{1_{0}}]\times[a_{2_{0}},b_{2_{0}}]\times...\times[a_{n_{0}},b_{n_{0}}]\subset\mathbb{R}^{n}$ verifying the Poincaré-Miranda condition on F. 2. 2.

We locate the center

[TABLE]

of $K_{0}$. 3. 3.

Generate a first $2^{n}$-refinement $Q_{1}$ through $c_{1}$. 4. 4.

If $Q_{1}$ verifies the Poincaré-Miranda condition, let $K_{1}=[a_{1_{1}},b_{1_{1}}]\times[a_{2_{1}},b_{2_{1}}]\times...\times[a_{n_{1}},b_{n_{1}}]$ be the quarter of $Q_{1}$ where the conditions of Theorem 1.1 are verified, we chose

[TABLE]

the center of $K_{1}$. If $Q_{1}$ does not verify the Poincaré-Miranda condition we preconditioning the system in $c_{1}$ setting

[TABLE]

and then we check again the sign conditions with the preconditioned system $\mbox{\bf{G}}_{1}(X)$ in $Q_{1}$.

This recursion is repeated while the Poincaré-Miranda condition are verified, generating a sequence of equivalent system

[TABLE]

and a decreasing cube sequence $K_{k}$, such that

[TABLE]

for each $j=1(1)n$ and where the length of the current interval $[a_{j_{k}},b_{j_{k}}]$ is a half of the last iteration,

[TABLE]

The root’s approximation after $k$-th iteration will be,

[TABLE]

and the method is stopped until the zero’s estimates gives sufficiently accuracy or until the Poincaré-Miranda condition leaves to maintain.

Remark 2.2**.**

It is easy to see that the $k$th-preconditioning system $\mbox{\bf{G}}_{k}(X)=(g_{1_{k}}(X),..,g_{n_{k}}(X))$ can be expressed as $D\mbox{\bf{F}}(c_{k})^{-1}\mbox{\bf{F}}(X)$, indeed, by induction suppose that it is true for $k-1$, then differencing and valuing in $c_{k}$ we have,

[TABLE]

[TABLE]

[TABLE]

Since we cannot always ensure that a refinement of a given cube will verify the Poincaré-Miranda condition, we cannot ensure the converge for any map that only has a sign change in a given initial cube. So, in case of a fail in the sign conditions in some step, we try to rebalance the system using preconditioning in the center of the current box recursion. The preconditioning allows us to increase the chances to be more often in the sign conditions and therefore keep going with the quadrisection procedure in order to get a better root’s approximation. In [3], J. B. Kioustelidis found sufficient conditions for the validity of the Poincaré-Miranda Miranda condition for preconditioning system, there it was proved that the sign conditions are always valid if the center of the cube $K$ is close enough to some root of F. So, if we start the multivariate bisection algorithm with an initial guess close to some root, Kioustelidis’s theorem will guarantee the validity of Poincaré-Miranda in each step of our method allowing the local convergence of it.

In the next theorem we will prove the local convergence for the multivariate bisection algorithm when we preconditioning in each step.

Theorem 2.3**.**

Let $\mbox{\bf{F}}=(f_{1},...,f_{n}):K_{0}\rightarrow\mathbb{R}^{n}$ be a $C^{2}$ map defined on the cube $K_{0}=\{x\in\mathbb{R}^{n}:\|x-c_{1}\|_{\infty}\leq\rho\}=[a_{1_{0}},b_{1_{0}}]\times..\times[a_{n_{0}},b_{n_{0}}]$ with $\rho$ small enough verifying the Poincaré-Miranda sign condition; assume that $D\mbox{\bf{F}}(X)$ is invertible for all $X\in K_{0}$, furthermore suppose that we perform the preconditioning in each step then the multivariate bisection algorithm generates a sequence $c_{k}$ such that

Starting at $K_{0}$, $c_{k}\stackrel{{\scriptstyle\|.\|}}{{\longrightarrow}}r$ with $\mbox{\bf{F}}(r)=0$. 2. 2.

$\|c_{k}-r\|_{2}\leq\dfrac{\sum^{n}_{j=1}b_{j_{0}}-a_{j_{0}}}{2^{k}}$.

Proof.

The Poincaré-Miranda sign conditions guarantee the existence of a root inside $K_{0}$ and given a refinement $Q_{1}$ of $K_{0}$ since $\rho$ is small enough Item c of Theorem 2 in [3] guarantees the validity of Poincaré-Miranda sign conditions for a member of $Q_{1}$. Performing successive refinements we will always find a member $K_{k}$ of the refinement $Q_{k}$ verifying the sign conditions for the preconditioned system $\mbox{\bf{G}}_{k}(X)$. For each $j=1(1)n$ the sequences $(a_{j_{k}})_{k},(b_{j_{k}})_{k}$ are monotones and bounded and therefore they converge. From equation 2.4 we have for each $j=1(1)n$,

[TABLE]

and from the border conditions,

[TABLE]

where $\widehat{K_{k}}^{j}$ means that the $j$ coordinate of $K_{k}$ is omitted. Since the diameter of $K_{k}$ tends to zero by Cantor’s intersection theorem the intersection of the $K_{k}$ contains exactly one point,

[TABLE]

and the equations 2.5 guarantee that $p=(r_{1},..,r_{n})$. Then, we can evaluate equations 2.6 in $p=(r_{1},..,r_{n})$ getting,

[TABLE]

It is clear that,

[TABLE]

then by the continuity of $D\mbox{\bf{F}}$ and the continuity of the inversion in the $\infty$-matrix norm we have

[TABLE]

Let $G(X)=D\mbox{\bf{F}}(r)^{-1}F(X)=(g_{1}(X),..,g_{n}(X))$, since

[TABLE]

for each $X\in K_{0}$ then we get the punctual convergence for each coordinate function

[TABLE]

From equations 2.5 we have,

[TABLE]

[TABLE]

therefore taking limit in equations 2.7 we get

[TABLE]

and finally it is clear that $\mbox{\bf{F}}(r_{1},..,r_{n})=0$. 2. 2.

Let $(c_{j_{k}})_{k}\ j=1(1)n$, be the coordinates of the sequence $(c_{k})_{k}$ we have following estimation

[TABLE]

Indeed, since the sequences $(a_{j_{k}})$ and $(b_{j_{k}})$ are monotones and bounded by $r_{j}$ we get for each $j=1(1)n$,

[TABLE]

On the other hand,

[TABLE]

Therefore,

[TABLE]

∎

As in the classical one dimensional bisection algorithm Item 2 of theorem 2.3 gives a way to determine the number of iterations that the bisection method would need to converge to a root to within a certain tolerance. The number of iterations needed, $k$, to achieve the given tolerance $\delta$ is given by,

[TABLE]

The following example shows an infinite application of the bisection algorithm in $\mathbb{R}^{2}$ with non-preconditioning.

Example 2.4**.**

Consider the map $\mbox{\bf{F}}(x,y)=(y+x-1,y-e^{-x^{2}})$, we start checking the Poincaré-Miranda condition on $K_{0}=[0,1]\times[0,1]$,

[TABLE]

then if we considerate the $K^{3}$ quarter of each quadrisection, we can check that it always will verify the Poincaré-Miranda condition and therefore it is not necessary preconditioning in each step. Let $a_{1_{k}}=0$, $b_{1_{k}}=\frac{1}{2^{k}}$, $a_{2_{k}}=1-\frac{1}{2^{k}}$ and $b_{2_{k}}=1$ be the coordinates of the $k$-th $K^{3}$ quarter rectangle, we have to prove that

[TABLE]

[TABLE]

Indeed, the first inequality is clear and follows directly from the domain of $y$, the second follows from the fact that the domain of $x$ implies that

[TABLE]

getting

[TABLE]

The $k$-th root’s approximation is,

[TABLE]

and the error verifies

[TABLE]

3 Implementation, performance and testing

Throughout this section we will focus in the implementation and performance of the bisection algorithm in $\mathbb{R}^{2}$. The bisection algorithm was developed in Matlab in a set of functions running from a main function. In order to check the P.M. conditions for the function $\mbox{F}=(f_{1},f_{2})$ we need to compute the intervals $f_{i}(F_{i}^{+}),f_{i}(F_{i}^{-})$ ($i=1,2$) and one way to achieve this is by using Interval Analysis (IA). IA was marked by the appearance of the book Interval Analysis by Ramon E. Moore in 1966 [4] and it gives a fast way to find an enclosure for the range of the functions. A disadvantage of IA is the well known overestimation. If intervals $f_{i}(F_{i}^{+}),f_{i}(F_{i}^{-})$ are available then the P.M. follows from the condition

[TABLE]

Interval-Valued Extensions of Real Functions gives a way to find an enclosure of the range of a given real-valued function. Most generally, if we note by $[\mathbb{R}]$ the set of all finite intervals, we say that $[f]:[\mathbb{R}]^{n}\rightarrow[\mathbb{R}]$ is an interval extension of $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ if

[TABLE]

where $X=(X_{1},...,X_{n})$ represents a vector of intervals. There are different kind of interval functional extensions; if we have the formula of a real-valued function $f$ then the natural interval extension is achieved by replacing the real variable $x$ with an interval variable $X$ and the real arithmetic operations with the corresponding interval operations. Another useful interval extension is the mean value form. Let $m=m(X)$ be the center of the interval vector $X$ and let $[\frac{\partial f_{i}}{\partial x_{i}}]$ be an interval extension of $\frac{\partial f_{i}}{\partial x_{i}}$ by the mean value theorem we have

[TABLE]

$[f_{mv}](X)$ is the mean value extension of $f$.

Let $[f_{i}](F_{i}^{+}),[f_{i}](F_{i}^{-})$ be an interval extension of $f_{i}(F_{i}^{+}),f_{i}(F_{i}^{-})$, then it is clear that if

[TABLE]

equations 3.9 and 3.10 are also true. So, in order to check the P.M. conditions along the edges we will compute equations 3.11 and 3.12.

Various interval-based software packages for Matlab are available, we have chosen the well known INTLAB toolbox [9]. The toolbox has several interval class constructor for intervals, affine arithmetic, gradients, hessians, slopes and more. Ordinary interval arithmetic has sometimes problems with dependencies and wrapping effect given large enclosures of the range and therefore overestimating the sign behaviour. A way to fight with this is affine arithmetic. In affine arithmetic an interval is stored as a midpoint $X_{0}$ together with error terms $E_{1},...,E_{k}$ and it represents

[TABLE]

where $U_{1},...,U_{k}$ are parameters independently varying within $[-1,1]$. In case of get wrong signs for $f_{i}(F_{i}^{-})$ and $f_{i}(F_{i}^{+})$ and dismiss the possibility of a not very sharp estimation of IA we also compute the interval extension but now using the affine arithmetic.

Other way to improve the enclosure of the range and get sharper lower and upper bounds is throughout subdivision or refinements. In this methodology we perform subdivision of the domain and then we take the union of interval extensions over the elements of the subdivision; this procedure is called a refinement of $[f]$ over $X$. Let $N$ be a positive integer we define

[TABLE]

We have $X_{i}=\cup^{N}_{j=1}X_{i,j}$ and $w(X_{i,j})=\dfrac{w(X_{i})}{N}$ and furthermore,

[TABLE]

The interval quantity

[TABLE]

is the refinement of $[f]$ over $X$.

The algorithms that we have performed to compute equations 3.11 and 3.12 combines all the above methodologies and were adapted from [5]. In the following steps we summarize the routines that we have performed. The mean value extension was implemented using an approximation of $[\frac{\partial f_{i}}{\partial x_{i}}]$ throughout the central finite difference of the natural interval extension of $f$, that is

[TABLE]

Algorithm 1 shows the routine for the mean value extension.

The refinement procedure was implemented twice, one for the case of mean value extension and other for the affine arithmetic implementation. Algorithm 2 computes the mean value extension over an uniform refinement of the interval $X$ with $N$ subintervals and Algorithm 3 computes the natural extension using affine arithmetic.

Now we are ready to compute equations 3.11 and 3.12 using the above algorithms. Let $K^{l}=[a_{1_{l}},b_{1_{l}}]\times[a_{2_{l}},b_{2_{l}}]=I_{1_{l}}\times I_{2_{l}}$ be a member of the refinement $Q$ and let

[TABLE]

[TABLE]

be the coordinate functions on the edges of $K^{l}$ ($l=1...4$), Algorithm 4 summarizes the routine that we have performed using IA in order to compute the sign along the edges.

Algorithm 5 summarizes the implementation of the Bisection Algorithm that we have performed in Matlab.

In order to check the accuracy and performance of the algorithm, we test it throughout different systems of equations. We start testing the algorithm in the system given in Example 2.4. We took as our starting guess the rectangle $K_{0}=[0,1]\times[0,1]$, in Table 1 we show the behaviour of the sequence throughout different tolerance levels and in Figure 3 we illustrate the procedure for tolerance level $10^{-15}$. We have chosen to use 10 digits in the mantissa representation for the root’s approximation and its evaluation and 5 digits for the norm evaluation notation. Since the system always verifies the P.M. the algorithm never performs preconditioning.

In the following steps we test the algorithm in more difficult problems, we will see that in some systems the algorithm needs to preconditioning in order to guarantee the P.M. conditions throughout the refinement. Let,

[TABLE]

be the testing maps. In Table 2 we show the numerical performance for the testing maps and in Figure 4 we illustrate the algorithm behaviour with the refinement procedure. The systems of equations and their successive possibles preconditioning are represented by a zero contour level on an mesh on the initial guess $K_{0}$ and the refinement procedure was illustrated using the rectangle Matlab’s functions. The method was implemented setting the tolerance level in $\delta=10^{-15}$ and the interval analysis refinement in N$=3$.

4 Conclusion

In this work we have clarified how a multidimensional bisection algorithm should be performed extending the idea of the classic one dimensional bisection algorithm. Due by the preconditioning in each step we could prove a local convergence theorem and we also found an error estimation. Interval Analysis allowed a fast and reliable way of computing the Poincaré-Miranda conditions and the numerical implementation showed that the method has a very good accuracy similar with the classic methods like Newton or continuous optimization.