Logo ROOT  
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
36
38
39#include "Math/IFunctionfwd.h"
40
41
42
43
45
46#include <vector>
47
48
49
50namespace ROOT {
51namespace 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
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 assignment 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;
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 */
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
int type
Definition: TGX11.cxx:121
Wrapper class to the gsl_function C structure.
Class for performing numerical integration of a function in one dimension.
Definition: GSLIntegrator.h:90
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options
Integration::GKRule fRule
void SetIntegrationRule(Integration::GKRule)
set the integration rule (Gauss-Kronrod rule).
GSLIntegrator & operator=(const GSLIntegrator &)
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...
IntegrationOneDim::Type GetType() const
get type name
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,...
GSLIntegrationWorkspace * fWorkspace
GSLFunctionWrapper * fFunction
int Status() const
return the Error Status of the last Integral calculation
int NEval() const
return number of function evaluations in calculating the integral
double Error() const
return the estimate of the absolute Error of the last Integral calculation
void SetRelTolerance(double relTolerance)
set the desired relative Error
GSLIntegrator(double absTol=1.E-9, double relTol=1E-6, size_t size=1000)
Default constructor of GSL Integrator for Adaptive Singular integration.
const char * GetTypeName() const
return the name
Integration::Type fType
double Integral()
evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with G...
double Result() const
return the Result of the last Integral calculation
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Numerical one dimensional integration options.
Interface (abstract) class for 1D numerical integration It must be implemented by the concrete Integr...
GKRule
enumeration specifying the Gauss-KronRod integration rule for ADAPTIVE integration type
Type
enumeration specifying the integration types.
Namespace for new Math classes and functions.
double(* GSLFuncPointer)(double, void *)
Function pointer corresponding to gsl_function signature.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:96
auto * a
Definition: textangle.C:12