Integration

One-dimensional integration

Several classes integrate arbitrary one-dimensional functions

All of these classes are children of o2scl::inte, and the relative and absolute tolerances are stored in o2scl::inte::tol_rel and o2scl::inte::tol_abs, respectively.

There are two competing factors that can slow down an adaptive integration algorithm: (1) each evaluation of the integrand can be numerically expensive, depending on how the function is defined, and (2) the process of subdividing regions and recalculating values is almost always numerically expensive in its own right. For integrands that are very smooth (e.g., analytic functions), a high-order Gauss-Kronrod rule (e.g., 61-point) will achieve the desired error tolerance with relatively few subdivisions. For integrands with discontinuities or singular derivatives, a low-order rule (e.g., 15-point) is often more efficient.

GSL-based integration routines

For the GSL-based integration routines, the variables o2scl::inte::tol_abs and o2scl::inte::tol_rel have the same role as the quantities usually denoted in the GSL integration routines by epsabs and epsrel. In particular, the integration classes attempt to ensure that

\[ |\mathrm{result}-I| \leq \mathrm{Max}(\mathrm{tol\_abs}, \mathrm{tol\_rel}|I|) \]

and returns an error to attempt to ensure that

\[ |\mathrm{result}-I| \leq \mathrm{abserr} \leq \mathrm{Max}(\mathrm{tol\_abs},\mathrm{tol\_rel}|I|) \]

where I is the integral to be evaluated. Even when the corresponding descendant of o2scl::inte::integ() returns success, these inequalities may fail for sufficiently difficult functions. All of the GSL integration routines except for o2scl::inte_qng_gsl use a workspace given in o2scl::inte_workspace_gsl which holds the results of the various subdivisions of the original interval.

The GSL routines were originally based on QUADPACK, which is available at http://www.netlib.org/quadpack .

For adaptive GSL integration classes, the type of Gauss-Kronrod quadrature rule that is used to approximate the integral and estimate the error of a subinterval is set by o2scl::inte_kronrod_gsl::set_rule().

The number of subdivisions of the integration region is limited by the size of the workspace, set in o2scl::inte_kronrod_gsl::set_limit(). The number of subdivisions required for the most recent call to o2scl::inte::integ() or o2scl::inte::integ_err() is given in o2scl::inte::last_iter. This number will always be less than or equal to the workspace size.

Note
The GSL integration routines can sometimes lose precision if the integrand is everywhere much smaller than unity. Some rescaling may be required in these cases.

GSL-based integration error messages

The error messages given by the adaptive GSL integration routines tend to follow a standard form and are documented here. There are several error messages which indicate improper usage and cause the error handler to be called regardless of the value of o2scl::inte::err_nonconv :

There are also convergence errors which will call the error handler unless o2scl::inte::err_nonconv is false (see What is an error?) for more discussion on convergence errors versus fatal errors):

Multi-dimensional integration routines

O2scl reimplements the Cubature library for multi-dimensional integration. The h-adaptive and p-adaptive integration methods are implemented in o2scl::inte_hcubature and o2scl::inte_pcubature . See also the Monte Carlo integration routines in Monte Carlo Integration .

One-dimensional integration example

This example computes the integral $ \int_{-\infty}^{\infty} e^{-x^2} ~dx $ with o2scl::inte_qagi_gsl, the integral $ \int_0^{\infty} e^{-x^2} ~dx $ with o2scl::inte_qagiu_gsl, the integral $ \int_{-\infty}^{0} e^{-x^2} ~dx $ with o2scl::inte_qagil_gsl, and the integral $ \int_0^1 \left[ \sin (2 x) + \frac{1}{2} \right]~dx $ with both o2scl::inte_qag_gsl and o2scl::inte_adapt_cern, and compares the computed results with the exact results.

/* Example: ex_inte.cpp
-------------------------------------------------------------------
An example to demonstrate numerical integration.
*/
#include <cmath>
#include <o2scl/test_mgr.h>
#include <o2scl/constants.h>
#include <o2scl/funct.h>
#include <o2scl/inte_qag_gsl.h>
#include <o2scl/inte_qagi_gsl.h>
#include <o2scl/inte_qagiu_gsl.h>
#include <o2scl/inte_qagil_gsl.h>
#include <o2scl/inte_adapt_cern.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_const;
class cl {
public:
// We'll use this to count the number of function
// evaulations required by the integration routines
int nf;
// A function to be integrated
double integrand(double x) {
nf++;
return exp(-x*x);
}
// Another function to be integrated
double integrand2(double x) {
nf++;
return sin(2.0*x)+0.5;
}
};
int main(void) {
cl acl;
funct11 f1=std::bind(std::mem_fn<double(double)>
(&cl::integrand),&acl,std::placeholders::_1);
funct11 f2=std::bind(std::mem_fn<double(double)>
(&cl::integrand2),&acl,std::placeholders::_1);
// We don't need to specify the function type in the integration
// objects, because we're using the default function type (type
// funct).
// The result and the uncertainty
double res, err;
// An integral from -infinity to +infinity (the limits are ignored)
acl.nf=0;
int ret1=gi.integ_err(f1,0.0,0.0,res,err);
cout << "inte_qagi_gsl: " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Result: " << res << " Uncertainty: " << err << endl;
cout << "Number of iterations: " << gi.last_iter << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(res,sqrt(pi),1.0e-8,"inte 1");
// An integral from 0 to +infinity (the second limit argument is
// ignored in the line below)
acl.nf=0;
gu.integ_err(f1,0.0,0.0,res,err);
cout << "inte_qagiu_gsl: " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Result: " << res << " Uncertainty: " << err << endl;
cout << "Number of iterations: " << gu.last_iter << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(res,sqrt(pi)/2.0,1.0e-8,"inte 2");
// An integral from -infinity to zero (the first limit argument is
// ignored in the line below)
acl.nf=0;
gl.integ_err(f1,0.0,0.0,res,err);
cout << "inte_qagil_gsl: " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Result: " << res << " Uncertainty: " << err << endl;
cout << "Number of iterations: " << gl.last_iter << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(res,sqrt(pi)/2.0,1.0e-8,"inte 3");
// An integral from 0 to 1 with the GSL integrator
acl.nf=0;
g.integ_err(f2,0.0,1.0,res,err);
cout << "inte_qag_gsl: " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Result: " << res << " Uncertainty: " << err << endl;
cout << "Number of iterations: " << g.last_iter << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(res,0.5+sin(1.0)*sin(1.0),1.0e-8,"inte 4");
// The same integral with the CERNLIB integrator
acl.nf=0;
ca.integ_err(f2,0.0,1.0,res,err);
cout << "inte_adapt_cern: " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Result: " << res << " Uncertainty: " << err << endl;
cout << "Number of iterations: " << ca.last_iter << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(res,0.5+sin(1.0)*sin(1.0),1.0e-8,"inte 5");
t.report();
return 0;
}
// End of example

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).