Biblioteca de C ++ para la integración numérica (cuadratura)

Tengo mi propia subrutina para la integración numérica (cuadratura), que es una adaptación en C ++ de un programa ALGOL publicado por Bulirsch & Stoer en 1967 (Numerische Mathematik, 9, 271-278).

Me gustaría actualizar a un algoritmo más moderno (adaptativo) y preguntarme si hay bibliotecas C ++ (gratuitas) que lo proporcionen. Tenía un aspecto como GSL (que es C), pero eso viene con una API horrible (aunque los números pueden ser buenos). ¿Hay algo mas?

Una API útil se vería así:

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

Echa un vistazo a Odeint . Ahora es parte de Boost e incluye el algoritmo Bulirsch-Stoer entre otros. Para comenzar puedes ver aquí un ejemplo muy simple.

MFEM [1] tiene funciones de cuadratura fáciles de usar (tanto para elementos superficiales como volumétricos). Pudimos usarlos para diversas tareas.

[1] http://mfem.org/

Puede escribir fácilmente un envoltorio delgado de C ++ alrededor de las funciones de cuadratura GSL. Lo siguiente necesita 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';
}

Salida

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

He tenido éxito con la biblioteca Cubature (aunque está escrita en C). Está dirigido a la integración multidimensional con un número relativamente bajo de dimensiones.

La biblioteca HIntLib está escrita en C ++ y tiene rutinas para la cuadratura adaptativa (cubature).

Echa un vistazo a https://github.com/tbs1980/NumericalIntegration . Se basa en QUADPACK (en el que también se basa GSL), y tiene algunas características modernas, por ejemplo, basadas en Eigen, soporte de multiprecisión.