C++ library for numerical integration

libIntegrate is a collection of methods for numerical integration, including discretized data.

Features:

Integrate one and two dimensional discretized functions, with Riemann, Trapezoid, and Simpson rules.
Simple, clean, uniform interface. Integration methods are implemented as classes with operator()(...) methods.
Handles several common random access, continuous memory container interfacess (std::vector, Eigen::Vector, etc). Containers that provide element access with operator[](int) or operator()(int) and a method for getting the size (.size(), .length(), .rows(), etc) can be integrated.
Support for non-uniform discretization.
Apply lazy transformations to discrete data before integrating. Useful for computing weighted integrals.

Description
Usage

Description

It is often necessary to integrate functions that have been discretized when doing physics simulations or working with experimental data. For example, we might use a finite-difference method to solve Schrodinger's equation and compute the time-evolution of a wavefunction describing some system. If we then want to compute probabilities, we will need to integrate the wavefunction (actually, we will need to integrate the modulus squared of the wavefunction). Or, we might have a measured spectrum for a broadband source that we need to integrate to compute the total power emitted. In either case, we will have a discrete set of points to use for the integration.

Most libraries for numerical integration provide routines for integrating functions that can be evaluated, not for functions that have already been discretized, and are therefore not suitable for this. This library provides methods that are.

If the function you need to integrate is already discretized, there are not many options to choose from. There is the Riemann sum, the Trapezoid rule, and perhaps Simpson's rule. Each of these methods basically interpolate the function over some interval and use the integral of the interpolating function to estimate the integral of the data. Riemann sum uses a constant (zero'th order polynomial), the Trapezoid rule uses a linear function (first order polynomial), and the 1/3 Simpson's rule uses a quadratic (second degree polynomial). Higher-order interpolating functions are of course possible, but may not be appropriate, just as fitting high order polynomials to data sets is usually not appropriate.

Quickly writing a Riemann sum loop to compute the integral of some discretized function is pretty easy. Writing the Trapezoid rule isn't difficult either, but you need to be careful not to index out of range at the end of the loop. Simpson's rule is simple enough for uniform data, but again requires care at the end of the loop. For non-uniform data, it is a bit more involved. I noticed that, when needing to integrate some discretized function, I would often just use a Riemann sum and say to myself "I'll come back and use a better method after I get this working", but then wouldn't. So this library was created to make using "better" methods on discretized data easy.

While this library does provide methods for integrating functions that can be evaluated too, they are not nearly as as sophisticated (actually, not very sophisticated at all), or complete (actually, not very complete at all), as the Gnu Scientific Library (https://www.gnu.org/software/gsl/), Numerical Recipes (http://numerical.recipes/), QUADPACK (http://www.netlib.org/quadpack/), or Boost (https://www.boost.org/doc/libs/1_75_0/libs/math/doc/html/quadrature.html). If you need to integrate a known function with some requirement on accuracy, you should use one of these libraries instead.

If you need to integrate discretized data and occasionally some simple functions that can be evaluated, this library may be for you. It is released under a permissive license, so you are free to do whatever you want with it.

Methods currently implemented include:

1D
- Discretized Functions
  - Riemann sum (the one that every undergrad physics major writes)
  - Trapezoid rule
  - Simpson's rule (1/3)
- Callable Functions
  - Riemann sum
  - Trapezoid rule
  - Simpson's rule (1/3)
  - Gauss-Legendre Quadrature of order 8, 16, 32, and 64.
2D
- Discretized Functions
  - Riemann sum
  - Trapezoid rule
  - Simpson's rule (1/3)
- Callable Functions
  - Riemann sum
  - Trapezoid rule
  - Simpson's rule (1/3)
  - Gauss-Legendre Quadrature of order 8, 16, 32, and 64.

Note that the library depends on Boost, and does provide some (incomplete) wrappers to the Boost.Math quadrature functions.

1D
- Gauss-Kronrod Quadrature (arbitrary order)

Usage

Installing

libIntegrate is a header-only library that uses CMake for building unit tests and installing. To install the library, clone this and checkout the latest release, or download one of the release tarballs. Then,

$ cd libIntegrate
$ mkdir build
$ cd build
$ cmake ..
$ cmake --build . --target install

Conan packages are also provided. To install the library with conan, add the remote https://cdc3.jfrog.io/artifactory/api/conan/default-conan and include libintegrate/VERSION@cd3/devel to you conanfile.txt or conanfile.py.

$ conan remote add cd3 https://cdc3.jfrog.io/artifactory/api/conan/default-conan
$ conan search -r cd3 # list all packages

Overview

There are two basic types of numerical integration. In the first case, the value of the function to be integrated can be calculated for any argument value. In the second case, the value of the function to be integrated is only know at some discrete set of argument values. Most numerical integration libraries handle the first case, and there many different methods for handling a wide variety of conditions, including integrating over infinite or semi-infinite bounds. Typically, the user should know something about the function to be integrated and carefully select the best method to use.

Any method that can handle the second case can be used to handle the first case as well by first discrediting the function to be integrated. Methods that can only handle the first case (for example, Gauss-Legendre Quadrature requires the function to be evaluated at specific argument values) can be used on discretized functions, but require the discretized function to be interpolated to the argument values required by the method first. The libInterpolate library is useful for this.

Examples

Discretized Functions

In the simplest case, you have a discretized function with the argument and function values stored in two vectors.

std::vector<double> x(100), f(100);
double dx = M_PI/(x.size()-1);
for(int i = 0; i < x.size(); i++)
{
  x[i] = i*dx;
  f[i] = sin(x[i]);
}

To integrate the function, create an integrator object and call it on the discretized function.

#include <libIntegrate/Integrate.hpp>

_1D::RiemannRule<double> integrage;

double I = integrate(x,f); // I \approx 2

Note: Discretized functions must be sorted with respect to the independent variable (x in this example).

Each method is implemented in a class. To use a different method, instantiate the corresponding class.

_1D::TrapezoidRule<double> trap_integrage;
_1D::SimpsonRule<double> simp_integrage;

Now, in practice, you would not need to manually discretize a function that can be evaluated. All of the methods support integrating a callable with a specified number of intervals. The example above is equivalent to

#include <libIntegrate/Integrate.hpp>

_1D::RiemannRule<double> integrage;

double I = integrate([](double x){return sin(x);}, 0, M_PI, 99);

However, when the discretized function is computed by solving a differential equation, or is data that was measured in an experiment, the discretized integration is required.

Two-dimensional discretized functions can also be integrated. The 2D integrators that support discretized functions take the arguments. The first two are 1D containers of the x and y coordinates, the third is a 2D container of the function values. For example:

  _2D::TrapezoidRule<double> integrate;
  double                     I;

  // create a discretized version of sin(x)*sin(y) to integrate.
  // probably not the best way to create an initialize the function...
  std::vector<double>              x(100), y(200);
  std::vector<std::vector<double>> f(100);
  for(size_t i = 0; i < x.size(); i++) {
      x[i] = i * (M_PI / 2) / (x.size() - 1);
  }
  for(size_t i = 0; i < y.size(); i++) {
      y[i] = i * (M_PI / 2) / (y.size() - 1);
  }

  for(size_t i = 0; i < f.size(); i++) {
      f[i] = std::vector<double>(200);
      for(size_t j = 0; j < f[i].size(); j++) {
        f[i][j] = sin(x[i]) * sin(y[j]);
      }
  }

  // now integrate
  I = integrate(x, y, f);

We could also use more efficient 2D container here, like an E

LibIntegrate

Install / Use

README