Implementation

The algorithm we have implemented is a direct derivation of the general purpose solving algorithm with Jacobian and Hessian in which these values are automatically derived. To isolate the roots we use the Kantorovitch theorem (which may also optionally be used during the resolution, see section 3.1.2). To eliminate boxes during the bisection process we use the safe Budan-Fourier method (see section 5.5.2).

 
int Solve_UP_JH_Interval(int Degree,VECTOR Coeff, 
            INTERVAL & TheDomain, 
            int Order,int M,int Stop,
            double epsilon,double epsilonf,
            INTERVAL_VECTOR & Solution,
            INTEGER_VECTOR & IsKanto,int NbSolution);

with:

• Degree: degree of the polynomial
• Coeff: the Degree+1 coefficients of the polynomial in increasing degree
• TheDomain: the interval in which we are looking for roots
• Order: the type of order which is used to store the intervals created during the bisection process. This order may be either MAX_FUNCTION_ORDER or MAX_MIDDLE_FUNCTION_ORDER. See the note on the order 2.3.4.4.
• M: the maximum number of boxes which may be stored. See the note 2.5.2.2
• Stop: the possible values are 0,1,2
  • 0: the algorithm will look for every solution in TheDomain
  • 1: the algorithm will stop as soon as 1 solution has been found
  • 2: the algorithm will stop as soon as Nb solutions have been found
• epsilon: the maximal width of the box, see the note 2.3.4.6
• epsilonf: the maximal width of the equation intervals, see the note 2.3.4.6
• Solution: an interval matrix of size (Nb,m) which will contained the solution intervals.
• IsKanto: an integer vector of dimension Nb. A value of 1 for IsKanto(i) indicate that Newton method (see section 2.9) with as estimate the center of some solution interval Solution(i) has been used and has converged toward the unique solution Solution(i) which lie within this solution intervals. Note that the interval which contain the solution may be retrieved in the interval vector Interval_Solution_UP.
• NbSolution: the maximum number of solution we are looking for.

The procedure will return the number of solution(s) or -3 if the order is not 0 or 1, -2 if the number of equations or unknowns is equal or lower to 0 and -1 if the number of iteration is too low. There is an alternate form of this procedure in the case where we are looking for all the roots of the polynomial.

 
int Solve_UP_JH_Interval(int Degree,VECTOR Coeff,
            int Order,int M,int Stop,
            double epsilon,double epsilonf,
            INTERVAL_VECTOR & Solution,
            INTEGER_VECTOR & Is_Kanto,int NbSolution);

There are two alternate forms of this procedure in the case where we are looking for the positive or negative roots of the polynomial.

 
int Solve_UP_JH_Positive_Interval
int Solve_UP_JH_Negative_Interval

In the three previous procedures there is no TestDomain as it is automatically determined by the procedure. If there was a failure in the determination of the domain (for the reasons explained in section 5.2) the procedures will return -1.

The previous procedures are numerically safe in the sense that we take into account rounding errors in the evaluation of the polynomial and its gradient. For well conditioned polynomials you may use faster procedures whose name has the prefix Fast. For example Fast_Solve_UP_JH_Interval is the general procedure for finding the roots of a polynomial.

Clearly this procedure is not intended to be used as substitute to more classical algorithms.

It makes use of a specific Krawczyk procedure for polynomials:

 
int Krawczyk_UP(int Degree,INTERVAL_VECTOR &Coeff, 
             INTERVAL_VECTOR &CoeffG,INTERVAL &Input)