ROOT  6.06/09
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  Begin_Latex
96  I = #int^{b}_{a} f(x)dx
97  End_Latex
98  and define
99  Begin_Latex
100  r(a,b) = #frac{#||{g_{16}(a,b)-g_{8}(a,b)}}{1+#||{g_{16}(a,b)}}
101  End_Latex
102  Then,
103  Begin_Latex
104  G = #sum_{i=1}^{k}g_{16}(x_{i-1},x_{i})
105  End_Latex
106  where, starting with x0 = A and finishing with xk = B,
107  the subdivision points xi(i=1,2,...) are given by
108  Begin_Latex
109  x_{i} = x_{i-1} + #lambda(B-x_{i-1})
110  End_Latex
111  Begin_Latex
112  #lambda
113  End_Latex
114  is equal to the first member of the
115  sequence 1,1/2,1/4,... for which r(xi-1, xi) < EPS.
116  If, at any stage in the process of subdivision, the ratio
117  Begin_Latex
118  q = #||{#frac{x_{i}-x_{i-1}}{B-A}}
119  End_Latex
120  is so small that 1+0.005q is indistinguishable from 1 to
121  machine accuracy, an error exit occurs with the function value
122  set equal to zero.
123 
124  Accuracy:
125  The user provides absolute and relative error bounds (epsrel and epsabs) and the
126  algorithm will stop when the estimated error is less than the epsabs OR is less
127  than |I| * epsrel.
128  Unless there is severe cancellation of positive and negative values of
129  f(x) over the interval [A,B], the relative error may be considered as
130  specifying a bound on the <I>relative</I> error of I in the case
131  |I|&gt;1, and a bound on the absolute error in the case |I|&lt;1. More
132  precisely, if k is the number of sub-intervals contributing to the
133  approximation (see Method), and if
134  Begin_Latex
135  I_{abs} = #int^{B}_{A} #||{f(x)}dx
136  End_Latex
137  then the relation
138  Begin_Latex
139  #frac{#||{G-I}}{I_{abs}+k} < EPS
140  End_Latex
141  will nearly always be true, provided the routine terminates without
142  printing an error message. For functions f having no singularities in
143  the closed interval [A,B] the accuracy will usually be much higher than
144  this.
145 
146  Error handling:
147  The requested accuracy cannot be obtained (see Method).
148  The function value is set equal to zero.
149 
150  Note 1:
151  Values of the function f(x) at the interval end-points A and B are not
152  required. The subprogram may therefore be used when these values are
153  undefined
154  */
155  double Integral (double a, double b);
156 
157  /** Returns Integral of function on an infinite interval.
158  This function computes, to an attempted specified accuracy, the value of the integral:
159  Begin_Latex
160  I = #int^{#infinity}_{-#infinity} f(x)dx
161  End_Latex
162  Usage:
163  In any arithmetic expression, this function has the approximate value
164  of the integral I.
165 
166  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
167  */
168  double Integral ();
169 
170  /** Returns Integral of function on an upper semi-infinite interval.
171  This function computes, to an attempted specified accuracy, the value of the integral:
172  Begin_Latex
173  I = #int^{#infinity}_{A} f(x)dx
174  End_Latex
175  Usage:
176  In any arithmetic expression, this function has the approximate value
177  of the integral I.
178  - A: lower end-point of integration interval.
179 
180  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
181  */
182  double IntegralUp (double a);
183 
184  /** Returns Integral of function on a lower semi-infinite interval.
185  This function computes, to an attempted specified accuracy, the value of the integral:
186  Begin_Latex
187  I = #int^{B}_{#infinity} f(x)dx
188  End_Latex
189  Usage:
190  In any arithmetic expression, this function has the approximate value
191  of the integral I.
192  - B: upper end-point of integration interval.
193 
194  The integral is mapped onto [0,1] using a transformation then integral computation is surrogated to DoIntegral.
195  */
196  double IntegralLow (double b);
197 
198 
199  /** Set integration function (flag control if function must be copied inside).
200  \@param f Function to be used in the calculations.
201  */
202  void SetFunction (const IGenFunction &);
203 
204  /** This method is not implemented. */
205  double Integral (const std::vector< double > &pts);
206 
207  /** This method is not implemented. */
208  double IntegralCauchy (double a, double b, double c);
209 
210  /// get the option used for the integration
212 
213  // set the options
214  virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt);
215 
216 private:
217 
218  /**
219  Integration surrugate method. Return integral of passed function in interval [a,b]
220  Derived class (like GaussLegendreIntegrator) can re-implement this method to modify to use
221  an improved algorithm
222  */
223  virtual double DoIntegral (double a, double b, const IGenFunction* func);
224 
225 protected:
226 
227  static bool fgAbsValue; // AbsValue used for the calculation of the integral
228  double fEpsRel; // Relative error.
229  double fEpsAbs; // Absolute error.
230  bool fUsedOnce; // Bool value to check if the function was at least called once.
231  double fLastResult; // Result from the last stimation.
232  double fLastError; // Error from the last stimation.
233  const IGenFunction* fFunction; // Pointer to function used.
234 
235 };
236 
237 /**
238  Auxillary inner class for mapping infinite and semi-infinite integrals
239 */
241 public:
242  enum ESemiInfinitySign {kMinus = -1, kPlus = +1};
243  IntegrandTransform(const IGenFunction* integrand);
244  IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand);
245  double operator()(double x) const;
246  double DoEval(double x) const;
247  IGenFunction* Clone() const;
248 private:
251  double fBoundary;
253  double DoEval(double x, double boundary, int sign) const;
254 };
255 
256 
257 
258 } // end namespace Math
259 
260 } // end namespace ROOT
261 
262 #endif /* ROOT_Math_GaussIntegrator */
double Error() const
Return the estimate of the absolute Error of the last Integral calculation.
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.
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
const double absTol
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
double Integral()
Returns Integral of function on an infinite interval.
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 surrugate 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.
double operator()(double x) const
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 DoEval(double x) const
implementation of the evaluation function. Must be implemented by derived classes ...
void SetFunction(const IGenFunction &)
Set integration function (flag control if function must be copied inside).
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Auxillary inner class for mapping infinite and semi-infinite integrals.
Namespace for new Math classes and functions.
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
IGenFunction * Clone() const
Clone a function.
int Status() const
return the status of the last integration - 0 in case of success
Int_t sign(Double_t x)
Definition: CsgOps.cxx:89
double Result() const
Returns the result of the last Integral calculation.
IntegrandTransform(const IGenFunction *integrand)
const IGenFunction * fFunction