Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:21

0001 // @(#)root/mathcore:$Id$
0002 // Authors: David Gonzalez Maline    01/2008
0003 
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2006 , LCG ROOT MathLib Team                         *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010 
0011 // Header file for RichardsonDerivator
0012 //
0013 // Created by: David Gonzalez Maline  : Mon Feb 4 2008
0014 //
0015 
0016 #ifndef ROOT_Math_RichardsonDerivator
0017 #define ROOT_Math_RichardsonDerivator
0018 
0019 #include <Math/IFunction.h>
0020 
0021 /**
0022     @defgroup Deriv Numerical Differentiation
0023     Classes for numerical differentiation
0024     @ingroup NumAlgo
0025 */
0026 
0027 
0028 namespace ROOT {
0029 namespace Math {
0030 
0031 //___________________________________________________________________________________________
0032 /**
0033    User class for calculating the derivatives of a function. It can calculate first (method Derivative1),
0034    second (method Derivative2) and third (method Derivative3) of a function.
0035 
0036    It uses the Richardson extrapolation method for function derivation in a given interval.
0037    The method use 2 derivative estimates (one computed with step h and one computed with step h/2)
0038    to compute a third, more accurate estimation. It is equivalent to the
0039    <a href = http://en.wikipedia.org/wiki/Five-point_stencil>5-point method</a>,
0040    which can be obtained with a Taylor expansion.
0041    A step size should be given, depending on x and f(x).
0042    An optimal step size value minimizes the truncation error of the expansion and the rounding
0043    error in evaluating x+h and f(x+h). A too small h will yield a too large rounding error while a too large
0044    h will give a large truncation error in the derivative approximation.
0045    A good discussion can be found in discussed in
0046    <a href=http://www.nrbook.com/a/bookcpdf/c5-7.pdf>Chapter 5.7</a>  of Numerical Recipes in C.
0047    By default a value of 0.001 is uses, acceptable in many cases.
0048 
0049    This class is implemented using code previously in  TF1::Derivate{,2,3}(). Now TF1 uses this class.
0050 
0051    @ingroup Deriv
0052 
0053  */
0054 
0055 class RichardsonDerivator {
0056 public:
0057 
0058    /** Destructor: Removes function if needed. */
0059    ~RichardsonDerivator();
0060 
0061    /** Default Constructor.
0062        Give optionally the step size for derivation. By default is 0.001, which is fine for x ~ 1
0063        Increase if x is in average larger or decrease if x is smaller
0064     */
0065    RichardsonDerivator(double h = 0.001);
0066 
0067    /** Construct from function and step size
0068     */
0069    RichardsonDerivator(const ROOT::Math::IGenFunction & f, double h = 0.001, bool copyFunc = false);
0070 
0071    /**
0072       Copy constructor
0073     */
0074    RichardsonDerivator(const RichardsonDerivator & rhs);
0075 
0076    /**
0077       Assignment operator
0078     */
0079    RichardsonDerivator & operator= ( const RichardsonDerivator & rhs);
0080 
0081 
0082    /** Returns the estimate of the absolute Error of the last derivative calculation. */
0083    double Error () const {  return fLastError; }
0084 
0085 
0086    /**
0087       Returns the first derivative of the function at point x,
0088       computed by Richardson's extrapolation method (use 2 derivative estimates
0089       to compute a third, more accurate estimation)
0090       first, derivatives with steps h and h/2 are computed by central difference formulas
0091      \f[
0092       D(h) = \frac{f(x+h) - f(x-h)}{2h}
0093      \f]
0094       the final estimate
0095      \f[
0096       D = \frac{4D(h/2) - D(h)}{3}
0097      \f]
0098        "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
0099 
0100       the argument eps may be specified to control the step size (precision).
0101       the step size is taken as eps*(xmax-xmin).
0102       the default value (0.001) should be good enough for the vast majority
0103       of functions. Give a smaller value if your function has many changes
0104       of the second derivative in the function range.
0105 
0106       Getting the error via TF1::DerivativeError:
0107         (total error = roundoff error + interpolation error)
0108       the estimate of the roundoff error is taken as follows:
0109      \f[
0110          err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}},
0111      \f]
0112       where k is the double precision, ai are coefficients used in
0113       central difference formulas
0114       interpolation error is decreased by making the step size h smaller.
0115    */
0116    double Derivative1 (double x) { return Derivative1(*fFunction,x,fStepSize); }
0117    double operator() (double x) { return Derivative1(*fFunction,x,fStepSize); }
0118 
0119    /**
0120       First Derivative calculation passing function object and step-size
0121     */
0122    double Derivative1(const IGenFunction & f, double x, double h);
0123 
0124    /// Computation of the first derivative using a forward formula
0125    double DerivativeForward(double x) {
0126       return DerivativeForward(*fFunction, x, fStepSize);
0127    }
0128 
0129    /// Computation of the first derivative using a forward formula
0130    double DerivativeForward(const IGenFunction &f, double x, double h);
0131 
0132       /// Computation of the first derivative using a backward formula
0133    double DerivativeBackward(double x) {
0134       return DerivativeBackward(*fFunction, x, fStepSize);
0135    }
0136 
0137    /// Computation of the first derivative using a forward formula
0138    double DerivativeBackward(const IGenFunction &f, double x, double h) {
0139       return DerivativeForward(f, x, -h);
0140    }
0141 
0142    /**
0143       Returns the second derivative of the function at point x,
0144       computed by Richardson's extrapolation method (use 2 derivative estimates
0145       to compute a third, more accurate estimation)
0146       first, derivatives with steps h and h/2 are computed by central difference formulas
0147      \f[
0148          D(h) = \frac{f(x+h) - 2f(x) + f(x-h)}{h^{2}}
0149      \f]
0150       the final estimate
0151      \f[
0152          D = \frac{4D(h/2) - D(h)}{3}
0153      \f]
0154        "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
0155 
0156       the argument eps may be specified to control the step size (precision).
0157       the step size is taken as eps*(xmax-xmin).
0158       the default value (0.001) should be good enough for the vast majority
0159       of functions. Give a smaller value if your function has many changes
0160       of the second derivative in the function range.
0161 
0162       Getting the error via TF1::DerivativeError:
0163         (total error = roundoff error + interpolation error)
0164       the estimate of the roundoff error is taken as follows:
0165      \f[
0166          err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}},
0167      \f]
0168       where k is the double precision, ai are coefficients used in
0169       central difference formulas
0170       interpolation error is decreased by making the step size h smaller.
0171    */
0172    double Derivative2 (double x) {
0173       return Derivative2( *fFunction, x, fStepSize);
0174    }
0175 
0176    /**
0177       Second Derivative calculation passing function and step-size
0178     */
0179    double Derivative2(const IGenFunction & f, double x, double h);
0180 
0181    /**
0182       Returns the third derivative of the function at point x,
0183       computed by Richardson's extrapolation method (use 2 derivative estimates
0184       to compute a third, more accurate estimation)
0185       first, derivatives with steps h and h/2 are computed by central difference formulas
0186      \f[
0187          D(h) = \frac{f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2h)}{2h^{3}}
0188      \f]
0189       the final estimate
0190      \f[
0191          D = \frac{4D(h/2) - D(h)}{3}
0192      \f]
0193        "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
0194 
0195       the argument eps may be specified to control the step size (precision).
0196       the step size is taken as eps*(xmax-xmin).
0197       the default value (0.001) should be good enough for the vast majority
0198       of functions. Give a smaller value if your function has many changes
0199       of the second derivative in the function range.
0200 
0201       Getting the error via TF1::DerivativeError:
0202         (total error = roundoff error + interpolation error)
0203       the estimate of the roundoff error is taken as follows:
0204      \f[
0205          err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}},
0206      \f]
0207       where k is the double precision, ai are coefficients used in
0208       central difference formulas
0209       interpolation error is decreased by making the step size h smaller.
0210    */
0211    double Derivative3 (double x) {
0212       return Derivative3(*fFunction, x, fStepSize);
0213    }
0214 
0215    /**
0216       Third Derivative calculation passing function and step-size
0217     */
0218    double Derivative3(const IGenFunction & f, double x, double h);
0219 
0220    /** Set function for derivative calculation (copy the function if option has been enabled in the constructor)
0221 
0222        \@param f Function to be differentiated
0223    */
0224    void SetFunction (const IGenFunction & f);
0225 
0226    /** Set step size for derivative calculation
0227 
0228        \@param h step size for calculation
0229    */
0230    void SetStepSize (double h) { fStepSize = h; }
0231 
0232 protected:
0233 
0234    bool fFunctionCopied;     ///< flag to control if function is copied in the class
0235    double fStepSize;         ///< step size used for derivative calculation
0236    double fLastError;        ///<  error estimate of last derivative calculation
0237    const IGenFunction* fFunction;  ///< pointer to function
0238 
0239 };
0240 
0241 } // end namespace Math
0242 
0243 } // end namespace ROOT
0244 
0245 #endif /* ROOT_Math_RichardsonDerivator */
0246