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