root/Math/RootFinder.h

0001 // @(#)root/mathmore:$Id$
0002 // Authors: L. Moneta, A. Zsenei   08/2005
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2004  CERN                                           *
0007   * All rights reserved.                                               *
0008   *                                                                    *
0009   * For the licensing terms see $ROOTSYS/LICENSE.                      *
0010   * For the list of contributors see $ROOTSYS/README/CREDITS.          *
0011   *                                                                    *
0012   **********************************************************************/
0013
0014 // Header file for class RootFinder
0015 //
0016 // Created by: moneta  at Sun Nov 14 16:59:55 2004
0017 //
0018 // Last update: Sun Nov 14 16:59:55 2004
0019 //
0020 #ifndef ROOT_Math_RootFinder
0021 #define ROOT_Math_RootFinder
0022
0023
0024 #include "Math/IFunctionfwd.h"
0025
0026 #include "Math/IRootFinderMethod.h"
0027
0028
0029 /**
0030    @defgroup RootFinders One-dimensional Root-Finding
0031    Classes implementing algorithms for finding the roots of a one-dimensional function.
0032    Various implementations exist in MathCore and MathMore
0033    The user interacts with a proxy class ROOT::Math::RootFinder which creates behind
0034    the chosen algorithms which are implemented using the ROOT::Math::IRootFinderMethod interface
0035
0036    @ingroup NumAlgo
0037 */
0038
0039
0040 namespace ROOT {
0041    namespace Math {
0042
0043
0044 //_____________________________________________________________________________________
0045       /**
0046          User Class to find the Root of one dimensional functions.
0047          The GSL Methods are implemented in MathMore and they are loaded automatically
0048          via the plug-in manager
0049
0050          The possible types of Root-finding algorithms are:
0051          <ul>
0052          <li>Root Bracketing Algorithms which do not require function derivatives
0053          <ol>
0054          <li>RootFinder::kBRENT  (default method implemented in MathCore)
0055          <li>RootFinder::kGSL_BISECTION
0056          <li>RootFinder::kGSL_FALSE_POS
0057          <li>RootFinder::kGSL_BRENT
0058          </ol>
0059          <li>Root Finding Algorithms using Derivatives
0060          <ol>
0061          <li>RootFinder::kGSL_NEWTON
0062          <li>RootFinder::kGSL_SECANT
0063          <li>RootFinder::kGSL_STEFFENSON
0064          </ol>
0065          </ul>
0066
0067          This class does not cupport copying
0068
0069          @ingroup RootFinders
0070
0071       */
0072
0073       class RootFinder {
0074
0075       public:
0076
0077          enum EType { kBRENT,                                   // Methods from MathCore
0078                      kGSL_BISECTION, kGSL_FALSE_POS, kGSL_BRENT, // GSL Normal
0079                      kGSL_NEWTON, kGSL_SECANT, kGSL_STEFFENSON   // GSL Derivatives
0080          };
0081
0082          /**
0083             Construct a Root-Finder algorithm
0084          */
0085          RootFinder(RootFinder::EType type = RootFinder::kBRENT);
0086          virtual ~RootFinder();
0087
0088       private:
0089          // usually copying is non trivial, so we make this unaccessible
0090          RootFinder(const RootFinder & ) {}
0091          RootFinder & operator = (const RootFinder & rhs)
0092          {
0093             if (this == &rhs) return *this;  // time saving self-test
0094             return *this;
0095          }
0096
0097       public:
0098
0099          bool SetMethod(RootFinder::EType type = RootFinder::kBRENT);
0100
0101          /**
0102             Provide to the solver the function and the initial search interval [xlow, xup]
0103             for algorithms not using derivatives (bracketing algorithms)
0104             The templated function f must be of a type implementing the \a operator() method,
0105             <em>  double  operator() (  double  x ) </em>
0106             Returns non zero if interval is not valid (i.e. does not contains a root)
0107          */
0108
0109          bool SetFunction( const IGenFunction & f, double xlow, double xup) {
0110             return fSolver->SetFunction( f, xlow, xup);
0111          }
0112
0113
0114          /**
0115             Provide to the solver the function and an initial estimate of the root,
0116             for algorithms using derivatives.
0117             The templated function f must be of a type implementing the \a operator()
0118             and the \a Gradient() methods.
0119             <em>  double  operator() (  double  x ) </em>
0120             Returns non zero if starting point is not valid
0121          */
0122
0123          bool  SetFunction( const IGradFunction & f, double xstart) {
0124             return fSolver->SetFunction( f, xstart);
0125          }
0126
0127          template<class Function, class Derivative>
0128          bool Solve(Function &f, Derivative &d, double start,
0129                    int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10);
0130
0131          template<class Function>
0132          bool Solve(Function &f, double min, double max,
0133                    int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10);
0134
0135          /**
0136              Compute the roots iterating until the estimate of the Root is within the required tolerance returning
0137              the iteration Status
0138          */
0139          bool Solve( int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10) {
0140             return fSolver->Solve( maxIter, absTol, relTol );
0141          }
0142
0143          /**
0144              Return the number of iteration performed to find the Root.
0145          */
0146          int Iterations() const {
0147             return fSolver->Iterations();
0148          }
0149
0150          /**
0151             Perform a single iteration and return the Status
0152          */
0153          int Iterate() {
0154             return fSolver->Iterate();
0155          }
0156
0157          /**
0158             Return the current and latest estimate of the Root
0159          */
0160          double Root() const {
0161             return fSolver->Root();
0162          }
0163
0164          /**
0165             Return the status of the last estimate of the Root
0166             = 0 OK, not zero failure
0167          */
0168          int Status() const {
0169             return fSolver->Status();
0170          }
0171
0172
0173          /**
0174             Return the current and latest estimate of the lower value of the Root-finding interval (for bracketing algorithms)
0175          */
0176 /*   double XLower() const {  */
0177 /*     return fSolver->XLower();  */
0178 /*   } */
0179
0180          /**
0181             Return the current and latest estimate of the upper value of the Root-finding interval (for bracketing algorithms)
0182          */
0183 /*   double XUpper() const {  */
0184 /*     return  fSolver->XUpper();  */
0185 /*   } */
0186
0187          /**
0188             Get Name of the Root-finding solver algorithm
0189          */
0190          const char * Name() const {
0191             return fSolver->Name();
0192          }
0193
0194
0195       protected:
0196
0197
0198       private:
0199
0200          IRootFinderMethod* fSolver;   // type of algorithm to be used
0201
0202
0203       };
0204
0205    } // namespace Math
0206 } // namespace ROOT
0207
0208
0209 #include "Math/WrappedFunction.h"
0210
0211 #include "Math/Functor.h"
0212
0213 /**
0214  * Solve `f(x) = 0`, given a derivative `d`.
0215  * @param f Function whose root should be found.
0216  * @param d Derivative of the function.
0217  * @param start Starting point for iteration.
0218  * @param maxIter Maximum number of iterations, passed to Solve(int,double,double)
0219  * @param absTol Absolute tolerance, as in Solve(int,double,double)
0220  * @param relTol Relative tolerance, passed to Solve(int,double,double)
0221  * @return true if a root was found. Retrieve the result using Root().
0222  */
0223 template<class Function, class Derivative>
0224 bool ROOT::Math::RootFinder::Solve(Function &f, Derivative &d, double start,
0225                                   int maxIter, double absTol, double relTol)
0226 {
0227    if (!fSolver) return false;
0228    ROOT::Math::GradFunctor1D wf(f, d);
0229    bool ret = fSolver->SetFunction(wf, start);
0230    if (!ret) return false;
0231    return Solve(maxIter, absTol, relTol);
0232 }
0233
0234 /**
0235  * Solve `f(x) = 0` numerically.
0236  * @param f Function whose root should be found.
0237  * @param min Minimum allowed value of `x`.
0238  * @param max Maximum allowed value of `x`.
0239  * @param maxIter Maximum number of iterations, passed to Solve(int,double,double)
0240  * @param absTol Absolute tolerance, as in Solve(int,double,double)
0241  * @param relTol Relative tolerance, passed to Solve(int,double,double)
0242  * @return true if a root was found. Retrieve the result using Root().
0243  */
0244 template<class Function>
0245 bool ROOT::Math::RootFinder::Solve(Function &f, double min, double max,
0246                                   int maxIter, double absTol, double relTol)
0247 {
0248    if (!fSolver) return false;
0249    ROOT::Math::WrappedFunction<Function &> wf(f);
0250    bool ret = fSolver->SetFunction(wf, min, max);
0251    if (!ret) return false;
0252    return Solve(maxIter, absTol, relTol);
0253 }
0254
0255 #endif /* ROOT_Math_RootFinder */