C ++ - Bibliothek zur numerischen Integration (Quadratur)

10

Ich habe meine eigene kleine Subroutine für die numerische Integration (Quadratur), eine C ++ - Adaption eines ALGOL-Programms, das 1967 von Bulirsch & Stoer veröffentlicht wurde (Numerische Mathematik, 9, 271-278).

Ich möchte auf einen moderneren (adaptiven) Algorithmus upgraden und mich fragen, ob es (kostenlose) C ++ - Bibliotheken gibt, die solche anbieten. Ich habe nach GSL gesucht (das ist C), aber das kommt mit einer schrecklichen API (obwohl die Zahlen gut sein können). Gibt es noch etwas?

Eine nützliche API würde folgendermaßen aussehen:

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
quelle
7
Abgesehen davon werden Sie feststellen, dass viele der besten Implementierungen in der Computerwissenschaft "schlechte" APIs haben, einfach weil sie über Jahrzehnte entwickelt wurden und nicht über Monate oder Jahre anderer Software. Ich denke, es wäre akzeptabel und wahrscheinlich sehr nützlich für Sie, eine Wrapper-API zu schreiben und die weniger saubere API intern aufzurufen. Dies gibt Ihnen den Vorteil einer netten API in Ihren Primärcodes und ermöglicht es Ihnen, einfach zwischen verschiedenen Quadraturbibliotheken zu wechseln, indem Sie nur eine einzige Funktion neu schreiben.
Godric Seer
1
@GodricSeer Wenn es so einfach wäre, würde ich dazu. Dies ist jedoch nicht der Fall. Die GSL-API erfordert einen vorab zugewiesenen Puffer, von dem möglicherweise nichts verwendet wird, der jedoch möglicherweise zu klein ist (was einen weiteren Aufruf mit mehr Speicher erfordert). Eine ordnungsgemäße Implementierung wäre rekursiv, erfordert keine Zuordnung, behält alle Daten auf dem Stapel und stellt eine saubere API bereit.
Walter
1
@GodricSeer Ein weiteres ernstes Problem mit der GSL-API besteht darin, dass nur Funktionen ohne Status akzeptiert werden (da ein einfacher Funktionszeiger verwendet wird). Das Generieren einer threadsicheren API für Funktionen mit Status daraus ist notwendigerweise ineffizient.
Walter
2
Ich stimme Godric Seer zu, das Schreiben eines Wrappers ist die beste Option. Ich denke nicht, dass es richtig ist, dass "GSL nur Funktionen ohne Status akzeptiert": Hier in den Dokumenten heißt es, dass a gsl_functionein Funktionszeiger zusammen mit einem undurchsichtigen Datenzeiger ist, der Ihren Status enthalten kann. Zweitens gibt es einige Effizienzbedenken hinsichtlich der (Neu-) Zuweisung beliebig großer Arbeitspuffer, so dass dieser Teil zumindest eine gültige Rechtfertigung dafür hat.
Kirill
1
Ein weiterer Kommentar zum vorab zugewiesenen Puffer von GSL. Die Größe des Arbeitsbereichs wird anhand der maximalen Anzahl von Intervallen definiert. Da die Quadraturroutine ohnehin fehlschlagen soll, wenn zu viele adaptive Halbierungen erforderlich sind, legen Sie einfach die Größe des Arbeitsbereichs auf eine Obergrenze für die Anzahl der Halbierungen fest. Wenn Sie von einer "richtigen" Implementierung sprechen, macht GSL hier das "Richtige", es halbiert das Intervall mit dem derzeit größten Fehler, was bedeutet, dass es alle bisherigen Intervalle verfolgen muss. Wenn Sie alle Daten auf dem Stapel behalten, geht Ihnen möglicherweise der Stapelspeicher aus. Dies ist nicht wirklich besser.
Kirill

Antworten:

3

Schauen Sie sich Odeint an . Es ist jetzt Teil von Boost und enthält unter anderem den Bulirsch-Stoer-Algorithmus. Zu Beginn sehen Sie hier ein sehr einfaches Beispiel.

Zythos
quelle
3
Der erste Satz der Übersicht für Odeint lautet: "Odeint ist eine Bibliothek zur Lösung von Anfangswertproblemen (IVP) gewöhnlicher Differentialgleichungen." Soweit mir bekannt ist, kann diese Bibliothek nicht für die Quadratur einer bekannten Funktion verwendet werden. Haben Sie ein Beispiel, in dem es für die Quadratur verwendet wurde?
Bill Greene
1
Ich denke (ich benutze die Bibliothek nicht selbst), dass sie keine Algorithmen für Quadraturen wie in Newton-Cotes, Romberg oder Gaußsche Quadratur enthält, aber da die Frage die Gragg-Bulirsch-Stoer-Methode erwähnte, dachte ich, dass das Problem vorliegt war eine ODE-Integration.
Zythos
2

MFEM [1] verfügt über benutzerfreundliche Quadraturfunktionen (sowohl für Oberflächen- als auch für Volumenelemente). Wir konnten sie für verschiedene Aufgaben einsetzen.

[1] http://mfem.org/

BrunoLevy
quelle
2

Sie können leicht einen dünnen C ++ - Wrapper um die GSL-Quadraturfunktionen schreiben. Folgendes benötigt 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';
}

Ausgabe

\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
quelle
1

Ich hatte Erfolg mit der Cubature-Bibliothek (sie ist jedoch in C geschrieben). Es zielt auf eine mehrdimensionale Integration mit einer relativ geringen Anzahl von Dimensionen ab.

Die HIntLib-Bibliothek ist in C ++ geschrieben und verfügt über Routinen für die adaptive Quadratur (Kubatur).

Juan M. Bello-Rivas
quelle