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

fminbr.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 one-dimensional minimizer 
00008  *
00009  *           finds a local minimum of a single argument function
00010  *                        over a given range
00011  *
00012  * Input
00013  *      double fminbr(ax,bx,f,tol)
00014  *      const double ax                 a and b, a < b, specify the interval
00015  *      const double bx                 the minimum is to be sought in
00016  *      UnivariateFunctor& f            The function under consideration
00017  *      const double tol                Acceptable tolerance for the minimum
00018  *                                      location. It is an optional parameter
00019  *                                      with default value DBL_EPSILON
00020  *
00021  * Output
00022  *      Fminbr returns an estimate to the location of the minimum
00023  *      with accuracy 3*SQRT_EPSILON*abs(x) + tol.
00024  *      The procedure can only determine a local minimum, which coincides with
00025  *      the global one if and only if the function under investigation is
00026  *      unimodular.
00027  *      If a function being examined possesses no local minimum within
00028  *      the given interval, Fminbr returns either the left or the right end
00029  *      point of the interval, wherever the function value is smaller.
00030  *
00031  * Algorithm
00032  *      G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
00033  *      computations. M., Mir, 1980, p.202 of the Russian edition
00034  *
00035  * The function makes use of a "gold section" procedure combined with
00036  * a parabolic interpolation.
00037  * At each step the code operates three abscissae - x,v, and w.
00038  *      x - the last and the best approximation to the minimum location,
00039  *              i.e. f(x) <= f(a) or/and f(x) <= f(b)
00040  *          (if the function f has a local minimum in (a,b), then both
00041  *           conditions are met after one or two steps).
00042  *      v,w are previous approximations to the location of the minimum.
00043  *      They may coincide with a, b, or x (although the algorithm tries
00044  *      to make all u, v, and w distinct). 
00045  * Points x, v, and w are used to construct an interpolating parabola,
00046  * whose minimum is regarded as a new approximation to the minimum
00047  * of the function, provided the parabola's minimum falls within [a,b]
00048  * and reduces the current interval [a,b] to a larger extent than the
00049  * gold section procedure does.
00050  * When f(x) has a positive second derivative at the point of minimum
00051  * (which does not coincide with a or b) the procedure converges
00052  * superlinearly at a rate of about 1.324
00053  *
00054  * $Id: fminbr.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00055  *
00056  ************************************************************************
00057  */
00058 
00059 #include "math_num.h"
00060 
00061 namespace linalg 
00062 {
00063     using namespace linalg;
00064     
00065 
00066 double fminbr(                          // An estimate to the min location
00067         const double ax,                // Specify the interval the minimum
00068         const double bx,                // to be sought in
00069         UnivariateFunctor& f,           // Function under investigation
00070         const double tol)               // Acceptable tolerance
00071 {
00072   assure( tol > 0, "Tolerance must be positive");
00073   assure( bx > ax, 
00074          "Left end point of the interval should be strictly less than the "
00075          "right one" );
00076   
00077   static const double r = (3-sqrt(5.0))/2;      // The golden section ratio
00078   static const double sqrt_eps = sqrt(DBL_EPSILON);
00079   
00080   double a = ax, b = bx;                // Current interval
00081   double v = a + r*(b-a);               // First step - always gold section
00082   double fv = f(v);
00083   double x = v;                         // the last and the best approximation
00084   double fx = fv;
00085   double w = v;                         // a previous approx to the min
00086   double fw = fv;
00087 
00088   for(;;)               // Main iteration loop
00089   {
00090     const double range = b-a;           // Interval where the minimum
00091                                         // is searched in
00092     const double midpoint = (a+b)/2;
00093     const double tol_act =              // The effective tolerance
00094                 sqrt_eps*fabs(x) + tol/3;
00095 
00096        
00097 
00098     if( 2*fabs(x-midpoint) + range <= 4*tol_act )
00099       return x;                         // Acceptable approximation is found
00100 
00101                                         // Compute a new step with the gold
00102                                         // section
00103     double new_step = r * ( x < midpoint ? b-x : a-x );
00104 
00105 
00106                         // Decide on the interpolation  
00107     if( fabs(x-w) >= tol_act  )         // If x and w are distinct
00108     {                                   // interpolatiom may be tried
00109       register double p;                // Interpolation step is calcula-
00110       register double q;                // ted as p/q; division operation
00111                                         // is delayed until last moment
00112       register double t;
00113 
00114       t = (x-w) * (fx-fv);
00115       q = (x-v) * (fx-fw);
00116       p = (x-v)*q - (x-w)*t;
00117       q = 2*(q-t);
00118 
00119       if( q > 0 )                       // Formulas above computed new_step
00120         p = -p;                         // = p/q with a wrong sign (on purpose).
00121       else                              // Correct this, but in such a way so
00122         q = -q;                         // that q would be positive
00123 
00124       if( fabs(p) < fabs(new_step*q) && // If x+p/q falls in [a,b] and is not
00125          p > q*(a-x+2*tol_act) &&       // too close to a and b, and isn't
00126          p < q*(b-x-2*tol_act)  )       // too large, it is accepted
00127            new_step = p/q;
00128                                         // If p/q is too large then the
00129                                         // gold section procedure would
00130                                         // reduce [a,b] to larger extent
00131     }
00132 
00133     if( fabs(new_step) < tol_act )      // Adjust the step to be not less
00134       new_step =  new_step > 0 ?        // than tolerance
00135         tol_act : -tol_act;
00136 
00137                                 // Obtain the next approximation to min
00138                                 // and reduce the encompassing interval
00139     register const double t = x + new_step;  // Tentative point for the min
00140     register const double ft = f(t);
00141     if( ft <= fx )
00142     {                                   // t is a better approximation
00143       ( t < x ? b : a ) = x;            // Reduce the interval so that
00144                                         // t would fall within it
00145       v = w;  w = x;  x = t;            // Assign the best approx to x
00146       fv=fw;  fw=fx;  fx=ft;
00147     }
00148     else                                // x remains the better approx
00149     {
00150       ( t < x ? a : b ) = t;            // Reduce the interval encompassing x
00151       
00152       if( ft <= fw || w==x )
00153       {
00154         v = w;  w = t;
00155         fv=fw;  fw=ft;
00156       }
00157       else if( ft<=fv || v==x || v==w )
00158         v = t, fv = ft;
00159     }
00160   }             // ===== End of loop =====
00161 
00162 }
00163 }

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