RombergMethod Class

In numerical analysis, Romberg's method (Romberg 1955) is used to estimate the definite integral ∫_a^b(x) dx by applying Richardson extrapolation repeatedly on the trapezium rule or the rectangle rule (midpoint rule). The estimates generate a triangular array. Romberg's method is a Newton–Cotes formula – it evaluates the integrand at equally spaced points. The integrand must have continuous derivatives, though fairly good results may be obtained if only a few derivatives exist. If it is possible to evaluate the integrand at unequally spaced points, then other methods such as Gaussian quadrature and Clenshaw–Curtis quadrature are generally more accurate.

References:

• Wikipedia, The Free Encyclopedia. Romberg's method. Available on: http://en.wikipedia.org/wiki/Romberg's_method

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