root/Math/GaussIntegrator.h

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 GaussIntegrator
0012 //
0013 // Created by: David Gonzalez Maline  : Wed Jan 16 2008
0014 //
0015
0016 #ifndef ROOT_Math_GaussIntegrator
0017 #define ROOT_Math_GaussIntegrator
0018
0019 #include "Math/IFunction.h"
0020
0021 #include "Math/VirtualIntegrator.h"
0022
0023 #include <vector>
0024
0025 namespace ROOT {
0026 namespace Math {
0027
0028
0029 //___________________________________________________________________________________________
0030 /**
0031    User class for performing function integration.
0032
0033    It will use the Gauss Method for function integration in a given interval.
0034    This class is implemented from TF1::Integral().
0035
0036    @ingroup Integration
0037
0038  */
0039
0040 class GaussIntegrator: public VirtualIntegratorOneDim {
0041
0042
0043 public:
0044
0045    /** Destructor */
0046    ~GaussIntegrator() override;
0047
0048    /** Default Constructor.
0049        If the tolerance are not given, use default values specified in  ROOT::Math::IntegratorOneDimOptions
0050     */
0051    GaussIntegrator(double absTol = -1, double relTol = -1);
0052
0053
0054    /** Static function: set the fgAbsValue flag.
0055        By default TF1::Integral uses the original function value to compute the integral
0056        However, TF1::Moment, CentralMoment require to compute the integral
0057        using the absolute value of the function.
0058    */
0059    void AbsValue(bool flag);
0060
0061
0062    // Implementing VirtualIntegrator Interface
0063
0064    /** Set the desired relative Error. */
0065    void SetRelTolerance (double eps) override { fEpsRel = eps; }
0066
0067    /** This method is not implemented. */
0068    void SetAbsTolerance (double eps) override { fEpsAbs = eps; }
0069
0070    /** Returns the result of the last Integral calculation. */
0071    double Result () const override;
0072
0073    /** Return the estimate of the absolute Error of the last Integral calculation. */
0074    double Error () const override;
0075
0076    /** return the status of the last integration - 0 in case of success */
0077    int Status () const override;
0078
0079    // Implementing VirtualIntegratorOneDim Interface
0080
0081    /**
0082      Returns Integral of function between a and b.
0083      Based on original CERNLIB routine DGAUSS by Sigfried Kolbig
0084      converted to C++ by Rene Brun
0085
0086      This function computes, to an attempted specified accuracy, the value
0087      of the integral.
0088
0089     Method:
0090        For any interval [a,b] we define g8(a,b) and g16(a,b) to be the 8-point
0091        and 16-point Gaussian quadrature approximations to
0092    \f[
0093       I = \int^{b}_{a} f(x)dx
0094    \f]
0095       and define
0096    \f[
0097       r(a,b) = \frac{\left|g_{16}(a,b)-g_{8}(a,b)\right|}{1+\left|g_{16}(a,b)\right|}
0098    \f]
0099       Then,
0100    \f[
0101       G = \sum_{i=1}^{k}g_{16}(x_{i-1},x_{i})
0102    \f]
0103       where, starting with \f$x_{0} = A\f$ and finishing with \f$x_{k} = B\f$,
0104       the subdivision points \f$x_{i}(i=1,2,...)\f$ are given by
0105    \f[
0106       x_{i} = x_{i-1} + \lambda(B-x_{i-1})
0107    \f]
0108       \f$\lambda\f$ is equal to the first member of the
0109       sequence 1,1/2,1/4,... for which \f$r(x_{i-1}, x_{i}) < EPS\f$.
0110       If, at any stage in the process of subdivision, the ratio
0111   \f[
0112       q = \left|\frac{x_{i}-x_{i-1}}{B-A}\right|
0113   \f]
0114       is so small that 1+0.005q is indistinguishable from 1 to
0115       machine accuracy, an error exit occurs with the function value
0116       set equal to zero.
0117
0118    Accuracy:
0119       The user provides absolute and relative error bounds (epsrel and epsabs) and the
0120       algorithm will stop when the estimated error is less than the epsabs OR is less
0121       than |I| * epsrel.
0122       Unless there is severe cancellation of positive and negative values of
0123       f(x) over the interval [A,B], the relative error may be considered as
0124       specifying a bound on the <I>relative</I> error of I in the case
0125       |I|&gt;1, and a bound on the absolute error in the case |I|&lt;1. More
0126       precisely, if k is the number of sub-intervals contributing to the
0127       approximation (see Method), and if
0128    \f[
0129       I_{abs} = \int^{B}_{A} \left|f(x)\right|dx
0130    \f]
0131       then the relation
0132    \f[
0133       \frac{\left|G-I\right|}{I_{abs}+k} < EPS
0134    \f]
0135       will nearly always be true, provided the routine terminates without
0136       printing an error message. For functions f having no singularities in
0137       the closed interval [A,B] the accuracy will usually be much higher than
0138       this.
0139
0140     Error handling:
0141       The requested accuracy cannot be obtained (see Method).
0142       The function value is set equal to zero.
0143
0144     Note 1:
0145       Values of the function f(x) at the interval end-points A and B are not
0146       required. The subprogram may therefore be used when these values are
0147       undefined
0148    */
0149    double Integral (double a, double b) override;
0150
0151    /** Returns Integral of function on an infinite interval.
0152       This function computes, to an attempted specified accuracy, the value of the integral:
0153    \f[
0154       I = \int^{\infty}_{-\infty} f(x)dx
0155    \f]
0156       Usage:
0157         In any arithmetic expression, this function has the approximate value
0158         of the integral I.
0159
0160       The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
0161    */
0162    double Integral () override;
0163
0164    /** Returns Integral of function on an upper semi-infinite interval.
0165       This function computes, to an attempted specified accuracy, the value of the integral:
0166    \f[
0167       I = \int^{\infty}_{A} f(x)dx
0168    \f]
0169       Usage:
0170         In any arithmetic expression, this function has the approximate value
0171         of the integral I.
0172         - A: lower end-point of integration interval.
0173
0174       The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
0175    */
0176    double IntegralUp (double a) override;
0177
0178    /** Returns Integral of function on a lower semi-infinite interval.
0179        This function computes, to an attempted specified accuracy, the value of the integral:
0180    \f[
0181       I = \int^{B}_{-\infty} f(x)dx
0182    \f]
0183       Usage:
0184          In any arithmetic expression, this function has the approximate value
0185          of the integral I.
0186          - B: upper end-point of integration interval.
0187
0188       The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
0189    */
0190    double IntegralLow (double b) override;
0191
0192
0193    /** Set integration function (flag control if function must be copied inside).
0194        \@param f Function to be used in the calculations.
0195    */
0196    void SetFunction (const IGenFunction &) override;
0197
0198    /** This method is not implemented. */
0199    double Integral (const std::vector< double > &pts) override;
0200
0201    /** This method is not implemented. */
0202    double IntegralCauchy (double a, double b, double c) override;
0203
0204    ///  get the option used for the integration
0205    ROOT::Math::IntegratorOneDimOptions Options() const override;
0206
0207    // set the options
0208    void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt) override;
0209
0210 private:
0211
0212    /**
0213       Integration surrogate method. Return integral of passed function in  interval [a,b]
0214       Derived class (like GaussLegendreIntegrator)  can re-implement this method to modify to use
0215       an improved algorithm
0216    */
0217    virtual double DoIntegral (double a, double b, const IGenFunction* func);
0218
0219 protected:
0220
0221    static bool fgAbsValue;          ///< AbsValue used for the calculation of the integral
0222    double fEpsRel;                  ///< Relative error.
0223    double fEpsAbs;                  ///< Absolute error.
0224    bool fUsedOnce;                  ///< Bool value to check if the function was at least called once.
0225    double fLastResult;              ///< Result from the last estimation.
0226    double fLastError;               ///< Error from the last estimation.
0227    const IGenFunction* fFunction;   ///< Pointer to function used.
0228
0229 };
0230
0231 /**
0232    Auxiliary inner class for mapping infinite and semi-infinite integrals
0233 */
0234 class IntegrandTransform : public IGenFunction {
0235 public:
0236    enum ESemiInfinitySign {kMinus = -1, kPlus = +1};
0237    IntegrandTransform(const IGenFunction* integrand);
0238    IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand);
0239    double operator()(double x) const;
0240    double DoEval(double x) const override;
0241    IGenFunction* Clone() const override;
0242 private:
0243    ESemiInfinitySign fSign;
0244    const IGenFunction* fIntegrand;
0245    double fBoundary;
0246    bool fInfiniteInterval;
0247    double DoEval(double x, double boundary, int sign) const;
0248 };
0249
0250
0251
0252 } // end namespace Math
0253
0254 } // end namespace ROOT
0255
0256 #endif /* ROOT_Math_GaussIntegrator */