# C++ library for numerical intergration (quadrature)

I have my own little subroutine for numerical integration (quadrature), which is a C++ adaptation of an ALGOL program published by Bulirsch & Stoer in 1967 (Numerische Mathematik, 9, 271-278).

I would like to upgrade to a more modern (adaptive) algorithm and wonder whether there are any (free) C++ libraries that provide such. I had a look as GSL (which is C), but that comes with a horrible API (though the numerics may be good). Is there anything else?

A useful API would look like this:

double quadrature(double lower_integration_limit,
double upper_integration_limit,
std::function<double(double)> const&func,
double desired_error_bound_relative=1.e-12,
double desired_error_bound_absolute=0,
double*error_estimate=nullptr);

• Just as an aside, you will find that many of the best implementations in computational science have "bad" API's simply because they have been developed over decades, rather than the months or years of other software. I think it would be acceptable, and likely very useful for you to write a wrapper API, and internally call the less clean API. This gives you the advantage of a nice API in your primary codes, and also allows you to easily switch between different quadrature libraries with only rewriting a single function. Sep 21, 2015 at 13:34
• @GodricSeer If it were that simple, I would to that. However, it is not. The GSL API requires a pre-allocated buffer, of which possibly nothing is used, but which potentially may be too small (requiring another call with more memory). A proper implementation would be recursive, need no allocation, keep all data on the stack, and provide a clean API. Sep 21, 2015 at 14:21
• @GodricSeer Another serious problem with the GSL API is that it only accepts functions without state (because it uses a simple function pointer). Generating a threadsafe API for functions with state from this is necessarily inefficient. Sep 21, 2015 at 14:51
• I agree with Godric Seer, writing a wrapper is the best option. I don't think it's correct that "GSL only accepts functions without state": here in the docs it says a gsl_function is a function pointer together with some opaque data pointer, which can contain your state. Second, there are some efficiency concerns about (re-)allocating arbitrarily-large work buffers, so that part has at least some valid justification to it. Sep 21, 2015 at 21:30
• Another comment about GSL's pre-allocated buffer. The size of the workspace is defined in terms of maximum number of intervals - since you want the quadrature routine to fail anyway if it takes too many adaptive bisections, just set the size of the workspace to some upper limit on the number of bisections. When you talk about a "proper" implementation, GSL does the "right" thing here, it bisects the interval with the currently largest error, meaning it has to keep track of all intervals so far. If you keep all data on the stack, you might run of out stack memory, it's not really better. Sep 21, 2015 at 21:43

## 5 Answers

Take a look at Odeint. It is now part of Boost and it includes the Bulirsch-Stoer algorithm among others. To start you can see here a very simple example.

• The first sentence of the overview for odeint is: "odeint is a library for solving initial value problems (IVP) of ordinary differential equations." As far as I am aware, this library cannot be used for quadrature of a known function. Do you have an example where it has been used for quadrature? Aug 1, 2016 at 11:59
• I think (I do not use the library myself) that it does not include algorithms for quadratures such as in Newton-Cotes, Romberg or Gaussian quadrature but given that the question mentioned the Gragg-Bulirsch-Stoer method I thought that the problem at hand was an ODE integration. Aug 1, 2016 at 14:11

You can easily write a thin C++ wrapper around the GSL quadrature functions. The following needs C++11.

#include <iostream>
#include <cmath>

#include <functional>
#include <memory>
#include <utility>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

template < typename F >
class gsl_quad
{
F f;
int limit;
std::unique_ptr < gsl_integration_workspace,
std::function < void(gsl_integration_workspace*) >
> workspace;

static double gsl_wrapper(double x, void * p)
{
gsl_quad * t = reinterpret_cast<gsl_quad*>(p);
return t->f(x);
}

public:
gsl_quad(F f, int limit)
: f(f)
, limit(limit)
, workspace(gsl_integration_workspace_alloc(limit), gsl_integration_workspace_free)
{}

double integrate(double min, double max, double epsabs, double epsrel)
{
gsl_function gsl_f;
gsl_f.function = &gsl_wrapper;
gsl_f.params = this;

double result, error;
if ( !std::isinf(min) && !std::isinf(max) )
{
gsl_integration_qags ( &gsl_f, min, max,
epsabs, epsrel, limit,
workspace.get(), &result, &error );
}
else if ( std::isinf(min) && !std::isinf(max) )
{
gsl_integration_qagil( &gsl_f, max,
epsabs, epsrel, limit,
workspace.get(), &result, &error );
}
else if ( !std::isinf(min) && std::isinf(max) )
{
gsl_integration_qagiu( &gsl_f, min,
epsabs, epsrel, limit,
workspace.get(), &result, &error );
}
else
{
gsl_integration_qagi ( &gsl_f,
epsabs, epsrel, limit,
workspace.get(), &result, &error );
}

return result;
}
};

template < typename F >
double quad(F func,
std::pair<double,double> const& range,
double epsabs = 1.49e-8, double epsrel = 1.49e-8,
int limit = 50)
{
return gsl_quad<F>(func, limit).integrate(range.first, range.second, epsabs, epsrel);
}

int main()
{
std::cout << "\\int_0^1 x^2 dx = "
<< quad([](double x) { return x*x; }, {0,1}) << '\n'
<< "\\int_1^\\infty x^{-2} dx = "
<< quad([](double x) { return 1/(x*x); }, {1,INFINITY}) << '\n'
<< "\\int_{-\\infty}^\\infty \\exp(-x^2) dx = "
<< quad([](double x) { return std::exp(-x*x); }, {-INFINITY,INFINITY}) << '\n';
}


Output

\int_0^1 x^2 dx = 0.333333
\int_1^\infty x^{-2} dx = 1
\int_{-\infty}^\infty \exp(-x^2) dx = 1.77245


MFEM [1] has easy-to-use quadrature functions (both for surfacic and volumetric elements). We were able to use them for various tasks.

I have had success with the Cubature library (it is written in C, though). It is aimed at multidimensional integration with a relatively low number of dimensions.

The HIntLib library is written in C++ and it has routines for adaptive quadrature (cubature).

Check out https://github.com/tbs1980/NumericalIntegration. It is based on QUADPACK (which GSL is also based on), and has some neat modern features, e.g. based on Eigen, multiprecision support.