Численное интегрирование в C++

Вы можете использовать Научную библиотеку GNU, которая поддерживает множество функций "численного анализа", включая интеграцию .

Очень простой пример интеграции из руководства всего несколько строк кода:

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
  double alpha = *(double *) params;
  return log(alpha*x) / sqrt(x);
}

int
main (void)
{
  double result, error;
  double expected = -4.0;
  double alpha = 1.0;
  gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (1000);

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                        w, &result, &error); 

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("intervals =  %d\n", w->size);

  gsl_integration_workspace_free (w);

  return 0;
}

Насколько я знаю, в стандартной библиотеке нет функций того типа, который вы ищете. Вот одна из реализаций:

Расширенное правило трапеций.

Чтобы фиксированная функция f(x) интегрировалась между фиксированными пределами a и b, можно удвоить количество интервалов в расширенном правиле трапеций без потери преимуществ предыдущей работы. Самая грубая реализация правила трапеций состоит в усреднении функции в ее конечных точках a и b. Первый этап уточнения заключается в добавлении к этому среднему значению функции в средней точке. Второй этап уточнения заключается в добавлении значений в точках 1/4 и 3/4.

Ряд элементарных квадратурных алгоритмов включает в себя добавление последовательных этапов уточнения. Эту функцию удобно инкапсулировать в структуру Quadrature:

struct Quadrature
{  
   //Abstract base class for elementary quadrature algorithms.
   Int n; // Current level of refinement.

   virtual Doub next() = 0;
   //Returns the value of the integral at the nth stage of refinement. 
   //The function next() must be defined in the derived class.
};

Затем структура Trapzd получается из этого следующим образом:

template<class T> 
struct Trapzd: Quadrature
{
    Doub a, b, s; // Limits of integration and current value of integral.
    T &func;

    Trapzd() { };

    // func is function or functor to be integrated between limits: a and b 
    Trapzd(T &funcc, const Doub aa, const Doub bb)   
        : func(funcc), a(aa), b(bb)
    {
        n = 0;
    }

    // Returns the nth stage of refinement of the extended trapezoidal rule. 
    // On the first call (n = 1), the routine returns the crudest estimate  
    // of integral of f x / dx in [a,b]. Subsequent calls set n=2,3,... and
    // improve the accuracy by adding 2n - 2 additional interior points.
    Doub next()
    {
        Doub x, tnm, sum, del;
        Int it, j;
        n++;

        if (n == 1)
        {
            return (s = 0.5 * (b-a) * (func(a) + func(b)));
        } 
        else
        {
            for (it = 1, j = 1; j < n - 1; j++)
            {
                it <<= 1;
            }
            tnm = it;
            // This is the spacing of the points to be added.          
            del = (b - a) / tnm; 
            x = a + 0.5 * del;

            for (sum = 0.0,j = 0; j < it; j++, x += del)
            {
                sum += func(x);
            }
            // This replaces s by its refined value.  
            s = 0.5 * (s + (b - a) * sum / tnm); 
            return s;
        }
    }
};

Применение:

Структуру Trapzd можно использовать несколькими способами. Самый простой и грубый способ — это интегрировать функцию по расширенному правилу трапеций, когда вы заранее знаете, сколько шагов вам нужно. Если вы хотите 2^M + 1, вы можете сделать это с помощью фрагмента:

Ftor func;      // Functor func here has no parameters.
Trapzd<Ftor> s(func, a, b);
for(j = 1 ;j <= m + 1; j++) val = s.next();

с ответом, возвращенным как val. Здесь Ftor — функтор, содержащий интегрируемую функцию.

Бонус:

Гораздо лучше, конечно, уточнять правило трапеций до тех пор, пока не будет достигнута некоторая заданная степень точности. Функция для этого:

template<class T>
Doub qtrap(T &func, const Doub a, const Doub b, const Doub eps = 1.0e-10)
{
    // Returns the integral of the function or functor func from a to b. 
    // The constants EPS can be set to the desired fractional accuracy and    
    // JMAX so that 2 to the power JMAX-1 is the maximum allowed number of   
    // steps. Integration is performed by the trapezoidal rule.

    const Int JMAX = 20;
    Doub s, olds = 0.0; // Initial value of olds is arbitrary.

    Trapzd<T> t(func, a, b);

    for (Int j = 0; j < JMAX; j++) 
    {
        s = t.next();

        if (j > 5) // Avoid spurious early convergence.
        {
            if (abs(s - olds) < eps * abs(olds) || (s == 0.0 && olds == 0.0)) 
            {
                return s;
            }
        }
        olds = s;
    }
    throw("Too many steps in routine qtrap");
}

Типы.

typedef double Doub;    // 64 - bit floating point
typedef int Int;        // 32 - bit signed integer

Возможно, вы захотите проверить Boost Quadrature and Differentiation Library. В частности, они предоставляют версию правила трапеций:

https://www.boost.org/doc/libs/1_69_0/libs/math/doc/html/math_toolkit/trapezoidal.html

Библиотека квадратур/дифференциации очень хорошо написана и совместима с современным C++, поскольку можно просто использовать лямбда-выражение или объект функции для подынтегрального выражения. Я очень быстро с ним разобрался.

Вот пример аппроксимации pi путем аппроксимации интеграла от 4/(1 + x^2) от x = 0 до x = 1, помещая подынтегральное выражение в виде лямбда-выражения.

#include <boost/math/quadrature/trapezoidal.hpp>
#include <iostream>
using boost::math::quadrature::trapezoidal;
using std::cout;
using std::endl;

// Put inside a test function:
auto f = [](double x)
{
    return 4.0 / (1.0 + x * x);
};

double appPi = trapezoidal(f, 0.0, 1.0);

double tol = 1e-6;
int max_refinements = 20;
double appPi2 = trapezoidal(f, 0.0, 1.0, tol, max_refinements);

cout << "Trapezoid Rule results for computing pi by integration" << endl;
cout << "a) with defaults, and b) with tol and max_refine set : " << endl;
cout << appPi << ", " << appPi2 << endl << endl;

Я привел два примера, один с использованием настроек по умолчанию для дискретизации диапазона интегрирования и сходимости, а второй с использованием пользовательских настроек.

Мои результаты (просто скопировав вывод с экрана):

Trapezoid Rule results for computing pi by integration
a) with defaults, and b) with tol and max_refine set :
3.14159, 3.14159

Я бы рекомендовал квадратурные и линейные функции формы Гаусса:

http://www.mathworks.com/matlabcentral/fileexchange/9230-gaussian-quadrature-for-triangles

http://www.wolframalpha.com/input/?i=gaussian+quadrature+triangle

http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf