Biblioteka C ++ do integracji numerycznej (kwadratura)

10

Mam własną małą procedurę integracji numerycznej (kwadraturę), która jest adaptacją C ++ programu ALGOL opublikowanego przez Bulirsch & Stoer w 1967 r. (Numerische Mathematik, 9, 271-278).

Chciałbym uaktualnić do bardziej nowoczesnego (adaptacyjnego) algorytmu i zastanawiać się, czy istnieją (bezpłatne) biblioteki C ++, które to zapewniają. Wyglądałem jak GSL (który jest C), ale ma straszne API (chociaż cyfry mogą być dobre). Czy jest coś jeszcze?

Przydatny interfejs API wyglądałby tak:

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);
Walter
źródło
7
Nawiasem mówiąc, przekonasz się, że wiele najlepszych implementacji w informatyce ma „złe” API po prostu dlatego, że były rozwijane przez dziesięciolecia, a nie miesiące czy lata innego oprogramowania. Myślę, że byłoby do przyjęcia i prawdopodobnie bardzo przydatne, aby napisać otoki API i wewnętrznie wywołać mniej czysty API. Daje to zaletę dobrego interfejsu API w kodach podstawowych, a także pozwala łatwo przełączać się między różnymi bibliotekami kwadraturowymi, przepisując tylko jedną funkcję.
Godric Seer,
1
@GodricSeer Gdyby to było takie proste, zrobiłbym to. Jednak tak nie jest. Interfejs API GSL wymaga wstępnie przydzielonego bufora, z którego prawdopodobnie nic nie jest używane, ale który potencjalnie może być zbyt mały (wymagający kolejnego wywołania z większą ilością pamięci). Prawidłowa implementacja byłaby rekurencyjna, nie wymagałaby alokacji, trzymałaby wszystkie dane na stosie i zapewniała czysty interfejs API.
Walter,
1
@GodricSeer Kolejnym poważnym problemem związanym z interfejsem API GSL jest to, że akceptuje tylko funkcje bez stanu (ponieważ używa prostego wskaźnika funkcji). Generowanie API wątkowego dla funkcji z tym stanem jest z konieczności nieefektywne.
Walter,
2
Zgadzam się z Godric Seer, najlepszym rozwiązaniem jest napisanie opakowania. Nie sądzę, że to prawda, że ​​„GSL akceptuje tylko funkcje bez stanu”: tutaj w dokumentacji napisano, że a gsl_functionjest wskaźnikiem funkcji wraz z pewnym nieprzezroczystym wskaźnikiem danych, który może zawierać twój stan. Po drugie, istnieją pewne obawy związane z wydajnością związane z (ponownym) przydzielaniem dowolnie dużych buforów roboczych, tak że część ma co najmniej pewne uzasadnione uzasadnienie.
Kirill,
1
Kolejny komentarz na temat wstępnie przydzielonego bufora GSL. Rozmiar obszaru roboczego jest zdefiniowany w kategoriach maksymalnej liczby interwałów - ponieważ chcesz, aby procedura kwadratury i tak zawiodła, jeśli zajmie zbyt wiele adaptacyjnych bisection, po prostu ustaw rozmiar obszaru roboczego na pewien górny limit liczby bisection. Kiedy mówimy o „właściwej” implementacji, GSL robi tutaj „właściwą” rzecz, przecina interwał z największym obecnie błędem, co oznacza, że ​​musi śledzić wszystkie interwały do ​​tej pory. Jeśli wszystkie dane pozostaną na stosie, może zabraknąć pamięci stosu, nie jest to wcale lepsze.
Kirill,

Odpowiedzi:

3

Spójrz na Odeint . Jest teraz częścią Boost i obejmuje między innymi algorytm Bulirsch-Stoer. Na początek możesz zobaczyć tutaj bardzo prosty przykład.

Zythos
źródło
3
Pierwsze zdanie przeglądu odeint jest następujące: „odeint to biblioteka do rozwiązywania problemów z wartością początkową (IVP) równań różniczkowych zwyczajnych”. O ile mi wiadomo, ta biblioteka nie może być używana do kwadratury znanej funkcji. Czy masz przykład, w którym wykorzystano go do kwadratury?
Bill Greene,
1
Wydaje mi się (sam nie korzystam z biblioteki), że nie zawiera ona algorytmów dla kwadratur, takich jak w kwadraturze Newtona-Cotesa, Romberga lub Gaussa, ale biorąc pod uwagę, że w pytaniu wspomniano o metodzie Gragga-Bulirscha-Stoera, pomyślałem, że to problem była integracja ODE.
Zythos
2

MFEM [1] ma łatwe w użyciu funkcje kwadraturowe (zarówno dla elementów powierzchniowych, jak i objętościowych). Byliśmy w stanie wykorzystać je do różnych zadań.

[1] http://mfem.org/

BrunoLevy
źródło
2

Możesz łatwo napisać cienkie opakowanie C ++ wokół funkcji kwadraturowych GSL. Poniższe potrzeby wymagają 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';
}

Wynik

\int_0^1 x^2 dx = 0.333333
\int_1^\infty x^{-2} dx = 1
\int_{-\infty}^\infty \exp(-x^2) dx = 1.77245
Henri Menke
źródło
1

Odniosłem sukces w bibliotece Cubature (choć jest napisana w C). Ma on na celu integrację wielowymiarową przy stosunkowo małej liczbie wymiarów.

Biblioteka HIntLib jest napisany w języku C ++ i ma rutyny Adaptive kwadraturze (kubaturowego).

Juan M. Bello-Rivas
źródło