Function Integration

Class structure for Numerical Integration

The algorithms provided by ROOT for numerical integration are implemented following the hierarchy shown in the next image. There are three base interfaces from which every method should derive from. VirtualIntegrator defines the most basic functionality while the specific behaviours for one or multiple dimensions are implemented in the other two interfaces. Each of these interfaces offer a method to set the function to be integrated, that must be of the most basic type. The only difference between the interfaces resides in the dimensionality of that function and some overloading that will be seen afterwards for the one dimensional one.

The rest of the classes shown in the diagram are the specialized classes provided. Each one implements a different method that will be explained in detail. It is important to notice that the two grayed classes (the one which name starts by GSL) are part of the MathMore library. Next section shows in more detail the differences between the implementations:

Using an existing Integrator

Using the class directly

The next is a list of the classes in ROOT that perform numerical integration. Every method implemented inside de class is briefly commented followed by an example of its use when directly created.

• GaussIntegratorOneDim: It uses the most basic gaussian integrator with a general approximation of 12 points. Example:

  #include "TF1.h"
  #include "Math/WrappedTF1.h"
  #include "Math/GaussIntegrator.h"
   
  int NumericalIntegration()
  {
   
     TF1 f("Sin Function", "sin(x)", 0, TMath::Pi());
   
     ROOT::Math::WrappedTF1 wf1(f);
   
     ROOT::Math::GaussIntegrator ig;
   
     ig.SetFunction(wf1, false);
     ig.SetRelTolerance(0.001);
   
     cout << ig.Integral(0, TMath::PiOver2()) << endl;
   
     return 0;
  }

• GaussLegendreIntegrator: This class implementes the Gauss-Legendre quadrature formulas. This sort of numerical methods requieres that the user specifies the number of intermediate function points used in the calculation of the integral. It will automatically determine the coordinates and weights of such points before performing the integration. Example:

  #include "TF1.h"
  #include "Math/WrappedTF1.h"
  #include "Math/GaussLegendreIntegrator.h"
   
  int NumericalIntegration()
  {
   
     // Create the function and wrap it
     TF1 f("Sin Function", "sin(x)", 0, TMath::Pi());
     ROOT::Math::WrappedTF1 wf1(f);
   
     // Create the Integrator
     ROOT::Math::GaussLegendreIntegrator ig;
   
     // Set parameters of the integration
     ig.SetFunction(wf1, false);
     ig.SetRelTolerance(0.001);
     ig.SetNumberPoints(40);
   
     cout << ig.Integral(0, TMath::PiOver2()) << endl;
   
     return 0;
  }

• GSLIntegrator: This is a wrapper for the QUADPACK integrator implemented in the GSL library. It supports several integration methods that can be chosen in construction time. The default type is adaptive integration with singularity applying a Gauss-Kronrod 21-point integration rule. For a detail description of the GSL methods visit the GSLreference guide

  This class implements the best algorithms for numerical integration for one dimensional functions. We encourage the use it as the main option, bearing in mind that it uses code from the GSL library, with the respective license restriction it has. Example:

  #include "TF1.h"
  #include "Math/WrappedTF1.h"
  #include "Math/GSLIntegrator.h"
   
  int NumericalIntegration()
  {
   
     // Create the function and wrap it
     TF1 f("Sin Function", "sin(x)", 0, TMath::Pi());
     ROOT::Math::WrappedTF1 wf1(f);
   
     // Create the Integrator
     ROOT::Math::GSLIntegrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVE);
   
     // Set parameters of the integration
     ig.SetFunction(wf1, false);
     ig.SetRelTolerance(0.001);
   
     cout << ig.Integral(0, TMath::PiOver2()) << endl;
   
     return 0;
  }

• AdaptiveIntegratorMultiDim: This class implements an adaptive quadrature integration method for multi dimensional functions. It is described in [Genz]. Example:

  #include "TF2.h"
  #include "Math/WrappedMultiTF1.h"
  #include "Math/AdaptiveIntegratorMultiDim.h"
   
  int NumericalIntegration()
  {
   
     // Create the function and wrap it
     TF2 f("Sin Function", "sin(x) + y", 0, TMath::Pi(), 0, 1);
     ROOT::Math::WrappedMultiTF1 wf1(f);
   
     // Create the Integrator
     ROOT::Math::AdaptiveIntegratorMultiDim ig;
   
     // Set parameters of the integration
     ig.SetFunction(wf1);
     ig.SetRelTolerance(0.001);
   
     double xmin[] = {0,0};
     double xmax[] = {TMath::PiOver2(), 3};
   
     cout << ig.Integral(xmin, xmax) << endl;
   
     return 0;
  }

• GSLMCIntegrator: It is a class for performing numerical integration of a multidimensional function. It uses the numerical integration algorithms of GSL, which reimplements the algorithms used in the QUADPACK, a numerical integration package written in Fortran. Plain MC, MISER and VEGAS integration algorithms are supported for integration over finite (hypercubic) ranges. For a detail description of the GSL methods visit the GSL reference guide

  Again, this class implements the best algoritms for numerical integration for multi dimensional functions. Beware of the license restrictions of the GSL library if you are using this class. Example:

  int NumericalIntegration()
  {
   
     // Create the function and wrap it
     TF2 f("Sin Function", "sin(x) + y", 0, TMath::Pi(), 0 ,1);
     ROOT::Math::WrappedMultiTF1 wf1(f);
   
     // Create the Integrator
     ROOT::Math::GSLMCIntegrator ig(ROOT::Math::IntegrationMultiDim::kVEGAS);
   
     // Set parameters of the integration
     ig.SetFunction(wf1);
     ig.SetRelTolerance(0.001);
   
     double xmin[] = {0,0};
     double xmax[] = {TMath::PiOver2(), 3};
   
     cout << ig.Integral(xmin, xmax) << endl;
   
     return 0;
  }

Using the Plug-In Mananger

There is a second way of using any of the classes mentioned before. Through the Plug-In Manager, and in particular through the classes ROOT::Math::Integrator and ROOT::Math::IntegratorMultiDim, the user can initialize any of the classes we have seen before without dealing with them directly. These two classes provide an interface similar to the ones in VirtualIntegratorOneDim and VirtualIntegratorMultiDim, but with the possibility to choose with the constructor, which method will be used to perform the integration. Example:

#include "Math/Integrator.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/AllIntegrationTypes.h"
#include "Math/Functor.h"
#include "Math/GaussIntegrator.h"
 
#include <cmath>
 
const double ERRORLIMIT = 1E-3;
 
double f(double x) { 
   return x; 
} 
 
double f2(const double * x) { 
   return x[0] + x[1]; 
} 
 
 
int testIntegration1D() { 
 
   const double RESULT = 0.5;
   int status = 0;
 
   ROOT::Math::Functor1D wf(&f);
   ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR); 
   ig.SetFunction(wf);
   double val = ig.Integral(0,1);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   ROOT::Math::Integrator ig2(ROOT::Math::IntegrationOneDim::kNONADAPTIVE); 
   ig2.SetFunction(wf);
   val = ig2.Integral(0,1);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   ROOT::Math::Integrator ig3(wf, ROOT::Math::IntegrationOneDim::kADAPTIVE); 
   val = ig3.Integral(0,1);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   //ROOT::Math::GaussIntegrator ig4;
   ROOT::Math::Integrator ig4(ROOT::Math::IntegrationOneDim::kGAUSS); 
   ig4.SetFunction(wf);
   val = ig4.Integral(0,1);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   return status;
}
 
int testIntegrationMultiDim() { 
 
   const double RESULT = 1.0;
   int status = 0;
 
   ROOT::Math::Functor wf(&f2,2);
   double a[2] = {0,0};
   double b[2] = {1,1};
 
   ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE); 
   ig.SetFunction(wf);
   double val = ig.Integral(a,b);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   ROOT::Math::IntegratorMultiDim ig2(ROOT::Math::IntegrationMultiDim::kVEGAS); 
   ig2.SetFunction(wf);
   val = ig2.Integral(a,b);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   ROOT::Math::IntegratorMultiDim ig3(wf,ROOT::Math::IntegrationMultiDim::kPLAIN); 
   val = ig3.Integral(a,b);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   ROOT::Math::IntegratorMultiDim ig4(wf,ROOT::Math::IntegrationMultiDim::kMISER); 
   val = ig4.Integral(a,b);
   std::cout << "integral result is " << val << std::endl;
   status += std::fabs(val-RESULT) > ERRORLIMIT;
 
   return status;
}
 
int  main() { 
   int status = 0;
 
   status += testIntegration1D();
   status += testIntegrationMultiDim();
 
   return status;
}

How to implement your own Integrator

Explain where to Derive From. Show the code of some real ones done in MathCore!

Show also how to integrate the class with the Plug-In Manager!

Pending of confirmation to do!

References

• Genz, A.A. Malik, An adaptive algorithm for numerical integration over an N-dimensional rectangular region, J. Comput. Appl. Math. 6 (1980) 295-302.
• GSL Documentation for Numerical Integration.
• Numerical Integration with ROOT