// @(#)root/mathcore:$Id$ // Authors: David Gonzalez Maline 01/2008 /********************************************************************** * * * Copyright (c) 2006 , LCG ROOT MathLib Team * * * * * **********************************************************************/ // Header file for RichardsonDerivator // // Created by: David Gonzalez Maline : Mon Feb 4 2008 // #ifndef ROOT_Math_RichardsonDerivator #define ROOT_Math_RichardsonDerivator #include /** @defgroup Deriv Numerical Differentiation Classes for numerical differentiation @ingroup NumAlgo */ namespace ROOT { namespace Math { //___________________________________________________________________________________________ /** User class for calculating the derivatives of a function. It can calculate first (method Derivative1), second (method Derivative2) and third (method Derivative3) of a function. It uses the Richardson extrapolation method for function derivation in a given interval. The method use 2 derivative estimates (one computed with step h and one computed with step h/2) to compute a third, more accurate estimation. It is equivalent to the 5-point method, which can be obtained with a Taylor expansion. A step size should be given, depending on x and f(x). An optimal step size value minimizes the truncation error of the expansion and the rounding error in evaluating x+h and f(x+h). A too small h will yield a too large rounding error while a too large h will give a large truncation error in the derivative approximation. A good discussion can be found in discussed in Chapter 5.7 of Numerical Recipes in C. By default a value of 0.001 is uses, acceptable in many cases. This class is implemented using code previously in TF1::Derivate{,2,3}(). Now TF1 uses this class. @ingroup Deriv */ class RichardsonDerivator { public: /** Destructor: Removes function if needed. */ ~RichardsonDerivator(); /** Default Constructor. Give optionally the step size for derivation. By default is 0.001, which is fine for x ~ 1 Increase if x is in average larger or decrease if x is smaller */ RichardsonDerivator(double h = 0.001); /** Construct from function and step size */ RichardsonDerivator(const ROOT::Math::IGenFunction & f, double h = 0.001, bool copyFunc = false); /** Copy constructor */ RichardsonDerivator(const RichardsonDerivator & rhs); /** Assignment operator */ RichardsonDerivator & operator= ( const RichardsonDerivator & rhs); /** Returns the estimate of the absolute Error of the last derivative calculation. */ double Error () const { return fLastError; } /** Returns the first derivative of the function at point x, computed by Richardson's extrapolation method (use 2 derivative estimates to compute a third, more accurate estimation) first, derivatives with steps h and h/2 are computed by central difference formulas \f[ D(h) = \frac{f(x+h) - f(x-h)}{2h} \f] the final estimate \f[ D = \frac{4D(h/2) - D(h)}{3} \f] "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition" the argument eps may be specified to control the step size (precision). the step size is taken as eps*(xmax-xmin). the default value (0.001) should be good enough for the vast majority of functions. Give a smaller value if your function has many changes of the second derivative in the function range. Getting the error via TF1::DerivativeError: (total error = roundoff error + interpolation error) the estimate of the roundoff error is taken as follows: \f[ err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}}, \f] where k is the double precision, ai are coefficients used in central difference formulas interpolation error is decreased by making the step size h smaller. */ double Derivative1 (double x) { return Derivative1(*fFunction,x,fStepSize); } double operator() (double x) { return Derivative1(*fFunction,x,fStepSize); } /** First Derivative calculation passing function object and step-size */ double Derivative1(const IGenFunction & f, double x, double h); /// Computation of the first derivative using a forward formula double DerivativeForward(double x) { return DerivativeForward(*fFunction, x, fStepSize); } /// Computation of the first derivative using a forward formula double DerivativeForward(const IGenFunction &f, double x, double h); /// Computation of the first derivative using a backward formula double DerivativeBackward(double x) { return DerivativeBackward(*fFunction, x, fStepSize); } /// Computation of the first derivative using a forward formula double DerivativeBackward(const IGenFunction &f, double x, double h) { return DerivativeForward(f, x, -h); } /** Returns the second derivative of the function at point x, computed by Richardson's extrapolation method (use 2 derivative estimates to compute a third, more accurate estimation) first, derivatives with steps h and h/2 are computed by central difference formulas \f[ D(h) = \frac{f(x+h) - 2f(x) + f(x-h)}{h^{2}} \f] the final estimate \f[ D = \frac{4D(h/2) - D(h)}{3} \f] "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition" the argument eps may be specified to control the step size (precision). the step size is taken as eps*(xmax-xmin). the default value (0.001) should be good enough for the vast majority of functions. Give a smaller value if your function has many changes of the second derivative in the function range. Getting the error via TF1::DerivativeError: (total error = roundoff error + interpolation error) the estimate of the roundoff error is taken as follows: \f[ err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}}, \f] where k is the double precision, ai are coefficients used in central difference formulas interpolation error is decreased by making the step size h smaller. */ double Derivative2 (double x) { return Derivative2( *fFunction, x, fStepSize); } /** Second Derivative calculation passing function and step-size */ double Derivative2(const IGenFunction & f, double x, double h); /** Returns the third derivative of the function at point x, computed by Richardson's extrapolation method (use 2 derivative estimates to compute a third, more accurate estimation) first, derivatives with steps h and h/2 are computed by central difference formulas \f[ D(h) = \frac{f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2h)}{2h^{3}} \f] the final estimate \f[ D = \frac{4D(h/2) - D(h)}{3} \f] "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition" the argument eps may be specified to control the step size (precision). the step size is taken as eps*(xmax-xmin). the default value (0.001) should be good enough for the vast majority of functions. Give a smaller value if your function has many changes of the second derivative in the function range. Getting the error via TF1::DerivativeError: (total error = roundoff error + interpolation error) the estimate of the roundoff error is taken as follows: \f[ err = k\sqrt{f(x)^{2} + x^{2}deriv^{2}}\sqrt{\sum ai^{2}}, \f] where k is the double precision, ai are coefficients used in central difference formulas interpolation error is decreased by making the step size h smaller. */ double Derivative3 (double x) { return Derivative3(*fFunction, x, fStepSize); } /** Third Derivative calculation passing function and step-size */ double Derivative3(const IGenFunction & f, double x, double h); /** Set function for derivative calculation (copy the function if option has been enabled in the constructor) \@param f Function to be differentiated */ void SetFunction (const IGenFunction & f); /** Set step size for derivative calculation \@param h step size for calculation */ void SetStepSize (double h) { fStepSize = h; } protected: bool fFunctionCopied; ///< flag to control if function is copied in the class double fStepSize; ///< step size used for derivative calculation double fLastError; ///< error estimate of last derivative calculation const IGenFunction* fFunction; ///< pointer to function }; } // end namespace Math } // end namespace ROOT #endif /* ROOT_Math_RichardsonDerivator */