C ++ библиотека для числовой интеграции (квадратура)

10

У меня есть своя небольшая подпрограмма для численного интегрирования (квадратура), которая представляет собой адаптацию C ++ программы ALGOL, опубликованной Bulirsch & Stoer в 1967 году (Numerische Mathematik, 9, 271-278).

Я хотел бы перейти на более современный (адаптивный) алгоритм и задаться вопросом, есть ли какие-либо (бесплатные) библиотеки C ++, которые предоставляют такие. Я выглядел как GSL (это C), но это идет с ужасным API (хотя цифры могут быть хорошими). Что-нибудь еще?

Полезный API будет выглядеть так:

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);
Вальтер
источник
7
Кроме того, вы обнаружите, что многие из лучших реализаций в вычислительной науке имеют «плохие» API просто потому, что они разрабатывались десятилетиями, а не месяцами или годами другого программного обеспечения. Я думаю, что было бы приемлемо и, вероятно, очень полезно для вас написать API-оболочку и внутренне вызвать менее чистый API. Это дает вам преимущество хорошего API в ваших основных кодах, а также позволяет легко переключаться между различными квадратурными библиотеками, переписывая только одну функцию.
Годрик Провидец
1
@GodricSeer Если бы это было так просто, я бы к этому. Однако это не так. Для GSL API требуется предварительно выделенный буфер, из которого, возможно, ничего не используется, но который потенциально может быть слишком маленьким (требующий другого вызова с большим объемом памяти). Правильная реализация будет рекурсивной, не требует выделения, хранит все данные в стеке и обеспечивает чистый API.
Уолтер
1
@GodricSeer Еще одна серьезная проблема с GSL API заключается в том, что он принимает только функции без состояния (поскольку использует простой указатель на функцию). Генерирование API, обеспечивающего безопасность потоков, для функций с состоянием из этого обязательно неэффективно.
Уолтер
2
Я согласен с Годриком Провидцем, написание обертки - лучший вариант. Я не думаю, что это правильно, что «GSL принимает только функции без состояния»: здесь, в документе, говорится, что gsl_functionэто указатель на функцию вместе с некоторым непрозрачным указателем на данные, который может содержать ваше состояние. Во-вторых, существуют некоторые проблемы с эффективностью (пере) распределения произвольно больших рабочих буферов, так что эта часть имеет хотя бы какое-то обоснованное обоснование.
Кирилл
1
Еще один комментарий о предварительно выделенном буфере GSL. Размер рабочего пространства определяется в терминах максимального количества интервалов - поскольку вы хотите, чтобы квадратурная подпрограмма не работала в любом случае, если требуется слишком много адаптивных делений, просто установите для размера рабочего пространства верхний предел числа делений. Когда вы говорите о «правильной» реализации, GSL делает здесь «правильную» вещь, она делит пополам интервал с самой большой на данный момент ошибкой, что означает, что он должен отслеживать все интервалы до сих пор. Если вы храните все данные в стеке, вы можете использовать нехватку стековой памяти, это не совсем лучше.
Кирилл

Ответы:

3

Посмотрите на Odeint . Теперь он является частью Boost и включает в себя алгоритм Bulirsch-Stoer среди других. Для начала вы можете увидеть здесь очень простой пример.

Zythos
источник
3
Первое предложение обзора для odeint: «odeint - это библиотека для решения начальных задач (IVP) обыкновенных дифференциальных уравнений». Насколько мне известно, эту библиотеку нельзя использовать для квадратуры известной функции. У вас есть пример, где он использовался для квадратуры?
Билл Грин
1
Я думаю (я не использую библиотеку сам), что она не включает алгоритмы для квадратур, таких как в квадратурах Ньютона, Ромберга или Гаусса, но учитывая, что в вопросе упоминается метод Грэгга-Булирша-Стоера, я думал, что проблема под рукой была интеграция ODE.
Zythos
2

MFEM [1] имеет простые в использовании квадратурные функции (как для поверхностных, так и для объемных элементов). Мы смогли использовать их для различных задач.

[1] http://mfem.org/

BrunoLevy
источник
2

Вы можете легко написать тонкую оболочку C ++ вокруг квадратурных функций GSL. Следующие требования 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';
}

Вывод

\int_0^1 x^2 dx = 0.333333
\int_1^\infty x^{-2} dx = 1
\int_{-\infty}^\infty \exp(-x^2) dx = 1.77245
Анри Менке
источник
1

У меня был успех с библиотекой Cubature (хотя она написана на C). Он нацелен на многомерную интеграцию с относительно небольшим количеством измерений.

Библиотека HIntLib написана на C ++ и содержит подпрограммы для адаптивной квадратуры (кубатуры).

Хуан М. Белло-Ривас
источник
1

Проверьте https://github.com/tbs1980/NumericIntegration . Он основан на QUADPACK (на котором также основан GSL) и обладает некоторыми изящными современными функциями, например, основанными на Eigen, поддержкой мультиточности.

Джефф
источник