Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

zeroin.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                        Numerical Math Package
00006  *
00007  *                          Brent's root finder
00008  *             obtains a zero of a function of one variable
00009  *
00010  * Synopsis
00011  *      double zeroin(ax,bx,f,tol=EPSILON)
00012  *      const double ax                 The root is to be sought within
00013  *      const double bx                 the interval [ax,bx]
00014  *      UnivariateFunctor& f            The function under consideration
00015  *      const double tol                Acceptable tolerance for the root
00016  *                                      position. It is an optional parameter
00017  *                                      with default value DBL_EPSILON
00018  *
00019  *      Zeroin returns an approximate location for the root with accuracy
00020  *      4*DBL_EPSILON*abs(x) + tol
00021  *
00022  * Algorithm
00023  *      G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
00024  *      computations. M., Mir, 1980, p.180 of the Russian edition
00025  *
00026  * The function makes use of a bissection procedure combined with
00027  * a linear or quadratic inverse interpolation.
00028  * At each step the code operates three abscissae - a, b, and c:
00029  *      b - the last and the best approximation to the root
00030  *      a - the last but one approximation
00031  *      c - the last but one or even an earlier approximation such that
00032  *              1) |f(b)| <= |f(c)|
00033  *              2) f(b) and f(c) have opposite signs, i.e. b and c encompass
00034  *                 the root
00035  * Given these abscissae, the code computes two new approximations, one by the 
00036  * bissection procedure and the other one from interpolation (if a,b, and c
00037  * are all different the quadratic interpolation is used, linear otherwise).
00038  * If the approximation obtained by the interpolation looks
00039  * reasonable (i.e. falls within the current interval [b,c], not too close
00040  * to the end points of the interval), the point is accepted as a new
00041  * approximation to the root. Otherwise, the result of the bissection is used.
00042  * Therefore, the range of uncertainty is guaranteed to tighten at 
00043  * least by a factor of 1.6
00044  *
00045  * $Id: zeroin.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00046  *
00047  ************************************************************************
00048  */
00049 
00050 #ifdef __GNUC__
00051 #pragma implementation "math_num.h"
00052 #endif
00053 #include "math_num.h"
00054 
00055 namespace linalg 
00056 {
00057     using namespace linalg;
00058     
00059 
00060 double zeroin(                          // An estimate to the root
00061         const double ax,                // Specify the interval the root
00062         const double bx,                // to be sought in
00063         UnivariateFunctor& f,           // Function under investigation
00064         const double tol)               // Acceptable tolerance
00065 {
00066   assure( tol > 0, "Tolerance must be positive");
00067   assure( bx > ax, 
00068          "Left end point of the interval should be strictly less than the "
00069          "right one" );
00070 
00071   double b = bx;                // the last and the best approx to the root
00072   double fb = f(b);
00073   double a = ax;                // the last but one approximation
00074   double fa = f(a);
00075   double c = a;                 // the last but one or even an earlier approx
00076   double fc = fa;               // see the condition above
00077 
00078   for(;;)               // Main iteration loop
00079   {
00080     const double prev_step = b-a;       // Step from the previous iteration
00081    
00082     if( fabs(fc) < fabs(fb) )
00083     {                                   // Swap data so that b would be the
00084       a = b;  b = c;  c = a;            // best approximation found so far
00085       fa=fb;  fb=fc;  fc=fa;
00086     }
00087                                         // Estimate the effective tolerance
00088     const double tol_act = 2*DBL_EPSILON*fabs(b) + tol/2;
00089     double new_step = (c-b)/2;          // Bissection step for this iteration
00090 
00091     if( fabs(new_step) <= tol_act || fb == 0 )
00092       return b;                         // Acceptable approximation is found
00093 
00094                         // Figuring out if the interpolation can be tried
00095     if( fabs(prev_step) >= tol_act      // If prev_step was large enough
00096         && fabs(fa) > fabs(fb) )                // and was in true direction,
00097     {                                   // Interpolatiom may be tried
00098 
00099       double p;                         // Interpolation step is calcu-
00100       double q;                         // lated in the form p/q; divi-
00101                                         // sion operations is delayed
00102                                         // until the last moment
00103       const double cb = c-b;
00104 
00105       if( a==c )                        // If we've got only two distinct
00106       {                                 // points linear interpolation
00107         register const double t1 = fb/fa;       // can only be applied
00108         p = cb*t1;
00109         q = 1.0 - t1;
00110       }
00111       else                              // Quadratic inverse interpolation
00112       {
00113         register const double t1=fb/fc, t2=fb/fa;
00114         q = fa/fc;
00115         p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
00116         q = (q-1.0) * (t1-1.0) * (t2-1.0);
00117       }
00118 
00119       if( p > 0 )                       // Formulas above computed new_step
00120         q = -q;                         // = p/q with the wrong sign (on purpose).
00121       else                              // Correct this, but in such a way so
00122         p = -p;                         // that p would be positive
00123       
00124       if( 2*p < (1.5*cb*q-fabs(tol_act*q))      // If b+p/q falls in [b,c]
00125          && 2*p < fabs(prev_step*q) )   // and isn't too large
00126         new_step = p/q;                 // it is accepted
00127                                         // If p/q is too large then the
00128                                         // bissection procedure can
00129                                         // reduce [b,c] to a larger
00130                                         // extent
00131     }
00132 
00133     if( fabs(new_step) < tol_act )      // Adjust the step to be not less
00134       new_step =  new_step > 0 ?        // than the tolerance
00135          tol_act : -tol_act;
00136 
00137     a = b;  fa = fb;                    // Save the previous approximation
00138     b += new_step;  fb = f(b);          // Do step to a new approximation
00139     if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) )
00140     {                                   // Adjust c for it to have the sign
00141       c = a;  fc = fa;                  // opposite to that of b
00142     }
00143   }
00144 
00145 }
00146 }

Generated on Wed Dec 15 21:20:38 2004 for vuVolume by  doxygen 1.3.9.1