InfiniteAdaptiveGaussKronrod Class

In applied mathematics, adaptive quadrature is a process in which the integral of a function f(x) is approximated using static quadrature rules on adaptively refined subintervals of the integration domain. Generally, adaptive algorithms are just as efficient and effective as traditional algorithms for "well behaved" integrands, but are also effective for "badly behaved" integrands for which traditional algorithms fail.

The algorithm implemented by this class has been based on the original FORTRAN implementation from QUADPACK. The function implemented the Non-adaptive Gauss- Kronrod integration is qagi(f,bound,inf,epsabs,epsrel,result,abserr,neval, ier,limit,lenw,last,iwork,work). The original source code is in the public domain, but this version is under the LGPL. The original authors, as long as the original routine description, are listed below:

Robert Piessens, Elise de Doncker; Applied Mathematics and Programming Division, K.U.Leuven, Leuvenappl. This routine calculates an approximation result to a given integral i = integral of f over (bound,+infinity) or i = integral of f over (-infinity,bound) or i = integral of f over (-infinity,+infinity) hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).

References:

• Wikipedia, The Free Encyclopedia. Adaptive quadrature. Available on: http://en.wikipedia.org/wiki/Adaptive_quadrature
• Wikipedia, The Free Encyclopedia. QUADPACK. Available on: http://en.wikipedia.org/wiki/QUADPACK
• Robert Piessens, Elise de Doncker; Non-adaptive integration standard fortran subroutine (qng.f). Applied Mathematics and Programming Division, K.U.Leuven, Leuvenappl. Available at: http://www.netlib.no/netlib/quadpack/qagi.f

Let's say we would like to compute the definite integral of the function f(x) = cos(x) in the interval -1 to +1 using a variety of integration methods, including the TrapezoidalRule, RombergMethod and NonAdaptiveGaussKronrod. Those methods can compute definite integrals where the integration interval is finite:

// Declare the function we want to integrate
Func<double, double> f = (x) => Math.Cos(x);

// We would like to know its integral from -1 to +1
double a = -1, b = +1;

// Integrate!
double trapez  = TrapezoidalRule.Integrate(f, a, b, steps: 1000); // 1.6829414
double romberg = RombergMethod.Integrate(f, a, b);                // 1.6829419
double nagk    = NonAdaptiveGaussKronrod.Integrate(f, a, b);      // 1.6829419

Moreover, it is also possible to calculate the value of improper integrals (it is, integrals with infinite bounds) using InfiniteAdaptiveGaussKronrod, as shown below. Let's say we would like to compute the area under the Gaussian curve from -infinite to +infinite. While this function has infinite bounds, this function is known to integrate to 1.

// Declare the Normal distribution's density function (which is the Gaussian's bell curve)
Func<double, double> g = (x) => (1 / Math.Sqrt(2 * Math.PI)) * Math.Exp(-(x * x) / 2);

// Integrate!
double iagk = InfiniteAdaptiveGaussKronrod.Integrate(g,
    Double.NegativeInfinity, Double.PositiveInfinity);   // Result should be 0.99999...

Reference