Numerical integration is an important tool in many applications, this article is an introduction to implementing several integration algorithms for univariable functions. The article finishes with a description of how to implement the Gauss Quadrature in D, which is the standard method used to carry out numerical integration. C’s LAPACKE library is used to carry out eigenvalue decomposition in order to evaluate the coefficients and points of the quadrature. In the process the reader will observe that calling C libraries from D is simple. As usual, the article focuses on the implementation details, we do not compare the accuracy of the various methods, only show their outputs because it is apparent that as the standard method for integration, Gauss Quadrature is the most accurate.
Enjoy!
References for the algorithms used in this article can be found below:
These analytical functions and limits will be calculated using the various integration methods:
$$ \int_{-2}^{2}{x^{3}e^{x}}\mathrm{d}x = \bigg[e^{x}(x^{3} - 3x^{2} + 6x - 6)\bigg]_{-2}^{2} $$
$$ = 19.920852908526\ldots $$
$$ \int_{0}^{2}{e^{2x}}\sin{(3x)}\mathrm{d}x = \bigg[\frac{1}{13}e^{2x}(2\sin{(3x)} - 3\cos{(3x)}\bigg]_{0}^{2} $$
$$ = -14.2139771298625\ldots $$
First we import some functions and a implement simple template predicate:
import std.stdio: writeln;
import std.conv: to;
import std.format: format;
import core.stdc.math: exp = expf, exp, exp = expl,
cos = cosf, cos, cos = cosl,
pow = powf, pow, pow = powl,
sin = sinf, sin, sin = sinl;
template isFloatingPoint(T)
{
enum isFloatingPoint = is(T == float) || is(T == double) || is(T == real);
}
The functions to be integrated are given as:
auto f1 = (double x) => exp(x)*(x^^3);
auto f2 = (double x) => exp(2*x)*sin(3*x);
If we subdivide the range of \(x \in [a, b]\) into \(N\) segments, the rectangular method for integration is given by:
$$ I(f) = \int_{a}^{b} \! f(x) \, \mathrm{d}x \approx h \sum_{i = 1}^{N}{f(x_{i})} $$
The D code for implementing rectangular integration is given below:
T rectangular(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
T h = (b - a)/N;
T sum = fun(a, u); T x = a;
for(long i = 1; i < (N - 1); ++i)
{
x += h;
sum += fun(x, u);
}
x += h;
sum += fun(x, u);
return sum*h;
}
call:
writeln("Rectangular integral I approximation: ",
format("%.13f", rectangular!(f1)(-2.0, 2.0, 1000)));
writeln("Rectangular integral II approximation: ",
format("%.13f", rectangular!(f2)(0.0, 2.0, 1000)));
outputs:
Rectangular integral I approximation: 19.8006590182198
Rectangular integral II approximation: -14.1986803074543
Even by eye we can tell that 1000 subdivisions of the range gives a rather poor result. Note that this is 1000 evaluations of the target function, which is impractical for very expensive functions.
The midpoint method for integration is given by:
$$ I(f) = \int_{a}^{b} \! f(x) \, \mathrm{d}x \approx h \sum_{i = 1}^{N}{f\bigg(\frac{x_{i} - x_{i + 1}}{2}\bigg)} $$
The D code for implementing midpoint integration is given below:
T midpoint(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
T h = (b - a)/N;
T sum = 0; T x = a;
for(long i = 0; i < N; ++i)
{
sum += fun(x + (h/2), u);
x += h;
}
return sum*h;
}
the call:
writeln("Midpoint integral I approximation: ",
format("%.13f", midpoint!(f1)(-2.0, 2.0, 1000)));
writeln("Midpoint integral II approximation: ",
format("%.13f", midpoint!(f2)(0.0, 2.0, 1000)));
the output:
Midpoint integral I approximation: 19.9207548011967
Midpoint integral II approximation: -14.2139977564108
The trapezoidal method for integration is given by:
$$ I(f) = \int_{a}^{b} \! f(x) \, \mathrm{d}x \approx \frac{h}{2}\bigg[f(a) + f(b)\bigg] + h\sum_{i = 2}^{N}{f(x_{i})} $$
The D code for implementing trapezoidal integration is given below:
T trapezoidal(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
T h = (b - a)/N; T x = a;
T sum = (fun(a, u) + fun(b, u))/2;
for(long i = 1; i < N; ++i)
{
sum += fun(x += h, u);
}
return sum*h;
}
the call:
writeln("Trapezoidal integral I approximation: ",
format("%.13f", trapezoidal!(f1)(-2.0, 2.0, 1000)));
writeln("Trapezoidal integral II approximation: ",
format("%.13f", trapezoidal!(f2)(0.0, 2.0, 1000)));
output:
Trapezoidal integral I approximation: 19.9210492803345
Trapezoidal integral II approximation: -14.2139358767466
Simpson’s \(\frac{1}{3}\) method is given by:
$$ I(f) = \int_{a}^{b} \! f(x) \, \mathrm{d}x \approx \frac{h}{3}\bigg[f(a) + 2\sum_{j = 1}^{\frac{N}{2} - 1}{f(x_{2j})} + 4\sum_{j = 1}^{\frac{N}{2}}{f(x_{2j - 1})} + f(b)\bigg] $$
Note that \(N\) must be a multiple of \(2\).
The D code for implementing Simpson’s third integration is given below:
T simpsonsThird(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
debug assert(N%2 == 0, "N must be even!");
long m = (N/2) - 1;
T h = (b - a)/N; T x = a + h;
T sum = fun(a, u) + fun(b, u) + 4*fun(x, u);
for(long i = 0; i < m; ++i)
{
// Even
sum += 2*fun(x += h, u);
// Odd
sum += 4*fun(x += h, u);
}
return sum*h/3;
}
the call:
writeln("Simpson's Third integral I approximation: ",
format("%.13f", simpsonsThird!(f1)(-2.0, 2.0, 1000)));
writeln("Simpson's Third integral II approximation: ",
format("%.13f", simpsonsThird!(f2)(0.0, 2.0, 1000)));
output:
Simpson's Third integral I approximation: 19.9208529617569
Simpson's Third integral II approximation: -14.2139771297590
Simpson’s \(\frac{3}{8}\) method is given by:
$$ I(f) = \int_{a}^{b} \! f(x) \, \mathrm{d}x \approx \frac{3h}{8}\bigg[f(a) + 3\sum_{i \neq 3k}^{N - 1}{f(x_{i})} + 2\sum_{j = 1}^{\frac{N}{3} - 1}{f(x_{3j})} + f(b)\bigg] $$
Note that \(N\) must be a multiple of \(3\).
The D code for implementing Simpson’s three-eights integration is given below:
T simpsonsThreeEights(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
debug assert(N%3 == 0, "N must be divisible by three!");
long m = (N/3) - 1;
T h = (b - a)/N; T x = a + h;
T sum = fun(a, u) + fun(b, u) + 3*(fun(x, u) + fun(x += h, u));
for(long i = 0; i < m; ++i)
{
sum += 2*fun(x += h, u);
x += h;
sum += 3*(fun(x, u) + fun(x += h, u));
}
return (3*sum*h)/8;
}
the call:
writeln("Simpson's Three Eights integral I approximation: ",
format("%.13f", simpsonsThreeEights!(f1)(-2.0, 2.0, 999)));
writeln("Simpson's Three Eights integral II approximation: ",
format("%.13f", simpsonsThreeEights!(f2)(0.0, 2.0, 999)));
output:
Simpson's Three Eights integral I approximation: 19.9208529628933
Simpson's Three Eights integral II approximation: -14.2139771296285
Gauss Quadrature calculates integrals of the form:
$$ I(f) = \int_{-1}^{1}{f(x)} \,\mathrm{d}x $$
and does this using the approximation:
$$ I_{n}(f) = \sum_{i = 1}^{n}{w_{i}f(x_{i})} $$
Consider the a Tridiagonal matrix \(\{T_{n}\}\) of size \(n\), with diagonal \(\alpha_{i}\) and immediate off diagonal \(\beta_{i}\):
$$ \alpha_{i} = 0, \quad \beta_{i} = \frac{1}{2}(1 - (2i)^{-2})^{-1/2} $$
all other matrix entries are zero. Let \(T_{n} = VDV^{T}\) be an orthogonal diagonalization of \(T_{n}\) with \(D = diag(\lambda_{1}, \ldots, \lambda_{n})\), and \(V = |v_{1}|\ldots|v_{n}|\). The points \(x_{i} = \lambda_{i}\) and weights \(w_{i} = 2(v_{i})_{1}^{2}, \quad i = 1, \ldots, n\). Meaning that the points are given by the eigen values and the weights are given by the square of the first element of the eigen vector multiplied by two.
We can transform an integral in the range \([a, b]\) to the range \([-1, 1]\) by using the transformation:
$$ \int_{a}^{b}{f(x)}\mathrm{d}x = \int_{-1}^{1}{f\bigg(\frac{(b - a)t + b + a}{2}\bigg)}\frac{(b - a)}{2}\,\mathrm{d}t $$
$$ t = \frac{(2x - a - b)}{(b - a)} $$
$$ x = \frac{1}{2}[t(b - a) + a + b] $$
$$ \mathrm{d}x = \frac{1}{2}(b - a)\,\mathrm{d}t $$
Calling C from D is simple, we use extern (C) in a similar way to C++:
extern(C) @nogc nothrow{
int LAPACKE_sstev(int matrix_layout, char jobz, int n, float* d, float* e, float* z, int ldz);
int LAPACKE_dstev(int matrix_layout, char jobz, int n, double* d, double* e, double* z, int ldz);
}
The above function calls the relevant routine from LAPACK to calculate the eigenvalue decomposition for a tri-diagonal matrix. This link gives the documentation of the function.
To compile, you’ll need to add the relevant compilation flags, for example
dmd script.d -L-lopenblas -L-lpthread -L-llapacke -L-llapack -L-lm
The code below shows the rapper function to obtain the eigenvalues and eigen vectors from the C functions.
auto eigen(T)(T[] diag, T[] sub)
if(isFloatingPoint!(T))
{
int n = cast(int)diag.length;
T[] vec = new T[n*n];
int info = -1;
static if(is(T == float))
{
//102 for column major
info = LAPACKE_sstev(102, 'V', n, diag.ptr, sub.ptr, vec.ptr, n);
}else static if(is(T == double)){
info = LAPACKE_dstev(102, 'V', n, diag.ptr, sub.ptr, vec.ptr, n);
}else{
static assert(0, "No current implementation of eigen for type " ~ T.stringof ~ ".");
}
if(info != 0)
assert(0, "Something went wrong, info = " ~ to!(string)(info) ~ ".");
return [diag, vec];
}
This is the code for calculating the sub-diagonal of the tridiagonal matrix.
T[] subDiag(T)(size_t n)
if(isFloatingPoint!(T))
{
T[] sub = new T[n];
for(long i = 1; i <= n; ++i)
sub[i - 1] = 0.5*pow(1 - pow(2.0*i, -2.0), -0.5);
return sub;
}
Next the function for calculating the quadrature coefficients:
auto calcQuadratureCoefficients(T)(long n)
if(isFloatingPoint!(T))
{
T[] alpha = new T[n];
alpha[] = 0;
T[] beta = subDiag!(T)(n - 1);
auto ev = eigen(alpha, beta);
T[] eValues = new T[n];
for(long i = 0; i < n; ++i)
{
eValues[i] = ev[1][i*n];
}
return [ev[0], eValues];
}
Quick function for transforming the coordinate:
T[] transformX(T)(ref T[] x, T a, T b)
if(isFloatingPoint!(T))
{
T ba = b - a;
for(long i = 0; i < x.length; ++i)
x[i] = 0.5*((x[i]*ba) + a + b);
return x;
}
The D code for implementing Gauss Quadrature integration is given below:
T gaussQuadrature(alias fun, T, U...)(T a, T b, long N, U u)
if(isFloatingPoint!(T))
{
auto qCoeff = calcQuadratureCoefficients!(T)(N);
T[] x = transformX(qCoeff[0], a, b);
T[] weights = qCoeff[1];
T sum = 0;
for(long i = 0; i < N; ++i)
sum += 2*pow(weights[i], 2.0)*fun(x[i], u);
T ba = (b - a)/2;
return sum*ba;
}
the call:
writeln("Gauss Quadrature integral I approximation: ",
format("%.13f", gaussQuadrature!(f1)(-2.0, 2.0, 15)));
writeln("Gauss Quadrature integral II approximation: ",
format("%.13f", gaussQuadrature!(f2)(0.0, 2.0, 15)));
output:
Gauss Quadrature integral I approximation: 19.9208529608526
Gauss Quadrature integral II approximation: -14.2139771298625
Note that Gauss Quadrature only uses 15 evaluations of the target function in comparison with 1000 for the other methods and returns the most accurate result.
This article shows how simple it is to write numerical integration algorithms in D from scratch, and also how easy it is to call LAPACK/BLAS algorithms, basically just using a declaration as would be done in a C++ library, only there is no need to faff around with header files, and the relevant flags are easy to pass to the D compiler.
Thank you for reading!