next up previous contents Next: Integration Up: ALIAS-C++ Previous: Optimization

  • 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


    Continuation for one dimensional system

    Although mostly devoted to zero-dimensional systems ALIAS still offers some possibilities to deal with $n$-dimensional systems. The purpose of the algorithms here is to consider a system with $n$ parameter and to follow the branches of the solutions of the system when these parameters are varying between some bounds.


    Continuation 1D

    Mathematical background

    Let consider a system of $n$ equations in $m$ unknowns $F_i(x_1,\ldots,x_m)=0$ and assume that $x_1$ may vary in the range $[\underline{x_1},\overline{x_1}]$ and that it is considered as the parameter of this system. Using the solving algorithms of ALIAS (or any other method) we are able to determine the real roots of this system, if any, for a given value of the parameter, for example for $x_1=\underline{x_1}$. Let assume that we have found $k$ real roots $X_1,\ldots, X_k$. Assume now that we change the value of $x_1$ to $\underline{x_1}+\epsilon$. Using Kantorovitch theorem (see section 3.1.2) we are able to determine if Newton method with as initial estimate $X_i$ will converge toward a solution of the new system and the radius of convergence $\alpha_i$ such that the solution is unique in the ball $B_i =[X_i-\alpha_i, X_i+\alpha_i]$ (alternatively we may also use Moore theorem, see section 2.10). If this is the case and if the intersection of the $B_i$ is empty, then we are able to compute in a certified way the new $k$ solutions of the new system. Otherwise we will repeat the procedure with a smaller value for $\epsilon$ until the method succeed. A failure case will be when the system become singular at a point (for example when 2 branches are collapsing). If this happen we will change the value of $x_1$ until we find a new set of solutions which satisfy Kantorovitch theorem and start again the process.

    This method enable to follow all the branches for which initial points have been found when solving the initial system. We may then store these points in an array, for example at regular step for the value of the parameter.

    By default we use Kantorovitch theorem to follow the branches but by setting the global variable
    ALIAS_Use_Kraw_Continuation to 1 we will use Moore test.

    Implementation

    Certified Newton

    The procedure that allow to determine the next point on a branch is implemented as: 

     
int Certified_Newton(int Nb_Var,int Dimension_Eq,int Nb_Branch,int Branch,
                       INTERVAL_VECTOR (* Func)(int,int,INTERVAL_VECTOR &), 
                       INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
                       INTERVAL_MATRIX (* Hessian)(int, int, INTERVAL_VECTOR &), 
                       double Accuracy_Var,double Accuracy,
                       double *param,double delta_param,double min_delta_param,
                       int sens,MATRIX &Solution)
    This procedure determine if it possible to find a new point on the branch numbered Branch among the Nb_Branch that have been found. In Solution we have all the solutions that have been found for the branches when the parameter value was param. The change in the parameter value will be at most sens$\times$delta_param and on exit of the procedure param will be set to the new value. sens indicates the direction of change in the parameter value: +1 for an increase, -1 for a decrease.

    This procedure returns:

    Procedure for following branches

    This procedure takes as input a set of $n$ solutions of the system and will return points on the $n$ branches. The branches will be followed until a given value for the parameter is reached or if Kantorovitch theorem is no more satisfied for some value of the parameter. It is implemented as: 

     
int ALIAS_Start_Continuation(int m,int n,int NUM,
       INTERVAL_MATRIX &Solutions,
       INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), 
       INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
       INTERVAL_MATRIX (* IntervalHessian)(int,int,INTERVAL_VECTOR &),
       double epsilon,double epsilonf,
       double *z,double delta,double mindelta,
       INTERVAL &Rangez,
       int sens,
       MATRIX &BRANCH,int *NBBRANCH)
    the arguments being: The return code is: Even for a negative return code the points in BRANCH are valid: the procedure will return the points that have computed on the $n$ branches but is simply not able to determine the location of new branches that may have appeared.


    Full continuation procedure

    This procedure will determine initial points on the branches of the system for the value of the parameter within some bounds and then follow the branches. It is implemented as 
     
int ALIAS_Full_Continuation(int m,int n,
       INTERVAL_VECTOR (* IntervalFunction)(int,int,INTERVAL_VECTOR &), 
       INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &), 
       INTERVAL_MATRIX (* IntervalHessian)(int,int,INTERVAL_VECTOR &),
       INTERVAL_VECTOR &Domain,
       int M,
       double epsilon,double epsilonf,
       double *z,double delta,double mindelta,double mindz,
       INTERVAL &Rangez,
       int sens,
       MATRIX &BRANCH,int *NBBRANCH)
    The arguments are the same than for the previous procedure except for: To determine the starting points of the branches this procedure uses the solving procedure with the Hessian (see section 2.5). As soon as initial valid solutions are found the ALIAS_Start_Continuation procedure is called until either a bound of the parameter range is reached, or a singularity occur. In the later case the solving procedure is called with an increased value for the parameter. This algorithm cannot find isolated points and may miss branches (see next section).

    There is also another version of this program where you indicate just before Domain the solutions which have already been found. The syntax is 

     
       INTERVAL_MATRIX (* IntervalHessian)(int,int,INTERVAL_VECTOR &),
       int NUM,
       INTERVAL_MATRIX &Solutions,
       INTERVAL_VECTOR &Domain


    The return code for these procedures are:


    Finding the initial starting point with the accuracy mindz may be computer intensive. Hence the integer global variable ALIAS_Allow_Backtrack enable to disable this process if it is set to 0 (its default value is 1): in that case as soon as starting points have been found (hence at $z_0$+ $n$delta) we will start following the branches.


    In fact these procedures are special occurrences of another ALIAS procedure which has another argument right after the Hessian argument. Assume for example that you are considering a system which has one equation written as:

    \begin{displaymath} (x-1)^2+(y-1)^2=l_{25} \end{displaymath}

    where $l_{25}$ is the parameter of the system and $x,y$ the unknowns. When using the continuation method we have to define ranges for these unknowns in order to be able to solve the system of equations. Up to now we have indicated bounds that are constants but for the equation example it will be interesting to be able to specify that these bounds may change according to the value for the parameter using a simplification procedure. In our example clearly $x$ and $y$ cannot exceed $l_{25}+1$ and cannot be lower than $-l_{25}+1$ (this is an application of the concept of 2B-consistency, see section 2.17). Hence we may specify right after the Hessian argument the name of a procedure, for example Range, that is able to update the bounds on the unknowns according to the value of the parameter (or any other variable that may play a role). The syntax of the procedure Range is: 
    INTERVAL_VECTOR Range(double z, INTERVAL_VECTOR &Variable)
    where z is parameter value and Variable the current set of ranges for the unknowns. This procedure must return a set of ranges for the unknowns (be careful to check that the returned ranges are included in the initial ranges).

    Note also that the ALIAS-Maple package offers a procedure that uses the method described in section 2.14 for finding the initial starting points of the branches: this method is efficient if the equations include linear terms in the unknowns.

    Missed branches

    Clearly there is a major problem with the method we are proposing: we may miss some branches. For example imagine that a system has $n$ roots for the initial value of the parameter, but will have more than $n$ solutions for another value of the parameter, meaning that new branches appear: our algorithm will find only the $n$ initial branches. There are two methods that can be used to find the correct number of branches. The first one is simply to start following the branches with as initial value for the parameter the one among the highest or smallest value having led to the maximum number of solution. There is also another mechanism that enable to avoid missing branches. Assume that the solving procedure has determined for some initial value $\lambda_0$ of the parameter $n$ solutions to the system and that at some point $\lambda_1$ the continuation method has failed: the solving procedure is called and determine that the system has now $m>n$ solutions. This mean that for some parameter value between $\lambda_0$ and $\lambda_1$ we have missed $m-n$ branches. At such point, called problem point it would be interesting to backtrack i.e. to start again a continuation process with as initial point for the branches the $m$ solutions obtained for $\lambda_1$ and a value for Sens which is the opposite of the initial value. This is what is done by the procedure which may store up to 10 problem points. As this process may be computer intensive it is possible to disable it by setting the integer global variable ALIAS_Problem_Continuation to -1 (it's default value is 0, which mean that the process is enabled).

    Example

    Let consider the system:

    \begin{eqnarray*} &&x^2-(y-1)^2+z^2-1=0\\ &&x^2-y^3=0 \end{eqnarray*}

    The purpose of this example is to write the points of the branches of this system in files whose name is BRANCH followed by the branch number. The main file is written as: 
     
#include "Functions.h"
#include "Vector.h"
#include "IntervalVector.h"
#include "IntervalMatrix.h"
#include "IntegerMatrix.h"
#include "header_Solve_General_Interval.h"
#include "header_Utilities_Interval.h"
double z;
//F,J,H are the equation, jacobian and hessian procedures
#include "F.C"
#include "J.C"
#include "H.C"

int main()
{
MATRIX BRANCH;
INTERVAL RANGE;
int NUM,i,j,NB_BRANCH;
char texte[400];
double delta=0.05;
double min_delta=1.e-6;
double mindz=1.e-4;
ostream out;

int Nb_Var=2;
int Nb_Eq=2;
Domain(1)=INTERVAL(-100,100);
Domain(2)=INTERVAL(-100,100);
RANGE=INTERVAL(0.2);
NUM=ALIAS_Full_Continuation(Nb_Var,Nb_Eq,F,J,H,Domain,1.e-6,1.e-6,&z,
	delta,min_delta,mindz,RANGE,1,BRANCH,&NB_BRANCH);

//order the branch and write the result in the file
for(i =1;i<=NUM;i++)
     {
       sprintf(texte,"BRANCH%d",i);
       out.open(texte,ios::out);
       for(j =1;j<=NB_BRANCH;j++)
          {
           if(BRANCH(j,Nb_Var+2)!=i)continue;
           out<<BRANCH(j,1)<<" "<<BRANCH(j,2)<<" "<<BRANCH(j,3)<<endl;
          }
      out.close();
      }
}

    next up previous contents Next: Integration Up: ALIAS-C++ Previous: Optimization
  • 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