home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / zeroin.cc < prev   
Encoding:
Text File  |  1995-12-04  |  4.9 KB  |  142 lines  |  [TEXT/ttxt]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *
  7.  *                Brent's root finder
  8.  *           obtains a zero of a function of one variable
  9.  *
  10.  * Synopsis
  11.  *    double zeroin(ax,bx,f,tol=EPSILON)
  12.  *    const double ax         The root is to be sought within
  13.  *    const double bx          the interval [ax,bx]
  14.  *    double (*f)(const double x)    Ptr to the function under 
  15.  *                    consideration    
  16.  *    const double tol        Acceptable tolerance for the root
  17.  *                    position. It is an optional parameter
  18.  *                    with default value EPSILON
  19.  *
  20.  *    Zeroin returns an approximate location for the root with accuracy
  21.  *    4*EPSILON*abs(x) + tol
  22.  *
  23.  * Algorithm
  24.  *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
  25.  *    computations. M., Mir, 1980, p.180 of the Russian edition
  26.  *
  27.  * The function makes use of the bissection procedure combined with
  28.  * the linear or quadratic inverse interpolation.
  29.  * At every step program operates three abscissae - a, b, and c.
  30.  *    b - the last and the best approximation to the root
  31.  *    a - the last but one approximation
  32.  *    c - the last but one or even earlier approximation such that
  33.  *        1) |f(b)| <= |f(c)|
  34.  *        2) f(b) and f(c) have opposite signs, i.e. b and c encompass
  35.  *           the root
  36.  * At every step Zeroin computes two new approximations, one by the 
  37.  * bissection procedure and the other one from interpolation (if a,b, and c
  38.  * are all different the quadratic interpolation is used, linear otherwise).
  39.  * If the latter (i.e. obtained by the interpolation) point looks
  40.  * reasonable (i.e. falls within the current interval [b,c] not close
  41.  * to the end points of the interval), the point is accepted as a new
  42.  * approximation to the root. Otherwise, the bissection result is used.
  43.  * Therefore, the range of uncertainty is guaranteed to be reduced at 
  44.  * least by the factor of 1.6
  45.  *
  46.  ************************************************************************
  47.  */
  48.  
  49. #pragma implementation "math_num.h"
  50. #include "math_num.h"
  51.  
  52.  
  53. double zeroin(                // An estimate to the root
  54.     const double ax,        // Specify the interval the root
  55.     const double bx,        // to be sought in
  56.     double (*f)(const double x),    // Function under investigation
  57.     const double tol)        // Acceptable tolerance
  58. {
  59.   double a = ax, b = bx, c;        // Abscissae, see above
  60.   double fa;                // f(a)
  61.   double fb;                // f(b)
  62.   double fc;                // f(c)
  63.  
  64.   assure( tol > 0, "Tolerance must be positive");
  65.   assure( b > a, 
  66.      "Left end point of the interval should be strictly less than the "
  67.      "right one" );
  68.  
  69.   fa = (*f)(a);  fb = (*f)(b);
  70.   c = a;   fc = fa;
  71.  
  72.   for(;;)        // Main iteration loop
  73.   {
  74.     double prev_step = b-a;        // Step from the previous iteration
  75.    
  76.     if( fabs(fc) < fabs(fb) )
  77.     {                                 // Swap data so that b would be the
  78.       a = b;  b = c;  c = a;              // best approximation found so far
  79.       fa=fb;  fb=fc;  fc=fa;
  80.     }
  81.                     // Estimate the effective tolerance
  82.     const double tol_act = 2*EPSILON*fabs(b) + tol/2;
  83.     double new_step = (c-b)/2;        // Bissection step for this iteration
  84.  
  85.     if( fabs(new_step) <= tol_act || fb == 0 )
  86.       return b;                // Acceptable approximation is found
  87.  
  88.                 // Figuring out if the interpolation can be tried
  89.     if( fabs(prev_step) >= tol_act    // If prev_step was large enough
  90.     && fabs(fa) > fabs(fb) )    // and was in true direction,
  91.     {                    // Interpolatiom may be tried
  92.  
  93.       double p;                  // Interpolation step is calcu-
  94.       double q;                  // lated in the form p/q; divi-
  95.                       // sion operations is delayed
  96.                      // until the last moment
  97.       const double cb = c-b;
  98.  
  99.       if( a==c )            // If we've got only two distinct
  100.       {                    // points linear interpolation
  101.     register double t1 = fb/fa;    // can only be applied
  102.     p = cb*t1;
  103.     q = 1.0 - t1;
  104.       }
  105.       else                // Quadratic inverse interpolation
  106.       {
  107.     register double t1, t2;
  108.     q = fa/fc;  t1 = fb/fc;  t2 = fb/fa;
  109.     p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
  110.     q = (q-1.0) * (t1-1.0) * (t2-1.0);
  111.       }
  112.  
  113.       if( p > 0 )            // Formulas above computed new_step
  114.     q = -q;                // = p/q with wrong sign (on purpose).
  115.       else                // Correct this, but in such a way so
  116.     p = -p;                // that p would be positive
  117.       
  118.       if( p < (0.75*cb*q-fabs(tol_act*q)/2)    // If b+p/q falls in [b,c]
  119.      && p < fabs(prev_step*q/2) )    // and isn't too large
  120.     new_step = p/q;            // it is accepted
  121.                     // If p/q is too large then the
  122.                     // bissection procedure can
  123.                     // reduce [b,c] to a larger
  124.                     // extent
  125.     }
  126.  
  127.     if( fabs(new_step) < tol_act )    // Adjust the step to be not less
  128.       if( new_step > 0 )        // than the tolerance
  129.     new_step = tol_act;
  130.       else
  131.     new_step = -tol_act;
  132.  
  133.     a = b;  fa = fb;            // Save the previous approximation
  134.     b += new_step;  fb = (*f)(b);    // Do step to a new approximation
  135.     if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) )
  136.     {                             // Adjust c for it to have the sign
  137.       c = a;  fc = fa;                  // opposite to that of b
  138.     }
  139.   }
  140.  
  141. }
  142.