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          <li>RootFinder::kMODAB
0059          </ol>
0060          <li>Root Finding Algorithms using Derivatives
0061          <ol>
0062          <li>RootFinder::kGSL_NEWTON
0063          <li>RootFinder::kGSL_SECANT
0064          <li>RootFinder::kGSL_STEFFENSON
0065          </ol>
0066          </ul>
0067
0068          This class does not cupport copying
0069
0070          @ingroup RootFinders
0071
0072       */
0073
0074       class RootFinder {
0075
0076       public:
0077
0078          enum EType {kBRENT,                                     // Methods from MathCore
0079                      kGSL_BISECTION, kGSL_FALSE_POS, kGSL_BRENT, // GSL Normal
0080                      kGSL_NEWTON, kGSL_SECANT, kGSL_STEFFENSON,  // GSL Derivatives
0081                      kMODAB                                      // Modified A&B method added in MathCore
0082          };
0083
0084          /**
0085             Construct a Root-Finder algorithm
0086          */
0087          RootFinder(RootFinder::EType type = RootFinder::kBRENT);
0088          virtual ~RootFinder();
0089
0090          // usually copying is non trivial, so we delete this
0091          RootFinder(const RootFinder & ) = delete;
0092          RootFinder & operator = (const RootFinder & rhs) = delete;
0093          RootFinder(RootFinder && ) = delete;
0094          RootFinder & operator = (RootFinder && rhs) = delete;
0095
0096          bool SetMethod(RootFinder::EType type = RootFinder::kBRENT);
0097
0098          /**
0099             Provide to the solver the function and the initial search interval [xlow, xup]
0100             for algorithms not using derivatives (bracketing algorithms)
0101             The templated function f must be of a type implementing the \a operator() method,
0102             <em>  double  operator() (  double  x ) </em>
0103             Returns non zero if interval is not valid (i.e. does not contains a root)
0104          */
0105
0106          bool SetFunction( const IGenFunction & f, double xlow, double xup) {
0107             return fSolver->SetFunction( f, xlow, xup);
0108          }
0109
0110
0111          /**
0112             Provide to the solver the function and an initial estimate of the root,
0113             for algorithms using derivatives.
0114             The templated function f must be of a type implementing the \a operator()
0115             and the \a Gradient() methods.
0116             <em>  double  operator() (  double  x ) </em>
0117             Returns non zero if starting point is not valid
0118          */
0119
0120          bool  SetFunction( const IGradFunction & f, double xstart) {
0121             return fSolver->SetFunction( f, xstart);
0122          }
0123
0124          template<class Function, class Derivative>
0125          bool Solve(Function &f, Derivative &d, double start,
0126                    int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10);
0127
0128          template<class Function>
0129          bool Solve(Function &f, double min, double max,
0130                    int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10);
0131
0132          /**
0133              Compute the roots iterating until the estimate of the Root is within the required tolerance returning
0134              the iteration Status
0135          */
0136          bool Solve( int maxIter = 100, double absTol = 1E-8, double relTol = 1E-10) {
0137             return fSolver->Solve( maxIter, absTol, relTol );
0138          }
0139
0140          /**
0141              Return the number of iteration performed to find the Root.
0142          */
0143          int Iterations() const {
0144             return fSolver->Iterations();
0145          }
0146
0147          /**
0148             Perform a single iteration and return the Status
0149          */
0150          int Iterate() {
0151             return fSolver->Iterate();
0152          }
0153
0154          /**
0155             Return the current and latest estimate of the Root
0156          */
0157          double Root() const {
0158             return fSolver->Root();
0159          }
0160
0161          /**
0162             Return the status of the last estimate of the Root
0163             = 0 OK, not zero failure
0164          */
0165          int Status() const {
0166             return fSolver->Status();
0167          }
0168
0169
0170          /**
0171             Return the current and latest estimate of the lower value of the Root-finding interval (for bracketing algorithms)
0172          */
0173 /*   double XLower() const {  */
0174 /*     return fSolver->XLower();  */
0175 /*   } */
0176
0177          /**
0178             Return the current and latest estimate of the upper value of the Root-finding interval (for bracketing algorithms)
0179          */
0180 /*   double XUpper() const {  */
0181 /*     return  fSolver->XUpper();  */
0182 /*   } */
0183
0184          /**
0185             Get Name of the Root-finding solver algorithm
0186          */
0187          const char * Name() const {
0188             return fSolver->Name();
0189          }
0190
0191
0192       protected:
0193
0194
0195       private:
0196
0197          IRootFinderMethod* fSolver;   // type of algorithm to be used
0198
0199
0200       };
0201
0202    } // namespace Math
0203 } // namespace ROOT
0204
0205
0206 #include "Math/WrappedFunction.h"
0207
0208 #include "Math/Functor.h"
0209
0210 /**
0211  * Solve `f(x) = 0`, given a derivative `d`.
0212  * @param f Function whose root should be found.
0213  * @param d Derivative of the function.
0214  * @param start Starting point for iteration.
0215  * @param maxIter Maximum number of iterations, passed to Solve(int,double,double)
0216  * @param absTol Absolute tolerance, as in Solve(int,double,double)
0217  * @param relTol Relative tolerance, passed to Solve(int,double,double)
0218  * @return true if a root was found. Retrieve the result using Root().
0219  */
0220 template<class Function, class Derivative>
0221 bool ROOT::Math::RootFinder::Solve(Function &f, Derivative &d, double start,
0222                                   int maxIter, double absTol, double relTol)
0223 {
0224    if (!fSolver) return false;
0225    ROOT::Math::GradFunctor1D wf(f, d);
0226    bool ret = fSolver->SetFunction(wf, start);
0227    if (!ret) return false;
0228    return Solve(maxIter, absTol, relTol);
0229 }
0230
0231 /**
0232  * Solve `f(x) = 0` numerically.
0233  * @param f Function whose root should be found.
0234  * @param min Minimum allowed value of `x`.
0235  * @param max Maximum allowed value of `x`.
0236  * @param maxIter Maximum number of iterations, passed to Solve(int,double,double)
0237  * @param absTol Absolute tolerance, as in Solve(int,double,double)
0238  * @param relTol Relative tolerance, passed to Solve(int,double,double)
0239  * @return true if a root was found. Retrieve the result using Root().
0240  */
0241 template<class Function>
0242 bool ROOT::Math::RootFinder::Solve(Function &f, double min, double max,
0243                                   int maxIter, double absTol, double relTol)
0244 {
0245    if (!fSolver) return false;
0246    ROOT::Math::WrappedFunction<Function &> wf(f);
0247    bool ret = fSolver->SetFunction(wf, min, max);
0248    if (!ret) return false;
0249    return Solve(maxIter, absTol, relTol);
0250 }
0251
0252 #endif /* ROOT_Math_RootFinder */