next up previous contents Next: Linear algebra Up: ALIAS-C++ Previous: Analyzing univariate polynomials

  • La page de J-P. Merlet
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    Subsections


    Parametric polynomials and eigenvalues of parametric matrices

    A parametric polynomial is a polynomial whose coefficients are functions of a set of parameters (in other words it is a set of polynomials). A typical parametric polynomial is obtained when calculating the characteristic polynomial of a parametric matrix.

    In this chapter we propose some algorithms to deal with the real roots of a parametric polynomial and, in some cases, with the real part of these roots. If the considered polynomial is the characteristic polynomial of a matrix we may make use of the components of the matrix.

    Some of these algorithms use a primary and secondary algorithms. The secondary algorithm uses also a list of boxes which is stored in the interval matrix BoxUP.


    Minimal and maximal real roots of a parametric polynomial

    The purpose here is to determine the minimal and maximal values of the set of real roots of the the set of polynomials, when the parameters are bounded. There may be eventually also constraints on the parameters. The algorithm is implemented as:

     
    int ALIAS_Min_Max_EigenValues(int Degree,  
                                  int Nb_Parameter,
                                  INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), 
                                  int Nb_Constraints,
                                  INTEGER_VECTOR &Type_Eq,
                                  int (* TheMatrix)(INTERVAL_VECTOR &, INTERVAL_MATRIX &), 
                                  int Has_Matrix,
                                  INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), 
                                  int *Has_Gradient,
                                  INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
                                  INTERVAL & TheDomain, 
                                  INTERVAL_VECTOR & TheDomain_Parameter, 
                                  int Type,int Nb_Points,int Use_Solve,int rand,
                                  int M,
                                  double Accuracy_Variable,double Accuracy,double AccuracyM,
                                  INTERVAL &Lowest_Root,INTERVAL &Highest_Root,
                                  INTERVAL_MATRIX &Place,int  Stop,double *Seuil,
                                  int (* Solve_Poly)(double *, int *,double *), 
                                  int (* Simp_Proc)(INTERVAL_VECTOR &))
    
    where the arguments are: The confidence in this procedure is at the same level than the confidence in the numerical algorithm that solve a polynomial.

    The return code is:

    The possible bisection mode are:

    For mode 2 if the gradient is available and if the polynomial parameter has a value that is better than the current optimum, then this variable is bisected otherwise we use the smear function to determine the bisected variable.

    Possible parameters values for a given range on the real roots

    It may be of interest to determine what may be the possible values of the parameters of a parametric polynomial such that the real roots of the corresponding polynomial are all enclosed in a given interval. ALIAS provides two routines for that purpose.

    Approximation of the set of solutions

    Let consider the parameters space i.e. a $n$ dimensional space where each of the dimension corresponds to one of the $n$ parameters. A point in this space corresponds to a unique value for all the parameters and therefore to a specific polynomial. In the parameters space there are possibly a set ${\cal S}$ of regions such that for any point in the region(s) the corresponding polynomial has all its root within the given interval. The purpose of the following procedure is to determine an approximation of ${\cal S}$. This approximation ${\cal A}$ will be constituted of a set of $n$ dimensional boxes which are guaranteed to be included in ${\cal S}$ and that will be written in a file. During the calculation the boxes whose width is lower than a given threshold $\epsilon$ and for which the algorithm has been unable to determine if they are fully enclosed in ${\cal S}$ will be neglected. A possible index for measuring the quality of the approximation ${\cal A}$ is the ratio $V/V_n$ between the total volume $V$ of the boxes written into the file over the total volume $V_n$ of the boxes that have been neglected as the volume of ${\cal S}$ is lower or equal to $V+V_n$.

    The procedure is:

     
    int ALIAS_Min_Max_EigenValues_Area(int Degree,int Nb_Parameter,
              int Has_Interval,
              INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), 
              INTERVAL_VECTOR (* TheCoeffCentered)(INTERVAL_VECTOR &,double), 
              int Nb_Constraints,INTEGER_VECTOR &Type_Eq,
              int (* TheMatrix)(INTERVAL_VECTOR &, INTERVAL_MATRIX &), 
              int Has_Matrix,
              INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), 
              int Has_Gradient,
              INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
              INTERVAL & TheDomain,INTERVAL_VECTOR & TheDomain_Parameter, 
              int Nb_Points,int Use_Solve,int rand,int Strong,int Iteration,
              double Accuracy_Variable,double Accuracy,double AccuracyM,double AccuracyB,
              double *Volume_Result,double *Volume_Neglected,double Seuil,
              char *FileName,int Has_Input,char *File_Input,
              int (* Solve_Poly)(double *, int *,double *),int RealRoot,
              INTERVAL_VECTOR (* Evaluate_Complex)(int,int,INTERVAL_VECTOR &), 
              int (* Simp_Proc)(INTERVAL_VECTOR &))
    
    where the arguments are similar to the one of the previous procedure except for:

    This procedure returns the number of boxes written in the result file or a negative number if the calculation has failed. The possible negative return code are:

    Largest square enclosed in the regions

    This second procedure allows to compute the largest square (up to a pre-defined accuracy $\alpha$) that is enclosed in the region. This largest box can clearly be obtained from the result of the previous algorithm but this weaker procedure will be faster.

    The principle is to have a set of boxes whose first element is a range for the center of the box and second element is a possible value for the length of the half-edge of the square. The main algorithm will test if for a given box the center may be a candidate to be the center of a square of half-edge $R+\alpha$ (where $R$ is the current optimum) using some heuristics and if the answer is positive will compute the minimal and maximal root of the polynomials defined by this square using a secondary algorithm. If these values are compatible with the bound the current optimum will be updated.

     
    int ALIAS_Geometry_Carre(int Degree,int Nb_Parameter,
               INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), 
               int Nb_Constraints,INTEGER_VECTOR &Type_Eq,INTEGER_VECTOR &Imperatif,
               int (* TheMatrix)(INTERVAL_VECTOR &, INTERVAL_MATRIX &), 
               int Has_Matrix,
               INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &),  
               int *Has_Gradient,
               INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
               INTERVAL & TheDomain,INTERVAL_VECTOR & TheDomain_Parameter, 
               int Nb_Points,int Use_Solve,int rand,
               int Iteration_Geometry,int Iteration_Polynom,
               double Accuracy_Variable,double Accuracy,double Accuracy_Geometry,
               double Accuracy_Polynom,INTERVAL_VECTOR &Solution,double *Seuil,
               int (* Solve_Poly)(double *, int *,double *), 
               int (* Simp_Proc)(INTERVAL_VECTOR &),
               int (* Simp_Proc_Pol)(INTERVAL_VECTOR &))
    
    The arguments are the same than for the previous procedure except for:

    This procedure will return:

    During the calculation the flag ALIAS_Has_OptimumG will be set to 1 as soon as an optimum is found: the center of the current optimal geometry is the mid point of the interval vector ALIAS_Vector_OptimumG while is edge is in the interval ALIAS_OptimumG.

    The secondary algorithm uses the flag Single_BisectionG and Reverse_StrorageG that play the same role than Single_Bisection and Reverse_Storage in the general solving algorithms (see sections 2.3.1.3 and 2.3.1.2).

    Condition number

    The condition number of a polynomial may be defined either as the ratio lowest root over largest root or as the ratio ${\rm
Min}(\vert x_i\vert)$ over ${\rm Max}(\vert x_i\vert)$ where $x_i$ are the roots of the polynomial. In the later case the condition number has a value between 0 and 1. The minimal and maximal values of the condition number of a parametric polynomial in both form may be calculated using the procedure:

     
    int ALIAS_Min_Max_CN(int Degree,  
         int Nb_Parameter,INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), 
         int Nb_Constraints,INTEGER_VECTOR &Type_Eq,
         int (* TheMatrix)(INTERVAL_VECTOR &, INTERVAL_MATRIX &), 
         int Has_Matrix,
         INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), 
         int Has_Gradient,
         INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
         INTERVAL & TheDomain,INTERVAL_VECTOR & TheDomain_Parameter, int Type,
         int Nb_Points,int Absolute,
         int rand,int Iteration,
         double Accuracy_Variable,double Accuracy,double AccuracyM,
         INTERVAL &Lowest,INTERVAL &Highest,
         INTERVAL_MATRIX &Place,int  Stop, double *Seuil,
         int (* Solve_Poly)(double *, int *,double *), 
         int (* Simp_Proc)(INTERVAL_VECTOR &))
    
    where the arguments are identical than for the previous procedure except for: This procedure allows for the calculation of the minimal and maximal condition number of a matrix. Various bisection methods are available through the use of the integer Single_Bisection: A specific procedure is used when the gradient is available. If the interval for the root includes 0 we bisect it. Otherwise we use the smear function to decide which other variable should be bisected.

    Kharitonov polynomials

    Kharitonov polynomials are special polynomials that have constant values for their coefficients and are associated to a parametric polynomial. It can be shown that if all the Kharitonov polynomials have the real parts of their roots of the same sign, then all the polynomials in the set will have the real part of their roots of the same sign. For a polynomial $P= \sum C(i)x^{i-1}$ the four Kharitonov polynomials are:

    \begin{eqnarray*}
&&P-=Inf(C1)+Inf(C2)x+Sup(C3)x^2+Sup(C4)x^3+ Inf Inf Sup Sup\...
...&P++=Sup(C1)+Sup(C2)x+Inf(C3)x^2+Inf(C4)x^3+ Sup Sup Inf Inf \\
\end{eqnarray*}

    Implementation

    The following procedure allows to determine if a polynomial has a real root within a given range:
     
    int Kharitonov(int Degree,
           INTERVAL_VECTOR (* TheCoeff)(double a,INTERVAL_VECTOR &), 
           int (* Solve_Poly)(double *, int *,double *), 
           INTERVAL_VECTOR &Input)
    
    where: Note that there exists an ALIAS-Maple procedure that uses this procedure to generate a simplification procedure, see the ALIAS-Maple manual.

    Gerschgorin circles

    Mathematical background

    Let $A=((a_{ij}))$ and the Gerschgorin circles defined as the set of $z$ such that:

    \begin{displaymath}
C_i =\vert z-a_{ii}\vert \le \sum_j \vert a_{ji}\vert, j \not= i
\end{displaymath}

    The roots of the characteristic polynomial are enclosed in the union of the $C_i$. Furthermore if a circle has no intersection with the other circles, then this circle contains one root of the polynomial. This allow for a fast determination of bounds for the real roots of the polynomial.

    Implementation

    It is possible to get these bounds by using the procedure:

     
    int Gerschgorin(INTERVAL_MATRIX &A,int Size, int Type, INTERVAL &Bound)
    
    where: the procedure returns 1 if it has been possible to compute Bound.

    A more complete procedure allows to get all the Gerschgorin circles and eventually to adjust an interval that is supposed to contain the roots of the polynomial:

     
    int Gerschgorin_Simplification(INTERVAL_MATRIX &A,int Size,int Type,
                           INTERVAL &Input,INTERVAL_VECTOR &Circle)
    
    The arguments are the same than the previous procedure except for Input which is the interval supposed to contain roots of the polynomial and Circle which contain the projection of the Gerschgorin circle on the real axis. This procedure returns:


    Cassini ovals

    Mathematical background

    The Cassini ovals are another method to determine a bound for the eigenvalues of a matrix. Let $A=((a_{ij}))$ and the Cassini ovals defined as the set of $z$ such that:

    \begin{displaymath}
O_{jk} =\vert z-a_{jj}\vert\vert z-a{kk} \le (\sum_i \vert a_{ji}\vert, i \not= j)(\sum_i
\vert a_{ki}\vert, i \not= k)
\end{displaymath}

    Column based Cassini ovals may also be defined. The roots of the characteristic polynomial are enclosed in the union of the row-based and column-based $O_{jk}$. Although more complicate to calculate the bounds obtained with the Cassini ovals are usually tighter than the bounds obtained with the Gerschgorin circles.

    Implementation

    The Cassini bounds may be obtained with the procedure:

     
    int Cassini_Simplification(INTERVAL_MATRIX &A,int Size, 
                      int Type_Matrix, INTERVAL &Input)
    
    where The procedure will return -1 if there is no intersection between Input and the Cassini bounds, 1 if there Input has been improved, 0 otherwise.

    This procedure may also be called with:

     
    int Cassini_Simplification(int (* TheMatrix)(INTERVAL_VECTOR &, INTERVAL_MATRIX &),
          INTERVAL_VECTOR &Param,int Size, int Type_Matrix, INTERVAL &Input)
    
    where TheMatrix is a procedure that compute the elements of the matrix according to the interval value of the parameters Param.


    Routh

    The Routh algorithm allows to calculate the number of roots with positive real part of a polynomial being given the coefficients of the polynomial. It is implemented as:

     
    int Routh(int Degree,double *Coeff)
    
    that returns the number of roots with positive real part or -1 if it has not been possible to compute the Routh table (because the first element of a row of the Routh table is close to 0).

    A similar algorithm allows to deal with polynomial with interval coefficients:

     
    INTERVAL Routh(int Degree,INTERVAL_VECTOR &Coeff)
    
    This algorithm returns in its interval:
    -1,-1
    : the Routh table cannot be computed
    a,a+1
    : there is at least a roots with positive real part but the exact number of roots with positive real part cannot be calculated
    a,a
    : there is exactly a roots with positive real part
    A similar procedure may be used when the coefficients are functions of parameters:
     
    INTERVAL Routh(int Degree,INTERVAL_VECTOR (* TheCoeff)(int,int,INTERVAL_VECTOR &),
                   INTERVAL_VECTOR &Input)
    
    where: If the derivatives of the coefficients with respect to the parameters are known an even better procedure will be:
     
    INTERVAL Routh(int Degree,INTERVAL_VECTOR (*TheCoeff)(int,int,INTERVAL_VECTOR &),
                   INTERVAL_MATRIX (* TheCoeffG)(int,int,INTERVAL_VECTOR &),
                   INTERVAL_VECTOR &Input)
    
    where TheCoeffG is a procedure that allow to compute the derivatives of the coefficients with respect to the parameters (see note 2.4.2.2). This procedure allows to a certain amount to take into account the dependency between the coefficients.

    Note also that the procedure Routh of ALIAS-Maple allows an even better calculation of the Routh table when dealing with parametric polynomials as the elements of the Routh table are calculated symbolically.


    Weyl filter

    Mathematical background

    Let $P$ be a polynomial and be maxroot the maximal modulus of the root of $P$. From $P$ we may derive a the unitary polynomial $Q$ such that the roots of $Q$ have a modulus lower or equal to 1 and if $w$ is a root of $Q$ then maxroot$w$ is a root of $P$.

    Let $Q=sum_{i=0}^{i=n} a_i x^i$ which may also be written as $sum_{i=0}^{i=n} b_i(x-a)^i$ where $a$ is some fixed point.

    Let a range $[a,b]$ for $x$ and let $z_0$ be the mid point of the range. We consider the square in the complex plane centered at $z_0$ and whose edge length is $b-a$. Let $\delta$ be the length of the half-diagonal of this square. If

    \begin{displaymath}
\vert b_0\vert> \sum_{j=1}^{j=n} \vert b_j\vert\delta^j
\end{displaymath}

    then the polynomial has no root in the square [1]. This procedure may be used for univariate polynomial, with polynomial with interval coefficients or with parametric poynomial. For example it is very efficient for solging the Wilkinson polynomial.

    Implementation

    The Weyl filter is implemented to manage univariate polynomial (with numerical or interval coefficients) or parametric polynomial. A basic tool of this filter is

     
    int Weyl_Filter_Utility(int Degree, double maxroot,
                          INTERVAL_VECTOR &b,INTERVAL_VECTOR &Input)
    
    where Input(1) is the range for the polynom variable and b the $b_i$ coefficients. A more general implementation is
     
    int Weyl_Filter(int Degree, INTERVAL_VECTOR &Coeff,double maxroot,
                    int (* GetB)(int,INTERVAL_VECTOR &,INTERVAL &,INTERVAL_VECTOR &), 
    		INTERVAL_VECTOR &Input)
    
    Here GetB is a procedure that calculates the $b_i$ coefficients. Its arguments are the degree of the polynomial, its coefficients, $z_0$ and the $b_i$ coefficients. It should return 1 if the calculation has been successfull, 0 otherwise. An example of such procedure is Derive_Polynomial_Expansion. Degree is the degree of the polynomial, Coeff its coefficients and Input(1) is the range for the polynomial unknown.

    This procedure has several variants:

     
    int Weyl_Filter(int Degree, INTERVAL_VECTOR &Coeff,double max_root,
                    INTERVAL_VECTOR &Input)
    
    here the GetB procedure will be Derive_Polynomial_Expansion. In
     
    int Weyl_Filter(int Degree, INTERVAL_VECTOR &Coeff,
                    INTERVAL_VECTOR &Input)
    
    here the GetB procedure will be Derive_Polynomial_Expansion and maxroot will be computed by an ALIAS procedure. The same procedure will be used in
     
    int Weyl_Filter(int Degree, INTERVAL_VECTOR &Coeff,
                    int (* GetB)(int,INTERVAL_VECTOR &,INTERVAL &,INTERVAL_VECTOR &), 
                    INTERVAL_VECTOR &Input)
    

    Coefficient of the characteristic polynomial

    In some cases it may be difficult to obtain the analytical form of the coefficients of the characteristic polynomial of a matrix. The following procedure allows to compute these coefficients even if the matrix is an interval matrix:

     
    INTERVAL_VECTOR Coeff_CharPoly(int Size,INTERVAL_MATRIX &A)
    
    where Size is the dimension of the matrix and A the matrix. This procedures returns the coefficients in an interval vector $C$ and the polynomial is $\sum C(i)x^{i-1}$.


    next up previous contents Next: Linear algebra Up: ALIAS-C++ Previous: Analyzing univariate polynomials
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    jean-pierre merlet
    2018-07-25