10
$\begingroup$

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);
$\endgroup$
  • 6
    $\begingroup$ 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. $\endgroup$ – Godric Seer Sep 21 '15 at 13:34
  • 1
    $\begingroup$ @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. $\endgroup$ – Walter Sep 21 '15 at 14:21
  • 1
    $\begingroup$ @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. $\endgroup$ – Walter Sep 21 '15 at 14:51
  • 2
    $\begingroup$ 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. $\endgroup$ – Kirill Sep 21 '15 at 21:30
  • 1
    $\begingroup$ 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. $\endgroup$ – Kirill Sep 21 '15 at 21:43
3
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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? $\endgroup$ – Bill Greene Aug 1 '16 at 11:59
  • 1
    $\begingroup$ 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. $\endgroup$ – Zythos Aug 1 '16 at 14:11
2
$\begingroup$

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

[1] http://mfem.org/

$\endgroup$
2
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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).

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.