Logo ROOT   6.12/07
Reference Guide
GaussIntegrator.cxx
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 #include "Math/GaussIntegrator.h"
12 #include "Math/Error.h"
13 #include "Math/IntegratorOptions.h"
14 #include "Math/IFunction.h"
15 #include "Math/IFunctionfwd.h"
16 #include <cmath>
17 #include <algorithm>
18 
19 namespace ROOT {
20 namespace Math {
21 
22 
23 bool GaussIntegrator::fgAbsValue = false;
24 
25  GaussIntegrator::GaussIntegrator(double epsabs, double epsrel)
26 {
27 // Default Constructor. If tolerances are not given use default values from ROOT::Math::IntegratorOneDimOptions
28 
29  fEpsAbs = epsabs;
30  fEpsRel = epsrel;
32  if (epsrel < 0 || (epsabs == 0 && epsrel == 0)) fEpsRel = ROOT::Math::IntegratorOneDimOptions::DefaultRelTolerance();
33  if (std::max(fEpsRel,fEpsAbs) <= 0.0 ) {
34  fEpsRel = 1.E-9;
35  fEpsAbs = 1.E-9;
36  MATH_WARN_MSG("ROOT::Math::GausIntegrator", "Invalid tolerance given, use values of 1.E-9");
37  }
38 
39  fLastResult = fLastError = 0;
40  fUsedOnce = false;
41  fFunction = 0;
42 }
43 
45 {
46  // Destructor. (no - operations)
47 }
48 
50 { fgAbsValue = flag; }
51 
52 double GaussIntegrator::Integral(double a, double b) {
53  return DoIntegral(a, b, fFunction);
54 }
55 
57  IntegrandTransform it(this->fFunction);
58  return DoIntegral(0., 1., it.Clone());
59 }
60 
61 double GaussIntegrator::IntegralUp (double a) {
63  return DoIntegral(0., 1., it.Clone());
64 }
65 
68  return DoIntegral(0., 1., it.Clone());
69 }
70 
71 double GaussIntegrator::DoIntegral(double a, double b, const IGenFunction* function)
72 {
73  // Return Integral of function between a and b.
74 
75  if (fEpsRel <=0 || fEpsAbs <= 0) {
76  if (fEpsRel > 0) fEpsAbs = fEpsRel;
77  else if (fEpsAbs > 0) fEpsRel = fEpsAbs;
78  else {
79  MATH_INFO_MSG("ROOT::Math::GausIntegratorOneDim", "Invalid tolerance given - use default values");
82  }
83  }
84 
85  const double kHF = 0.5;
86  const double kCST = 5./1000;
87 
88  double x[12] = { 0.96028985649753623, 0.79666647741362674,
89  0.52553240991632899, 0.18343464249564980,
90  0.98940093499164993, 0.94457502307323258,
91  0.86563120238783174, 0.75540440835500303,
92  0.61787624440264375, 0.45801677765722739,
93  0.28160355077925891, 0.09501250983763744};
94 
95  double w[12] = { 0.10122853629037626, 0.22238103445337447,
96  0.31370664587788729, 0.36268378337836198,
97  0.02715245941175409, 0.06225352393864789,
98  0.09515851168249278, 0.12462897125553387,
99  0.14959598881657673, 0.16915651939500254,
100  0.18260341504492359, 0.18945061045506850};
101 
102  double h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2;
103  double xx[1];
104  int i;
105 
106  if ( fFunction == 0 )
107  {
108  MATH_ERROR_MSG("ROOT::Math::GausIntegratorOneDim", "A function must be set first!");
109  return 0.0;
110  }
111 
112  h = 0;
113  fUsedOnce = true;
114  if (b == a) return h;
115  aconst = kCST/std::abs(b-a);
116  bb = a;
117 CASE1:
118  aa = bb;
119  bb = b;
120 CASE2:
121  c1 = kHF*(bb+aa);
122  c2 = kHF*(bb-aa);
123  s8 = 0;
124  for (i=0;i<4;i++) {
125  u = c2*x[i];
126  xx[0] = c1+u;
127  f1 = (*function)(xx);
128  if (fgAbsValue) f1 = std::abs(f1);
129  xx[0] = c1-u;
130  f2 = (*function) (xx);
131  if (fgAbsValue) f2 = std::abs(f2);
132  s8 += w[i]*(f1 + f2);
133  }
134  s16 = 0;
135  for (i=4;i<12;i++) {
136  u = c2*x[i];
137  xx[0] = c1+u;
138  f1 = (*function) (xx);
139  if (fgAbsValue) f1 = std::abs(f1);
140  xx[0] = c1-u;
141  f2 = (*function) (xx);
142  if (fgAbsValue) f2 = std::abs(f2);
143  s16 += w[i]*(f1 + f2);
144  }
145  s16 = c2*s16;
146  //if (std::abs(s16-c2*s8) <= fEpsilon*(1. + std::abs(s16))) {
147  double error = std::abs(s16-c2*s8);
148  if (error <= fEpsAbs || error <= fEpsRel*std::abs(s16)) {
149  h += s16;
150  if(bb != b) goto CASE1;
151  } else {
152  bb = c1;
153  if(1. + aconst*std::abs(c2) != 1) goto CASE2;
154  double maxtol = std::max(fEpsRel, fEpsAbs);
155  MATH_WARN_MSGVAL("ROOT::Math::GausIntegrator", "Failed to reach the desired tolerance ",maxtol);
156  h = s8; //this is a crude approximation (cernlib function returned 0 !)
157  }
158 
159  fLastResult = h;
160  fLastError = std::abs(s16-c2*s8);
161 
162  return h;
163 }
164 
165 
166 double GaussIntegrator::Result () const
167 {
168  // Returns the result of the last Integral calculation.
169 
170  if (!fUsedOnce)
171  MATH_ERROR_MSG("ROOT::Math::GaussIntegrator", "You must calculate the result at least once!");
172 
173  return fLastResult;
174 }
175 
177 { return fLastError; }
178 
180 { return (fUsedOnce) ? 0 : -1; }
181 
183 {
184  // Set integration function
185  fFunction = &function;
186  // reset fUsedOne flag
187  fUsedOnce = false;
188 }
189 
190 double GaussIntegrator::Integral (const std::vector< double > &/*pts*/)
191 {
192  // not impl.
193  MATH_WARN_MSG("ROOT::Math::GaussIntegrator", "This method is not implemented in this class !");
194  return -1.0;
195 }
196 
197 double GaussIntegrator::IntegralCauchy (double /*a*/, double /*b*/, double /*c*/)
198 {
199  // not impl.
200  MATH_WARN_MSG("ROOT::Math::GaussIntegrator", "This method is not implemented in this class !");
201  return -1.0;
202 }
203 
205  // set options
208 }
209 
211  // retrieve options
213  opt.SetIntegrator("Gauss");
216  opt.SetWKSize(0);
217  opt.SetNPoints(0);
218  return opt;
219 }
220 
221 
222 
223 //implementation of IntegrandTransform class
224 
226  : fSign(kPlus), fIntegrand(integrand), fBoundary(0.), fInfiniteInterval(true) {
227 }
228 
229 IntegrandTransform::IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand)
230  : fSign(sign), fIntegrand(integrand), fBoundary(boundary), fInfiniteInterval(false) {
231 }
232 
233 double IntegrandTransform::DoEval(double x) const {
234  double result = DoEval(x, fBoundary, fSign);
235  return (result += (fInfiniteInterval ? DoEval(x, 0., -1) : 0.));
236 }
237 
238 double IntegrandTransform::DoEval(double x, double boundary, int sign) const {
239  double mappedX = 1. / x - 1.;
240  return (*fIntegrand)(boundary + sign * mappedX) * std::pow(mappedX + 1., 2);;
241 }
242 
243 double IntegrandTransform::operator()(double x) const {
244  return DoEval(x);
245 }
246 
249 }
250 
251 
252 } // end namespace Math
253 } // end namespace ROOT
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:135
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double Integral()
Returns Integral of function on an infinite interval.
return c1
Definition: legend1.C:41
void SetWKSize(unsigned int size)
set workspace size
IGenFunction * Clone() const
Clone a function.
TH1 * h
Definition: legend2.C:5
GaussIntegrator(double absTol=-1, double relTol=-1)
Default Constructor.
void SetNPoints(unsigned int n)
set number of points rule values of 1,2,3,4,5,6 corresponds to 15,21,31,41,51,61 and they are used in...
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrogate method.
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:47
double RelTolerance() const
absolute tolerance
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 pow(double, double)
double operator()(double x) const
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
int Status() const
return the status of the last integration - 0 in case of success
double DoEval(double x) const
implementation of the evaluation function. Must be implemented by derived classes ...
void AbsValue(bool flag)
Static function: set the fgAbsValue flag.
#define MATH_INFO_MSG(loc, str)
Definition: Error.h:44
const IGenFunction * fIntegrand
double AbsTolerance() const
non-static methods for retrivieng options
virtual ~GaussIntegrator()
Destructor.
auto * a
Definition: textangle.C:12
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 ...
double Result() const
Returns the result of the last Integral calculation.
return c2
Definition: legend2.C:14
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
#define MATH_WARN_MSGVAL(loc, str, x)
Definition: Error.h:67
Auxiliary inner class for mapping infinite and semi-infinite integrals.
Namespace for new Math classes and functions.
TF1 * f1
Definition: legend1.C:11
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
IntegrandTransform(const IGenFunction *integrand)
void SetIntegrator(const char *name)
set 1D integrator name
void SetRelTolerance(double tol)
set the relative tolerance
const IGenFunction * fFunction
void SetAbsTolerance(double tol)
non-static methods for setting options