Logo ROOT  
Reference Guide
GaussLegendreIntegrator.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"
13 #include "Math/IFunction.h"
14 #include "Math/IFunctionfwd.h"
15 #include "Math/IntegratorOptions.h"
16 #include <cmath>
17 #include <string.h>
18 #include <algorithm>
19 
20 namespace ROOT {
21 namespace Math {
22 
24  GaussIntegrator(eps, eps)
25 {
26  // Basic contructor
27  fNum = num;
28  fX = 0;
29  fW = 0;
30 
31  CalcGaussLegendreSamplingPoints();
32 }
33 
35 {
36  // Default Destructor
37 
38 
39  delete [] fX;
40  delete [] fW;
41 }
42 
44 {
45  // Set the number of points used in the calculation of the integral
46 
47  fNum = num;
49 }
50 
51 void GaussLegendreIntegrator::GetWeightVectors(double *x, double *w) const
52 {
53  // Returns the arrays x and w.
54 
55  std::copy(fX,fX+fNum, x);
56  std::copy(fW,fW+fNum, w);
57 }
58 
59 
60 double GaussLegendreIntegrator::DoIntegral(double a, double b, const IGenFunction* function)
61 {
62  // Gauss-Legendre integral, see CalcGaussLegendreSamplingPoints.
63 
64  if (fNum<=0 || fX == 0 || fW == 0)
65  return 0;
66 
67  fUsedOnce = true;
68 
69  const double a0 = (b + a)/2;
70  const double b0 = (b - a)/2;
71 
72  double xx[1];
73 
74  double result = 0.0;
75  for (int i=0; i<fNum; i++)
76  {
77  xx[0] = a0 + b0*fX[i];
78  result += fW[i] * (*function)(xx);
79  }
80 
81  fLastResult = result*b0;
82  return fLastResult;
83 }
84 
85 
87 {
88  // Set the desired relative Error.
89  fEpsRel = eps;
91 }
92 
94 { MATH_WARN_MSG("ROOT::Math::GaussLegendreIntegrator", "There is no Absolute Tolerance!"); }
95 
96 
97 
99 {
100  // Given the number of sampling points this routine fills the
101  // arrays x and w.
102 
103  if (fNum<=0 || fEpsRel<=0)
104  return;
105 
106  if ( fX )
107  delete [] fX;
108 
109  if ( fW )
110  delete [] fW;
111 
112  fX = new double[fNum];
113  fW = new double[fNum];
114 
115  // The roots of symmetric is the interval, so we only have to find half of them
116  const unsigned int m = (fNum+1)/2;
117 
118  double z, pp, p1,p2, p3;
119 
120  // Loop over the desired roots
121  for (unsigned int i=0; i<m; i++) {
122  z = std::cos(3.14159265358979323846*(i+0.75)/(fNum+0.5));
123 
124  // Starting with the above approximation to the i-th root, we enter
125  // the main loop of refinement by Newton's method
126  do {
127  p1=1.0;
128  p2=0.0;
129 
130  // Loop up the recurrence relation to get the Legendre
131  // polynomial evaluated at z
132  for (int j=0; j<fNum; j++)
133  {
134  p3 = p2;
135  p2 = p1;
136  p1 = ((2.0*j+1.0)*z*p2-j*p3)/(j+1.0);
137  }
138  // p1 is now the desired Legendre polynomial. We next compute pp, its
139  // derivative, by a standard relation involving also p2, the polynomial
140  // of one lower order
141  pp = fNum*(z*p1-p2)/(z*z-1.0);
142  // Newton's method
143  z -= p1/pp;
144 
145  } while (std::fabs(p1/pp) > fEpsRel);
146 
147  // Put root and its symmetric counterpart
148  fX[i] = -z;
149  fX[fNum-i-1] = z;
150 
151  // Compute the weight and put its symmetric counterpart
152  fW[i] = 2.0/((1.0-z*z)*pp*pp);
153  fW[fNum-i-1] = fW[i];
154  }
155 }
156 
159  opt.SetAbsTolerance(0);
161  opt.SetWKSize(0);
162  opt.SetNPoints(fNum);
163  opt.SetIntegrator("GaussLegendre");
164  return opt;
165 }
166 
168 {
169  // set integration options
170 // std::cout << "fEpsilon = " << fEpsilon << std::endl;
171 // std::cout << opt.RelTolerance() << " abs " << opt.AbsTolerance() << std::endl;
172  //double tol = opt.RelTolerance(); fEpsilon = tol;
173  fEpsRel = opt.RelTolerance();
174 // std::cout << "fEpsilon = " << fEpsilon << std::endl;
175  fNum = opt.NPoints();
176  if (fNum <= 7) MATH_WARN_MSGVAL("GaussLegendreIntegrator::SetOptions","setting a low number of points ",fNum);
178 }
179 
180 } // end namespace Math
181 } // end namespace ROOT
m
auto * m
Definition: textangle.C:8
IFunction.h
ROOT::Math::IntegratorOneDimOptions::SetNPoints
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...
Definition: IntegratorOptions.h:150
ROOT::Math::GaussLegendreIntegrator::SetNumberPoints
void SetNumberPoints(int num)
Set the number of points used in the calculation of the integral.
Definition: GaussLegendreIntegrator.cxx:53
ROOT::Math::GaussLegendreIntegrator::SetOptions
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options (should be re-implemented by derived classes -if more options than tolerance exist
Definition: GaussLegendreIntegrator.cxx:177
ROOT::Math::GaussLegendreIntegrator::GetWeightVectors
void GetWeightVectors(double *x, double *w) const
Returns the arrays x and w containing the abscissa and weight of the Gauss-Legendre n-point quadratur...
Definition: GaussLegendreIntegrator.cxx:61
ROOT::Math::GaussIntegrator::fLastResult
double fLastResult
Definition: GaussIntegrator.h:240
IFunctionfwd.h
ROOT::Math::IntegratorOneDimOptions::SetIntegrator
void SetIntegrator(const char *name)
set 1D integrator name
Definition: IntegratorOptions.cxx:208
x
Double_t x[n]
Definition: legend1.C:17
cos
double cos(double)
ROOT::Math::GaussLegendreIntegrator::fNum
int fNum
Definition: GaussLegendreIntegrator.h:128
b
#define b(i)
Definition: RSha256.hxx:118
MATH_WARN_MSG
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:80
ROOT::Math::BaseIntegratorOptions::SetWKSize
void SetWKSize(unsigned int size)
set workspace size
Definition: IntegratorOptions.h:93
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Math::GaussIntegrator::fEpsRel
double fEpsRel
Definition: GaussIntegrator.h:237
ROOT::Math::GaussLegendreIntegrator::fW
double * fW
Definition: GaussLegendreIntegrator.h:130
ROOT::Math::GaussLegendreIntegrator::DoIntegral
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrugate method.
Definition: GaussLegendreIntegrator.cxx:70
ROOT::Math::BaseIntegratorOptions::RelTolerance
double RelTolerance() const
absolute tolerance
Definition: IntegratorOptions.h:74
GaussLegendreIntegrator.h
ROOT::Math::GaussLegendreIntegrator::Options
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
Definition: GaussLegendreIntegrator.cxx:167
Error.h
a
auto * a
Definition: textangle.C:12
ROOT::Math::GaussLegendreIntegrator::SetRelTolerance
virtual void SetRelTolerance(double)
Set the desired relative Error.
Definition: GaussLegendreIntegrator.cxx:96
ROOT::Math::GaussLegendreIntegrator::GaussLegendreIntegrator
GaussLegendreIntegrator(int num=10, double eps=1e-12)
Basic contructor of GaussLegendreIntegrator.
Definition: GaussLegendreIntegrator.cxx:33
ROOT::Math::GaussLegendreIntegrator::~GaussLegendreIntegrator
virtual ~GaussLegendreIntegrator()
Default Destructor.
Definition: GaussLegendreIntegrator.cxx:44
ROOT::Math::GaussIntegrator::fUsedOnce
bool fUsedOnce
Definition: GaussIntegrator.h:239
ROOT::Math::IntegratorOneDimOptions
Numerical one dimensional integration options.
Definition: IntegratorOptions.h:123
ROOT::Math::BaseIntegratorOptions::SetRelTolerance
void SetRelTolerance(double tol)
set the relative tolerance
Definition: IntegratorOptions.h:90
IntegratorOptions.h
ROOT::Math::IBaseFunctionOneDim
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
ROOT::Math::BaseIntegratorOptions::SetAbsTolerance
void SetAbsTolerance(double tol)
non-static methods for setting options
Definition: IntegratorOptions.h:87
ROOT::Math::GaussLegendreIntegrator::SetAbsTolerance
virtual void SetAbsTolerance(double)
This method is not implemented.
Definition: GaussLegendreIntegrator.cxx:103
MATH_WARN_MSGVAL
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:105
ROOT::Math::GaussLegendreIntegrator::fX
double * fX
Definition: GaussLegendreIntegrator.h:129
ROOT::Math::GaussLegendreIntegrator::CalcGaussLegendreSamplingPoints
void CalcGaussLegendreSamplingPoints()
Type: unsafe but fast interface filling the arrays x and w (static method)
Definition: GaussLegendreIntegrator.cxx:108
ROOT::Math::IntegratorOneDimOptions::NPoints
unsigned int NPoints() const
maximum number of function calls
Definition: IntegratorOptions.h:153
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Math