 ROOT   6.19/01 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
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;
73
74  double result = 0.0;
75  for (int i=0; i<fNum; i++)
76  {
77  xx = 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 == 0 )
107  delete [] fX;
108
109  if ( fW == 0 )
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
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
auto * m
Definition: textangle.C:8
VSD Structures.
Definition: StringConv.hxx:21
void SetWKSize(unsigned int size)
set workspace size
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:104
GaussLegendreIntegrator(int num=10, double eps=1e-12)
Basic contructor of GaussLegendreIntegrator.
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...
void CalcGaussLegendreSamplingPoints()
Type: unsafe but fast interface filling the arrays x and w (static method)
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:79
double RelTolerance() const
absolute tolerance
double cos(double)
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrugate method.
Double_t x[n]
Definition: legend1.C:17
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options (should be re-implemented by derived classes -if more options than tolerance exist ...
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual ~GaussLegendreIntegrator()
Default Destructor.
virtual void SetRelTolerance(double)
Set the desired relative Error.
auto * a
Definition: textangle.C:12
Numerical one dimensional integration options.
User class for performing function integration.
unsigned int NPoints() const
maximum number of function calls
void SetNumberPoints(int num)
Set the number of points used in the calculation of the integral.
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual void SetAbsTolerance(double)
This method is not implemented.
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
void SetIntegrator(const char *name)
set 1D integrator name
void SetRelTolerance(double tol)
set the relative tolerance
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...
void SetAbsTolerance(double tol)
non-static methods for setting options