Logo ROOT   6.10/09
Reference Guide
GSLIntegrator.h
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Header file for class GSLIntegrator
26 //
27 // Created by: Lorenzo Moneta at Thu Nov 11 14:22:32 2004
28 //
29 // Last update: Thu Nov 11 14:22:32 2004
30 //
31 #ifndef ROOT_Math_GSLIntegrator
32 #define ROOT_Math_GSLIntegrator
33 
34 
35 #include "Math/VirtualIntegrator.h"
36 
37 #include "Math/IntegrationTypes.h"
38 
39 #include "Math/IFunctionfwd.h"
40 
41 
42 
43 
45 
46 #include <vector>
47 
48 
49 
50 namespace ROOT {
51 namespace Math {
52 
53 
54 
55  class GSLIntegrationWorkspace;
56  class GSLFunctionWrapper;
57 
58  //_________________________________________________________________________
59  /**
60 
61  Class for performing numerical integration of a function in one dimension.
62  It uses the numerical integration algorithms of GSL, which reimplements the
63  algorithms used in the QUADPACK, a numerical integration package written in Fortran.
64 
65  Various types of adaptive and non-adaptive integration are supported. These include
66  integration over infinite and semi-infinite ranges and singular integrals.
67 
68  The integration type is selected using the Integration::type enumeration
69  in the class constructor.
70  The default type is adaptive integration with singularity
71  (ADAPTIVESINGULAR or QAGS in the QUADPACK convention) applying a Gauss-Kronrod 21-point integration rule.
72  In the case of ADAPTIVE type, the integration rule can also be specified via the
73  Integration::GKRule. The default rule is 31 points.
74 
75  In the case of integration over infinite and semi-infinite ranges, the type used is always
76  ADAPTIVESINGULAR applying a transformation from the original interval into (0,1).
77 
78  The ADAPTIVESINGULAR type is the most sophicticated type. When performances are
79  important, it is then recommened to use the NONADAPTIVE type in case of smooth functions or
80  ADAPTIVE with a lower Gauss-Kronrod rule.
81 
82  For detailed description on GSL integration algorithms see the
83  <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html">GSL Manual</A>.
84 
85 
86  @ingroup Integration
87  */
88 
89 
91 
92  public:
93 
94 
95 
96  // constructors
97 
98 
99  /** Default constructor of GSL Integrator for Adaptive Singular integration
100 
101  @param absTol desired absolute Error
102  @param relTol desired relative Error
103  @param size maximum number of sub-intervals
104  */
105 
106  GSLIntegrator(double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
107 
108 
109 
110 
111  /** constructor of GSL Integrator. In the case of Adaptive integration the Gauss-Krond rule of 31 points is used
112 
113  @param type type of integration. The possible types are defined in the Integration::Type enumeration
114  @param absTol desired absolute Error
115  @param relTol desired relative Error
116  @param size maximum number of sub-intervals
117  */
118 
119 
120  GSLIntegrator(const Integration::Type type, double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
121 
122 
123  /**
124  generic constructor for GSL Integrator
125 
126  @param type type of integration. The possible types are defined in the Integration::Type enumeration
127  @param rule Gauss-Kronrod rule. It is used only for ADAPTIVE::Integration types. The possible rules are defined in the Integration::GKRule enumeration
128  @param absTol desired absolute Error
129  @param relTol desired relative Error
130  @param size maximum number of sub-intervals
131 
132  */
133 
134  GSLIntegrator(const Integration::Type type, const Integration::GKRule rule, double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
135 
136 
137  /** constructor of GSL Integrator. In the case of Adaptive integration the Gauss-Krond rule of 31 points is used
138  This is used by the plug-in manager (need a char * instead of enumerations)
139 
140  @param type type of integration. The possible types are defined in the Integration::Type enumeration
141  @param rule Gauss-Kronrod rule (from 1 to 6)
142  @param absTol desired absolute Error
143  @param relTol desired relative Error
144  @param size maximum number of sub-intervals
145  */
146  GSLIntegrator(const char * type, int rule, double absTol, double relTol, size_t size );
147 
148  virtual ~GSLIntegrator();
149  //~GSLIntegrator();
150 
151  // disable copy ctrs
152  private:
153 
154  GSLIntegrator(const GSLIntegrator &);
156 
157  public:
158 
159 
160  // template methods for generic functors
161 
162  /**
163  method to set the a generic integration function
164 
165  @param f integration function. The function type must implement the assigment operator, <em> double operator() ( double x ) </em>
166 
167  */
168 
169 
170  void SetFunction(const IGenFunction &f);
171 
172  /**
173  Set function from a GSL pointer function type
174  */
175  void SetFunction( GSLFuncPointer f, void * p = 0);
176 
177  // methods using IGenFunction
178 
179  /**
180  evaluate the Integral of a function f over the defined interval (a,b)
181  @param f integration function. The function type must implement the mathlib::IGenFunction interface
182  @param a lower value of the integration interval
183  @param b upper value of the integration interval
184  */
185 
186  double Integral(const IGenFunction & f, double a, double b);
187 
188 
189  /**
190  evaluate the Integral of a function f over the infinite interval (-inf,+inf)
191  @param f integration function. The function type must implement the mathlib::IGenFunction interface
192  */
193 
194  double Integral(const IGenFunction & f);
195 
196  /**
197  evaluate the Cauchy principal value of the integral of a previously defined function f over
198  the defined interval (a,b) with a singularity at c
199  @param a lower interval value
200  @param b lower interval value
201  @param c singular value of f
202  */
203  double IntegralCauchy(double a, double b, double c);
204 
205  /**
206  evaluate the Cauchy principal value of the integral of a function f over the defined interval (a,b)
207  with a singularity at c
208  @param f integration function. The function type must implement the mathlib::IGenFunction interface
209  @param a lower interval value
210  @param b lower interval value
211  @param c singular value of f
212  */
213  double IntegralCauchy(const IGenFunction & f, double a, double b, double c);
214 
215  /**
216  evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
217  @param f integration function. The function type must implement the mathlib::IGenFunction interface
218  @param a lower value of the integration interval
219 
220  */
221  double IntegralUp(const IGenFunction & f, double a );
222 
223  /**
224  evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
225  @param f integration function. The function type must implement the mathlib::IGenFunction interface
226  @param b upper value of the integration interval
227  */
228  double IntegralLow(const IGenFunction & f, double b );
229 
230  /**
231  evaluate the Integral of a function f with known singular points over the defined Integral (a,b)
232  @param f integration function. The function type must implement the mathlib::IGenFunction interface
233  @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
234 
235  */
236  double Integral(const IGenFunction & f, const std::vector<double> & pts );
237 
238  // evaluate using cached function
239 
240  /**
241  evaluate the Integral over the defined interval (a,b) using the function previously set with GSLIntegrator::SetFunction method
242  @param a lower value of the integration interval
243  @param b upper value of the integration interval
244  */
245 
246  double Integral(double a, double b);
247 
248 
249  /**
250  evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with GSLIntegrator::SetFunction method.
251  */
252  double Integral( );
253 
254  /**
255  evaluate the Integral of a function f over the semi-infinite interval (a,+inf) using the function previously set with GSLIntegrator::SetFunction method.
256  @param a lower value of the integration interval
257  */
258  double IntegralUp(double a );
259 
260  /**
261  evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b) using the function previously set with GSLIntegrator::SetFunction method.
262  @param b upper value of the integration interval
263  */
264  double IntegralLow( double b );
265 
266  /**
267  evaluate the Integral over the defined interval (a,b) using the function previously set with GSLIntegrator::SetFunction method. The function has known singular points.
268  @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
269 
270  */
271  double Integral( const std::vector<double> & pts);
272 
273  // evaluate using free function pointer (same GSL signature)
274 
275  /**
276  signature for function pointers used by GSL
277  */
278  //typedef double ( * GSLFuncPointer ) ( double, void * );
279 
280  /**
281  evaluate the Integral of of a function f over the defined interval (a,b) passing a free function pointer
282  The integration function must be a free function and have a signature consistent with GSL functions:
283 
284  <em>double my_function ( double x, void * p ) { ...... } </em>
285 
286  This method is the most efficient since no internal adapter to GSL function is created.
287  @param f pointer to the integration function
288  @param p pointer to the Parameters of the function
289  @param a lower value of the integration interval
290  @param b upper value of the integration interval
291 
292  */
293  double Integral(GSLFuncPointer f, void * p, double a, double b);
294 
295  /**
296  evaluate the Integral of a function f over the infinite interval (-inf,+inf) passing a free function pointer
297  */
298  double Integral(GSLFuncPointer f, void * p);
299 
300  /**
301  evaluate the Integral of a function f over the semi-infinite interval (a,+inf) passing a free function pointer
302  */
303  double IntegralUp(GSLFuncPointer f, void * p, double a);
304 
305  /**
306  evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b) passing a free function pointer
307  */
308  double IntegralLow(GSLFuncPointer f, void * p, double b);
309 
310  /**
311  evaluate the Integral of a function f with knows singular points over the over a defined interval passing a free function pointer
312  */
313  double Integral(GSLFuncPointer f, void * p, const std::vector<double> & pts);
314 
315  /**
316  return the Result of the last Integral calculation
317  */
318  double Result() const;
319 
320  /**
321  return the estimate of the absolute Error of the last Integral calculation
322  */
323  double Error() const;
324 
325  /**
326  return the Error Status of the last Integral calculation
327  */
328  int Status() const;
329 
330  /**
331  return number of function evaluations in calculating the integral
332  */
333  int NEval() const { return fNEval; }
334 
335  // setter for control Parameters (getters are not needed so far )
336 
337  /**
338  set the desired relative Error
339  */
340  void SetRelTolerance(double relTolerance);
341 
342 
343  /**
344  set the desired absolute Error
345  */
346  void SetAbsTolerance(double absTolerance);
347 
348  /**
349  set the integration rule (Gauss-Kronrod rule).
350  The possible rules are defined in the Integration::GKRule enumeration.
351  The integration rule can be modified only for ADAPTIVE type integrations
352  */
354 
355  /// set the options
356  virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt);
357 
358  /// get the option used for the integration
360 
361  /// get type name
363 
364  /**
365  return the name
366  */
367  const char * GetTypeName() const;
368 
369 
370  protected:
371 
372  // internal method to check validity of GSL function pointer
373  bool CheckFunction();
374 
375 
376  private:
377 
380  double fAbsTol;
381  double fRelTol;
382  size_t fSize;
384 
385  // cache Error, Result and Status of integration
386 
387  double fResult;
388  double fError;
389  int fStatus;
390  int fNEval;
391 
392  // GSLIntegrationAlgorithm * fAlgorithm;
393 
396 
397  };
398 
399 
400 
401 
402 
403 } // namespace Math
404 } // namespace ROOT
405 
406 
407 #endif /* ROOT_Math_GSLIntegrator */
double Result() const
return the Result of the last Integral calculation
GKRule
enumeration specifying the Gauss-KronRod integration rule for ADAPTIVE integration type ...
Interface (abstract) class for 1D numerical integration It must be implemented by the concrate Integr...
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:134
const double absTol
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
void SetIntegrationRule(Integration::GKRule)
set the integration rule (Gauss-Kronrod rule).
TArc * a
Definition: textangle.C:12
int Status() const
return the Error Status of the last Integral calculation
GSLFunctionWrapper * fFunction
IntegrationOneDim::Type GetType() const
get type name
int NEval() const
return number of function evaluations in calculating the integral
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options
Class for performing numerical integration of a function in one dimension.
Definition: GSLIntegrator.h:90
Integration::Type fType
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
double(* GSLFuncPointer)(double, void *)
Function pointer corresponding to gsl_function signature.
double Integral()
evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with G...
Numerical one dimensional integration options.
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
double Error() const
return the estimate of the absolute Error of the last Integral calculation
constexpr Double_t E()
Definition: TMath.h:74
Type
enumeration specifying the integration types.
GSLIntegrator(double absTol=1.E-9, double relTol=1E-6, size_t size=1000)
Default constructor of GSL Integrator for Adaptive Singular integration.
void SetRelTolerance(double relTolerance)
set the desired relative Error
double f(double x)
int type
Definition: TGX11.cxx:120
GSLIntegrationWorkspace * fWorkspace
double IntegralCauchy(double a, double b, double c)
evaluate the Cauchy principal value of the integral of a previously defined function f over the defin...
Namespace for new Math classes and functions.
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
const char * GetTypeName() const
return the name
GSLIntegrator & operator=(const GSLIntegrator &)
Wrapper class to the gsl_function C structure.
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
Integration::GKRule fRule