 ROOT   6.08/07 Reference Guide
GaussIntegrator.h
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: David Gonzalez Maline 01/2008
3
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10
11 // Header file for GaussIntegrator
12 //
13 // Created by: David Gonzalez Maline : Wed Jan 16 2008
14 //
15
16 #ifndef ROOT_Math_GaussIntegrator
17 #define ROOT_Math_GaussIntegrator
18
19 #ifndef ROOT_Math_IFunction
20 #include "Math/IFunction.h"
21 #endif
22
23 #ifndef ROOT_Math_VirtualIntegrator
24 #include "Math/VirtualIntegrator.h"
25 #endif
26
27
28 namespace ROOT {
29 namespace Math {
30
31
32 //___________________________________________________________________________________________
33 /**
34  User class for performing function integration.
35
36  It will use the Gauss Method for function integration in a given interval.
37  This class is implemented from TF1::Integral().
38
39  @ingroup Integration
40
41  */
42
44
45
46 public:
47
48  /** Destructor */
49  virtual ~GaussIntegrator();
50
51  /** Default Constructor.
52  If the tolerance are not given, use default values specified in ROOT::Math::IntegratorOneDimOptions
53  */
54  GaussIntegrator(double absTol = -1, double relTol = -1);
55
56
57  /** Static function: set the fgAbsValue flag.
58  By default TF1::Integral uses the original function value to compute the integral
59  However, TF1::Moment, CentralMoment require to compute the integral
60  using the absolute value of the function.
61  */
62  void AbsValue(bool flag);
63
64
65  // Implementing VirtualIntegrator Interface
66
67  /** Set the desired relative Error. */
68  virtual void SetRelTolerance (double eps) { fEpsRel = eps; }
69
70  /** This method is not implemented. */
71  virtual void SetAbsTolerance (double eps) { fEpsAbs = eps; }
72
73  /** Returns the result of the last Integral calculation. */
74  double Result () const;
75
76  /** Return the estimate of the absolute Error of the last Integral calculation. */
77  double Error () const;
78
79  /** return the status of the last integration - 0 in case of success */
80  int Status () const;
81
82  // Implementing VirtualIntegratorOneDim Interface
83
84  /**
85  Returns Integral of function between a and b.
86  Based on original CERNLIB routine DGAUSS by Sigfried Kolbig
87  converted to C++ by Rene Brun
88
89  This function computes, to an attempted specified accuracy, the value
90  of the integral.
91
92  Method:
93  For any interval [a,b] we define g8(a,b) and g16(a,b) to be the 8-point
94  and 16-point Gaussian quadrature approximations to
95  \f[
96  I = \int^{b}_{a} f(x)dx
97  \f]
98  and define
99  \f[
100  r(a,b) = \frac{\left|g_{16}(a,b)-g_{8}(a,b)\right|}{1+\left|g_{16}(a,b)\right|}
101  \f]
102  Then,
103  \f[
104  G = \sum_{i=1}^{k}g_{16}(x_{i-1},x_{i})
105  \f]
106  where, starting with \f$x_{0} = A\f$ and finishing with \f$x_{k} = B\f$,
107  the subdivision points \f$x_{i}(i=1,2,...)\f$ are given by
108  \f[
109  x_{i} = x_{i-1} + \lambda(B-x_{i-1})
110  \f]
111  \f$\lambda\f$ is equal to the first member of the
112  sequence 1,1/2,1/4,... for which \f$r(x_{i-1}, x_{i}) < EPS\f$.
113  If, at any stage in the process of subdivision, the ratio
114  \f[
115  q = \left|\frac{x_{i}-x_{i-1}}{B-A}\right|
116  \f]
117  is so small that 1+0.005q is indistinguishable from 1 to
118  machine accuracy, an error exit occurs with the function value
119  set equal to zero.
120
121  Accuracy:
122  The user provides absolute and relative error bounds (epsrel and epsabs) and the
123  algorithm will stop when the estimated error is less than the epsabs OR is less
124  than |I| * epsrel.
125  Unless there is severe cancellation of positive and negative values of
126  f(x) over the interval [A,B], the relative error may be considered as
127  specifying a bound on the <I>relative</I> error of I in the case
128  |I|&gt;1, and a bound on the absolute error in the case |I|&lt;1. More
129  precisely, if k is the number of sub-intervals contributing to the
130  approximation (see Method), and if
131  \f[
132  I_{abs} = \int^{B}_{A} \left|f(x)\right|dx
133  \f]
134  then the relation
135  \f[
136  \frac{\left|G-I\right|}{I_{abs}+k} < EPS
137  \f]
138  will nearly always be true, provided the routine terminates without
139  printing an error message. For functions f having no singularities in
140  the closed interval [A,B] the accuracy will usually be much higher than
141  this.
142
143  Error handling:
144  The requested accuracy cannot be obtained (see Method).
145  The function value is set equal to zero.
146
147  Note 1:
148  Values of the function f(x) at the interval end-points A and B are not
149  required. The subprogram may therefore be used when these values are
150  undefined
151  */
152  double Integral (double a, double b);
153
154  /** Returns Integral of function on an infinite interval.
155  This function computes, to an attempted specified accuracy, the value of the integral:
156  \f[
157  I = \int^{\infty}_{-\infty} f(x)dx
158  \f]
159  Usage:
160  In any arithmetic expression, this function has the approximate value
161  of the integral I.
162
163  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
164  */
165  double Integral ();
166
167  /** Returns Integral of function on an upper semi-infinite interval.
168  This function computes, to an attempted specified accuracy, the value of the integral:
169  \f[
170  I = \int^{\infty}_{A} f(x)dx
171  \f]
172  Usage:
173  In any arithmetic expression, this function has the approximate value
174  of the integral I.
175  - A: lower end-point of integration interval.
176
177  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
178  */
179  double IntegralUp (double a);
180
181  /** Returns Integral of function on a lower semi-infinite interval.
182  This function computes, to an attempted specified accuracy, the value of the integral:
183  \f[
184  I = \int^{B}_{-\infty} f(x)dx
185  \f]
186  Usage:
187  In any arithmetic expression, this function has the approximate value
188  of the integral I.
189  - B: upper end-point of integration interval.
190
191  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
192  */
193  double IntegralLow (double b);
194
195
196  /** Set integration function (flag control if function must be copied inside).
197  \@param f Function to be used in the calculations.
198  */
199  void SetFunction (const IGenFunction &);
200
201  /** This method is not implemented. */
202  double Integral (const std::vector< double > &pts);
203
204  /** This method is not implemented. */
205  double IntegralCauchy (double a, double b, double c);
206
207  /// get the option used for the integration
209
210  // set the options
211  virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt);
212
213 private:
214
215  /**
216  Integration surrogate method. Return integral of passed function in interval [a,b]
217  Derived class (like GaussLegendreIntegrator) can re-implement this method to modify to use
218  an improved algorithm
219  */
220  virtual double DoIntegral (double a, double b, const IGenFunction* func);
221
222 protected:
223
224  static bool fgAbsValue; // AbsValue used for the calculation of the integral
225  double fEpsRel; // Relative error.
226  double fEpsAbs; // Absolute error.
227  bool fUsedOnce; // Bool value to check if the function was at least called once.
228  double fLastResult; // Result from the last estimation.
229  double fLastError; // Error from the last estimation.
230  const IGenFunction* fFunction; // Pointer to function used.
231
232 };
233
234 /**
235  Auxiliary inner class for mapping infinite and semi-infinite integrals
236 */
238 public:
239  enum ESemiInfinitySign {kMinus = -1, kPlus = +1};
240  IntegrandTransform(const IGenFunction* integrand);
241  IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand);
242  double operator()(double x) const;
243  double DoEval(double x) const;
244  IGenFunction* Clone() const;
245 private:
248  double fBoundary;
250  double DoEval(double x, double boundary, int sign) const;
251 };
252
253
254
255 } // end namespace Math
256
257 } // end namespace ROOT
258
259 #endif /* ROOT_Math_GaussIntegrator */
Interface (abstract) class for 1D numerical integration It must be implemented by the concrate Integr...
virtual void SetAbsTolerance(double eps)
This method is not implemented.
double IntegralUp(double a)
Returns Integral of function on an upper semi-infinite interval.
double Error() const
Return the estimate of the absolute Error of the last Integral calculation.
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
const double absTol
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
double Integral()
Returns Integral of function on an infinite interval.
return c
GaussIntegrator(double absTol=-1, double relTol=-1)
Default Constructor.
TArc * a
Definition: textangle.C:12
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrogate method.
double IntegralLow(double b)
Returns Integral of function on a lower semi-infinite interval.
virtual void SetRelTolerance(double eps)
Set the desired relative Error.
Double_t x[n]
Definition: legend1.C:17
double IntegralCauchy(double a, double b, double c)
This method is not implemented.
int Status() const
return the status of the last integration - 0 in case of success
void AbsValue(bool flag)
Static function: set the fgAbsValue flag.
const IGenFunction * fIntegrand
virtual ~GaussIntegrator()
Destructor.
Numerical one dimensional integration options.
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options (should be re-implemented by derived classes -if more options than tolerance exist ...
User class for performing function integration.
double Result() const
Returns the result of the last Integral calculation.
void SetFunction(const IGenFunction &)
Set integration function (flag control if function must be copied inside).
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Auxiliary inner class for mapping infinite and semi-infinite integrals.
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
const IGenFunction * fFunction