ROOT   Reference Guide
Searching...
No Matches
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
12#include "Math/Error.h"
14#include "Math/IFunction.h"
15#include "Math/IFunctionfwd.h"
16#include <cmath>
17#include <algorithm>
18
19namespace ROOT {
20namespace Math {
21
22
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
40 fUsedOnce = false;
41 fFunction = 0;
42}
43
45{
46 // Destructor. (no - operations)
47}
48
50{ fgAbsValue = flag; }
51
52double GaussIntegrator::Integral(double a, double b) {
53 return DoIntegral(a, b, fFunction);
54}
55
58 return DoIntegral(0., 1., &it);
59}
60
63 return DoIntegral(0., 1., &it);
64}
65
68 return DoIntegral(0., 1., &it);
69}
70
71double 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;
117CASE1:
118 aa = bb;
119 bb = b;
120CASE2:
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
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
190double 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
197double 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
229IntegrandTransform::IntegrandTransform(const double boundary, ESemiInfinitySign sign, const IGenFunction* integrand)
230 : fSign(sign), fIntegrand(integrand), fBoundary(boundary), fInfiniteInterval(false) {
231}
232
233double IntegrandTransform::DoEval(double x) const {
234 double result = DoEval(x, fBoundary, fSign);
235 return (result += (fInfiniteInterval ? DoEval(x, 0., -1) : 0.));
236}
237
238double 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
243double IntegrandTransform::operator()(double x) const {
244 return DoEval(x);
245}
246
249}
250
251
252} // end namespace Math
253} // end namespace ROOT
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition Error.h:77
#define MATH_ERROR_MSG(loc, str)
Definition Error.h:83
#define MATH_WARN_MSG(loc, str)
Definition Error.h:80
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition Error.h:105
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
@ kPlus
Definition TAttMarker.h:49
void SetAbsTolerance(double tol)
non-static methods for setting options
double RelTolerance() const
absolute tolerance
void SetRelTolerance(double tol)
set the relative tolerance
double AbsTolerance() const
non-static methods for retrivieng options
void SetWKSize(unsigned int size)
set workspace size
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options (should be re-implemented by derived classes -if more options than tolerance exist
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
double Result() const
Returns the result of the last Integral calculation.
double Integral()
Returns Integral of function on an infinite interval.
double IntegralLow(double b)
Returns Integral of function on a lower semi-infinite interval.
double IntegralCauchy(double a, double b, double c)
This method is not implemented.
double Error() const
Return the estimate of the absolute Error of the last Integral calculation.
void SetFunction(const IGenFunction &)
Set integration function (flag control if function must be copied inside).
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrogate method.
GaussIntegrator(double absTol=-1, double relTol=-1)
Default Constructor.
virtual void SetAbsTolerance(double eps)
This method is not implemented.
virtual ~GaussIntegrator()
Destructor.
virtual void SetRelTolerance(double eps)
Set the desired relative Error.
const IGenFunction * fFunction
void AbsValue(bool flag)
Static function: set the fgAbsValue flag.
int Status() const
return the status of the last integration - 0 in case of success
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:135
Auxiliary inner class for mapping infinite and semi-infinite integrals.
double DoEval(double x) const
implementation of the evaluation function. Must be implemented by derived classes
IntegrandTransform(const IGenFunction *integrand)
const IGenFunction * fIntegrand
IGenFunction * Clone() const
Clone a function.
double operator()(double x) const
Numerical one dimensional integration options.
void SetIntegrator(const char *name)
set 1D integrator name
void SetNPoints(unsigned int n)
Set number of points for active integration rule.
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
TF1 * f1
Definition legend1.C:11
return c2
Definition legend2.C:14
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...