Logo ROOT   6.08/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 
18 namespace ROOT {
19 namespace Math {
20 
21 
22 bool GaussIntegrator::fgAbsValue = false;
23 
24  GaussIntegrator::GaussIntegrator(double epsabs, double epsrel)
25 {
26 // Default Constructor. If tolerances are not given use default values from ROOT::Math::IntegratorOneDimOptions
27 
28  fEpsAbs = epsabs;
29  fEpsRel = epsrel;
31  if (epsrel < 0 || (epsabs == 0 && epsrel == 0)) fEpsRel = ROOT::Math::IntegratorOneDimOptions::DefaultRelTolerance();
32  if (std::max(fEpsRel,fEpsAbs) <= 0.0 ) {
33  fEpsRel = 1.E-9;
34  fEpsAbs = 1.E-9;
35  MATH_WARN_MSG("ROOT::Math::GausIntegrator", "Invalid tolerance given, use values of 1.E-9");
36  }
37 
38  fLastResult = fLastError = 0;
39  fUsedOnce = false;
40  fFunction = 0;
41 }
42 
44 {
45  // Destructor. (no - operations)
46 }
47 
49 { fgAbsValue = flag; }
50 
51 double GaussIntegrator::Integral(double a, double b) {
52  return DoIntegral(a, b, fFunction);
53 }
54 
56  IntegrandTransform it(this->fFunction);
57  return DoIntegral(0., 1., it.Clone());
58 }
59 
60 double GaussIntegrator::IntegralUp (double a) {
62  return DoIntegral(0., 1., it.Clone());
63 }
64 
67  return DoIntegral(0., 1., it.Clone());
68 }
69 
70 double GaussIntegrator::DoIntegral(double a, double b, const IGenFunction* function)
71 {
72  // Return Integral of function between a and b.
73 
74  if (fEpsRel <=0 || fEpsAbs <= 0) {
75  if (fEpsRel > 0) fEpsAbs = fEpsRel;
76  else if (fEpsAbs > 0) fEpsRel = fEpsAbs;
77  else {
78  MATH_INFO_MSG("ROOT::Math::GausIntegratorOneDim", "Invalid tolerance given - use default values");
81  }
82  }
83 
84  const double kHF = 0.5;
85  const double kCST = 5./1000;
86 
87  double x[12] = { 0.96028985649753623, 0.79666647741362674,
88  0.52553240991632899, 0.18343464249564980,
89  0.98940093499164993, 0.94457502307323258,
90  0.86563120238783174, 0.75540440835500303,
91  0.61787624440264375, 0.45801677765722739,
92  0.28160355077925891, 0.09501250983763744};
93 
94  double w[12] = { 0.10122853629037626, 0.22238103445337447,
95  0.31370664587788729, 0.36268378337836198,
96  0.02715245941175409, 0.06225352393864789,
97  0.09515851168249278, 0.12462897125553387,
98  0.14959598881657673, 0.16915651939500254,
99  0.18260341504492359, 0.18945061045506850};
100 
101  double h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2;
102  double xx[1];
103  int i;
104 
105  if ( fFunction == 0 )
106  {
107  MATH_ERROR_MSG("ROOT::Math::GausIntegratorOneDim", "A function must be set first!");
108  return 0.0;
109  }
110 
111  h = 0;
112  fUsedOnce = true;
113  if (b == a) return h;
114  aconst = kCST/std::abs(b-a);
115  bb = a;
116 CASE1:
117  aa = bb;
118  bb = b;
119 CASE2:
120  c1 = kHF*(bb+aa);
121  c2 = kHF*(bb-aa);
122  s8 = 0;
123  for (i=0;i<4;i++) {
124  u = c2*x[i];
125  xx[0] = c1+u;
126  f1 = (*function)(xx);
127  if (fgAbsValue) f1 = std::abs(f1);
128  xx[0] = c1-u;
129  f2 = (*function) (xx);
130  if (fgAbsValue) f2 = std::abs(f2);
131  s8 += w[i]*(f1 + f2);
132  }
133  s16 = 0;
134  for (i=4;i<12;i++) {
135  u = c2*x[i];
136  xx[0] = c1+u;
137  f1 = (*function) (xx);
138  if (fgAbsValue) f1 = std::abs(f1);
139  xx[0] = c1-u;
140  f2 = (*function) (xx);
141  if (fgAbsValue) f2 = std::abs(f2);
142  s16 += w[i]*(f1 + f2);
143  }
144  s16 = c2*s16;
145  //if (std::abs(s16-c2*s8) <= fEpsilon*(1. + std::abs(s16))) {
146  double error = std::abs(s16-c2*s8);
147  if (error <= fEpsAbs || error <= fEpsRel*std::abs(s16)) {
148  h += s16;
149  if(bb != b) goto CASE1;
150  } else {
151  bb = c1;
152  if(1. + aconst*std::abs(c2) != 1) goto CASE2;
153  double maxtol = std::max(fEpsRel, fEpsAbs);
154  MATH_WARN_MSGVAL("ROOT::Math::GausIntegrator", "Failed to reach the desired tolerance ",maxtol);
155  h = s8; //this is a crude approximation (cernlib function returned 0 !)
156  }
157 
158  fLastResult = h;
159  fLastError = std::abs(s16-c2*s8);
160 
161  return h;
162 }
163 
164 
165 double GaussIntegrator::Result () const
166 {
167  // Returns the result of the last Integral calculation.
168 
169  if (!fUsedOnce)
170  MATH_ERROR_MSG("ROOT::Math::GaussIntegrator", "You must calculate the result at least once!");
171 
172  return fLastResult;
173 }
174 
176 { return fLastError; }
177 
179 { return (fUsedOnce) ? 0 : -1; }
180 
182 {
183  // Set integration function
184  fFunction = &function;
185  // reset fUsedOne flag
186  fUsedOnce = false;
187 }
188 
189 double GaussIntegrator::Integral (const std::vector< double > &/*pts*/)
190 {
191  // not impl.
192  MATH_WARN_MSG("ROOT::Math::GaussIntegrator", "This method is not implemented in this class !");
193  return -1.0;
194 }
195 
196 double GaussIntegrator::IntegralCauchy (double /*a*/, double /*b*/, double /*c*/)
197 {
198  // not impl.
199  MATH_WARN_MSG("ROOT::Math::GaussIntegrator", "This method is not implemented in this class !");
200  return -1.0;
201 }
202 
204  // set options
207 }
208 
210  // retrieve options
212  opt.SetIntegrator("Gauss");
215  opt.SetWKSize(0);
216  opt.SetNPoints(0);
217  return opt;
218 }
219 
220 
221 
222 //implementation of IntegrandTransform class
223 
225  : fSign(kPlus), fIntegrand(integrand), fBoundary(0.), fInfiniteInterval(true) {
226 }
227 
228 IntegrandTransform::IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand)
229  : fSign(sign), fIntegrand(integrand), fBoundary(boundary), fInfiniteInterval(false) {
230 }
231 
232 double IntegrandTransform::DoEval(double x) const {
233  double result = DoEval(x, fBoundary, fSign);
234  return (result += (fInfiniteInterval ? DoEval(x, 0., -1) : 0.));
235 }
236 
237 double IntegrandTransform::DoEval(double x, double boundary, int sign) const {
238  double mappedX = 1. / x - 1.;
239  return (*fIntegrand)(boundary + sign * mappedX) * std::pow(mappedX + 1., 2);;
240 }
241 
242 double IntegrandTransform::operator()(double x) const {
243  return DoEval(x);
244 }
245 
248 }
249 
250 
251 } // end namespace Math
252 } // 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:133
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 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.
TArc * a
Definition: textangle.C:12
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.
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.
double f2(const double *x)
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
double result[121]
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