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"
16#include <cmath>
17#include <string.h>
18#include <algorithm>
19
20namespace ROOT {
21namespace 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
51void 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
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
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:79
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:104
#define b(i)
Definition: RSha256.hxx:100
double cos(double)
void SetAbsTolerance(double tol)
non-static methods for setting options
double RelTolerance() const
absolute tolerance
void SetRelTolerance(double tol)
set the relative tolerance
void SetWKSize(unsigned int size)
set workspace size
User class for performing function integration.
virtual double DoIntegral(double a, double b, const IGenFunction *func)
Integration surrugate method.
virtual void SetAbsTolerance(double)
This method is not implemented.
virtual ~GaussLegendreIntegrator()
Default Destructor.
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options (should be re-implemented by derived classes -if more options than tolerance exist
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...
GaussLegendreIntegrator(int num=10, double eps=1e-12)
Basic contructor of GaussLegendreIntegrator.
void CalcGaussLegendreSamplingPoints()
Type: unsafe but fast interface filling the arrays x and w (static method)
virtual void SetRelTolerance(double)
Set the desired relative Error.
void SetNumberPoints(int num)
Set the number of points used in the calculation of the integral.
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Numerical one dimensional integration options.
void SetIntegrator(const char *name)
set 1D integrator name
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...
unsigned int NPoints() const
maximum number of function calls
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12