Note that a generic analyzer based on the `ALIAS` parser has been
developed and will be presented in the chapter devoted to the
parser. This generic analyzer enable to analyze almost any type of
system in which at least one equation is algebraic in at least one of
the unknowns.

Moore theorem

and let be the norm of the matrix . If

then there is a unique solution [16] of in . This unique solution can be found using Krawczyk solving method (see section 2.10).

The previous test is implemented for and . The procedure is implemented as:

int Krawczyk_Analyzer(int m,int n, INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* J)(int, int, INTERVAL_VECTOR &),INTERVAL_VECTOR &Input)with

`m`: the number of unknown`n`: the number of equations`IntervalFunction`: a function which return the interval vector evaluation of the equations, see note 2.3.4.3`J`: a function which calculate an interval evaluation of the elements of the jacobian of the equations, see note 2.4.2.2`Input`: the ranges for the variables

Kantorovitch theorem

each being at least . Let be a point and , the norm being . Assume that is such that:

- the Jacobian matrix of the system has an inverse at such that
- for and
- the constants satisfy

int Kantorovitch(int m,VECTOR (* TheFunction)(VECTOR &),MATRIX (* Gradient)(VECTOR &), INTERVAL_MATRIX (* Hessian)(int, int, INTERVAL_VECTOR &),VECTOR &Input,double *eps)

`m`: number of variables and unknowns`TheFunction`: a procedure to compute the value of the equations for given values of the unknowns. This procedure has one arguments which is the value of the unknowns in vector form`Gradient`: a procedure to compute the Jacobian matrix of the system in matrix form. This procedure has one arguments which is the value of the unknowns in vector form`Hessian`: a procedure to compute the Hessian for the equation for interval value input. This procedure compute the`m X n`,`n`Hessian matrix in interval matrix form.This procedure has 3 arguments`l1,l2,X`. The function should return the value of the Hessian of the equations from`l1`to`l2`The Hessian of the first equation is stored in hess(1..`n`,1...`n`), the Hessian of the second equation in hess(`n`+1..2`n`,1..`n`) and so on`Input`: the value of the variables which constitute the center of the convergence ball

int Kantorovitch(int m, INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* Gradient)(int, int, INTERVAL_VECTOR &), INTERVAL_MATRIX (* Hessian)(int, int, INTERVAL_VECTOR &), VECTOR &Input,double *eps)There is also an implementation of Kantorovitch theorem for univariate polynomial, see section 5.2.12.

- = -1: the Jacobian matrix has no inverse at the mid-point of
`Input` - =0: the procedure has failed to determine a convergence ball
centered at
`Input` - =1: the procedure has found that an unique solution exist
within the interval [
`Input-eps,Input+eps`]

Rouche theorem

We will denote by the matrix of the derivatives of order of with respect to the variable. Let us consider a point and define as

and . If is strictly lower than , then has a single root in a ball centered at with radius and the Newton scheme with initial guess will converge to the solution.

The most difficult part for using this theorem is to determine
. For algebraic equations it is easy to determine a value
, that we will call the *order* of Rouche theorem, such that
and consequently may be
obtained by computing

for all in and taking as the Sup of all .

For non algebraic finding requires an analysis of the system.

Rouche theorem may be more efficient than Moore or Kantorovitvh theorems. For example when combined with a polynomial deflation (see section 5.9.6) it allows one to solve Wilkinson polynomial of order up to 18 with the C++ arithmetic on a PC, while stand solving procedure fails for order 13.

Rouche theorem is implemented in the following way:

- Rouche theorem is checked with respect to the mid-point of a box
- if Roucche theorem is satisfied, then a limited number of Newton iteration is performed to check if Newton indeed converge. If this is the case a ball that include a single solution has been determined
- if a ball has been determined, then, optionaly an inflation procedure )see section 3.1.6) is used to try to enlarge the ball

int Rouche(int DimensionEq,int DimVar,int order, INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_VECTOR (* Jacobian)(int, int, INTERVAL_VECTOR &), INTERVAL_MATRIX (* Gradient)(int, int, INTERVAL_VECTOR &), INTERVAL_VECTOR (* OtherDerivatives)(int, int, INTERVAL_VECTOR &), double Accuracy, int MaxIter, INTERVAL_VECTOR &Input, INTERVAL_VECTOR &UnicityBox)where

`DimensionEq`: number of equations`DimVar`: number of variables`order`: the order for Rouche theorem minus 1`TheIntervalFunction`: a procedure in`MakeF`format for computing an interval evaluation of the equations`Jacobian`: a procedure in`MakeF`format that computes the jacobian row by row`Gradient`: a procedure that compute the jacobian in`MakeJ`format`OtherDerivatives`: a procedure in`MakeF`format that computes the derivative of order larger or equal to 2, row by row. This procedure returns an interval vector of dimension

If a ball with a single solution has been found it will be returned in
`UnicityBox` and the procedure returns 1, otherwise it returns 0.

If the flag
`ALIAS_Always_Use_Inflation`
is set to 1, then an
inflation procedure is used to try to enlarge the box up to the
accuracy `ALIAS_Eps_Inflation`.

Interval Newton

The classical interval Newton method is embedded in the procedure
`GradientSolve` and `HessianSolve` but may also be useful in
other procedures. Furthermore this method relies on the use of the
product
where is the Jacobian of the system of
equations and the inverse of computed at some
particular point . In the classical method this product is
cimputed numerically and this does not take into account that the
element of
are functions of the same parameters. For example if the first column of
is where is some parameter with interval value, the
first element of
will be computer as

where are the elements of . Clearly the double occurence of in the numerical evaluation of the elements may lead to an overestimation of the elements: this element should be written as which is optimal in term of interval evaluation. Furthermore it may also be interesting to have the derivatives of each element of the product in order to improve the interval evalation of the matrix product. Indeed the interval evaluation of plays a very important role in the interval Newton method either for filtering a box for possible solution or for determining that a box includes a solution of the system.

The procedure `IntervalNewton` is a sophisticated interval Newton
algorithm that allows one to introduce knowledge on the product
in the classical scheme. Its syntax is:

int IntervalNewton(int Dim,INTERVAL_VECTOR &P,INTERVAL_VECTOR &FDIM, INTERVAL_MATRIX &Grad,MATRIX &GradMid, MATRIX &InvGradMid, int hasBGrad, INTERVAL_VECTOR (* BgradFunc)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* BgradJFunc)(int, int,INTERVAL_VECTOR &), int grad1, int grad3B1)where

`Dim`: the size of the system`P`: an interval vector that describes the range for the unknowns`FDIM`: interval value of the equation at the mid-point of`P``Grad`: interval jacobian at`P``GradMid`: jacobian at the mid-point of`P``hasBgrad`: a flag that indicates how will be calculated:- 0: will be calculated numerically
- 1: will be calculated using the procedure
`BgradFunc` - 2: will be calculated using the procedure
`BgradFunc`and the derivatives of the elements of available through the procedure`BgradJFunc`

`BgradFunc`: a user-provided procedure in`MakeF`format that calculate the element of , row by row`BgradJFunc`: a user-provided procedure in`MakeJ`format that calculate the derivatives of the elements of`grad1`: if`hasBgrad`is 1 or 2 we use the procedure`BgradFunc`to evaluate when we are in the 3B filter if`grad1`is set to 1. If set to 0 we use`BgradFunc`only when dealing with the full box`grad3B1`: if 1 we use the procedure that evaluates through`BgradFunc`even if we are in the 3B case. If 2 we use both`BgradFunc`and`BgradJFunc`

Various variants of `IntervalNewton` are available:

int IntervalNewton(int Dim,INTERVAL_VECTOR &P,INTERVAL_VECTOR &FMID, INTERVAL_MATRIX &Grad,MATRIX &GradMid,MATRIX &InvGradMid)is the classical interval Newton method with

int IntervalNewton(int Dim,INTERVAL_VECTOR &P,int DimVar,int DimEq, int TypeGradMid,MATRIX &InvGradMid, INTERVAL_VECTOR (*TheIntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* Gradient)(int, int, INTERVAL_VECTOR &))is also the classical interval Newton method for a system having

int IntervalNewton(int Dim,INTERVAL_VECTOR &P,int DimEq,int DimVar, int has_BGrad, INTERVAL_VECTOR (* BgradFunc)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* BgradJFunc)(int, int,INTERVAL_VECTOR &), int grad1,int grad3B1, int TypeGradMid, MATRIX &GradFuncMid, MATRIX &InvGradFuncMid, INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* Gradient)(int, int, INTERVAL_VECTOR &))Here the the mid jacobian

Miranda theorem

Miranda theorem provides a simple way to determine if there is one, or more, solution of a system of equations in a given box. It has the advantage of not requiring the derivatives of the equations but the drawback of not provinding the proof of the existence of a single solution in the box.

Let a system with . Let us consider a ball for and define

for in . If

or if

then has at least one zero in [15].

The simplest implementation of the Miranda theorem is

int Miranda(int Dim,INTERVAL_VECTOR (* F)(int,int,INTERVAL_VECTOR &), INTERVAL_VECTOR &Input)where

Another implementation uses the derivatives for improving the interval evaluation:

\begin{verbatim} int Miranda(int Dim, INTERVAL_VECTOR (* F)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* J)(int,int,INTERVAL_VECTOR &), INTERVAL_VECTOR &Input)

Inflation

Let a system , the Jacobian matrix of this system and a solution of the system. The purpose of the inflation method is to build a box that will contain only this solution. Let be a ball centered at : if for any point in is not singular, then the ball contains only one solution of the system.

The problem now is to determine a ball such for any point in the ball
the Jacobian is regular. Let be the matrix
whose
components are intervals. Let be the diagonal element of H having
the lowest absolute value, let be the maximum of the absolute
value of the sum of the elements at row of , discarding the
diagonal element of the row and let be the maximum of the
's. If , then the matrix is denoted *diagonally dominant*
and all the matrices are regular [19].

Let be a small constant: we will build incrementally the ball by using an iterative scheme defined as:

that will be repeated until is no more diagonally dominant. Note that in some cases (see for example section 2.16) it is possible to calculate directly the largest possible so that all the matrices in are regular without relying on the iterative scheme.

int ALIAS_Epsilon_Inflation(int Dimension,int Dimension_Eq, INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), INTERVAL_MATRIX (* Gradient)(int, int, INTERVAL_VECTOR &), INTERVAL_MATRIX (* Hessian)(int, int, INTERVAL_VECTOR &), VECTOR &X0, INTERVAL_VECTOR &B)Note that the

2018-07-25