Numerical integration in C++

Question

I need to integrate a function (of two variables). I know I can do it by using Fubini theorem to integrate one variable functions, then using numerical methods such as the Rectangle method or the Trapezoidal rule.

But are there any pre-built functions to do that in C++? I need to integrate over the unit R2 triangle ((0,0), (1,0), (0,1)).

Answer 1

You can use the GNU Scientific Library, which supports many "Numerical analysis" functions including integration.

A very simple example of integration from the manual is just a few lines of code:

#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;
}

Answer 2

As far as I know there are no functions of the type you are searching for in the standard library. Here is one implementation of the:

The extended trapezoidal rule.

For a fixed function f(x) to be integrated between fixed limits a and b, one can double the number of intervals in the extended trapezoidal rule without losing the benefit of previous work. The coarsest implementation of the trapezoidal rule is to average the function at its endpoints a and b. The first stage of refinement is to add to this average the value of the function at the halfway point. The second stage of refinement is to add the values at the 1/4 and 3/4 points.

A number of elementary quadrature algorithms involve adding successive stages of refinement. It is convenient to encapsulate this feature in a Quadrature structure:

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.
};

Then the Trapzd structure is derived from this as follows:

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

Usage:

The Trapzd structure could be used in several ways. The simplest and crudest is to integrate a function by the extended trapezoidal rule where you know in advance the number of steps you want. If you want 2^M + 1, you can accomplish this by the fragment:

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

with the answer returned as val. Here Ftor is a functor containing the function to be integrated.

Bonus:

Much better, of course, is to refine the trapezoidal rule until some specified degree of accuracy has been achieved. A function for this is:

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

Typedefs.

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

Answer 3

You might wish to check out the Boost Quadrature and Differentiation Library. In particular, they provide a version of the Trapezoid Rule:

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

The Quadrature/Differentiation Library is very well written and compatible with modern C++, in that one can just use a lambda expression or a function object for the integrand. I was up and running with it very quickly.

Here's an example of approximating pi, by approximating the integral of 4/(1 + x^2), from x = 0 to x = 1, putting the integrand as a lambda expression.

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

I provided two examples, one using the default settings for the discretization of the range of integration and convergence, and the second using custom settings.

My results (just taking a copy of the output from the screen) are:

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

Answer 4

I would recommend Gaussian quadrature and linear shape functions:

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