Logo ROOT   6.14/05
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[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 == 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
static double p3(double t, double a, double b, double c, double d)
Namespace for new ROOT classes and functions.
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
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:47
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
static double p2(double t, double a, double b, double c)
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
static double p1(double t, double a, double b)
void SetNumberPoints(int num)
Set the number of points used in the calculation of the integral.
#define MATH_WARN_MSGVAL(loc, str, x)
Definition: Error.h:67
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