Romberg Integration

This documentation was generated with the assistance of AI. Please report any inaccuracies.

Romberg integration generalizes Simpson’s rule from a single Richardson-extrapolation step to a full polynomial extrapolation of several successive refinements to the continuum limit of zero step size. It is described in Press et al., 2007, Section 4.3 (page 166), generalized to open-interval quadrature rules in Section 4.4 (pages 168–169), and implemented as the generic RombergIntegrator base class (and RombergMatrixIntegrator), specialized to each quadrature rule in the package.

From Simpson’s rule to Romberg’s method

Simpson’s Rule Integration combines two successive trapezoidal estimates to cancel the leading error term, using a single fixed formula. Romberg’s method is the natural generalization: use the results of successive refinements of the underlying quadrature rule and remove all the error terms up to (but not including) the one of order . Simpson’s rule is exactly the case.

This is an instance of a general technique called Richardson’s deferred approach to the limit: run a numerical algorithm for several values of a step-size parameter , then extrapolate the sequence of results to the limit . The combination used by Simpson’s rule is a special case of polynomial extrapolation with only two points; Romberg’s method instead fits a polynomial of degree through the last points — using Neville’s algorithm — and evaluates that polynomial at to estimate the value in the continuum limit, along with an error estimate for the extrapolation itself.

Algorithm

flowchart TD A["h[0] = 1"] --> B["s[j] = quadrature.next() (stage j+1 estimate)"] B --> C{"j >= K (5 points collected) ?"} C -->|No| D["h[j] = h[j-1] * decay factor"] --> B C -->|Yes| E["Fit a degree-(K-1) polynomial through\nthe last K points (h_i, s_i) via Neville's algorithm"] E --> F["Extrapolate the polynomial to h = 0"] F --> G{"extrapolation error <= eps * |estimate| ?"} G -->|No| D G -->|Yes| H["Return extrapolated estimate"]

The step-size sequence is relative, not absolute — what matters is its ratio between stages, which must match how the underlying quadrature rule actually refines itself from one stage to the next, so that the polynomial extrapolation is fit against the correct powers of the true step size:

  • Quadrature rules that double their number of evaluations each stage (like Trapezoidal Rule Integration) use a decay factor of per stage — this is what RombergTrapezoidalQuadratureIntegrator overrides the base class to do.

  • Every other quadrature rule in this package uses `RombergIntegrator’s default decay factor of per stage, matching the step-tripling refinement of Midpoint Rule Integration for Improper Integrals and its variants.

Library implementation

RombergIntegrator<T extends Quadrature> drives any of the package’s quadrature rules through this extrapolation loop, reusing the library’s own PolynomialInterpolator (from the interpolation package) to perform the Neville-algorithm extrapolation, rather than reimplementing it inline. Convergence is checked from the 5th stage onward, up to a maximum of 14 stages (RombergTrapezoidalQuadratureIntegrator uses its own, larger maximum of 20 stages, matching its different decay factor), with a default relative tolerance eps of 3.0e-9 (1e-10 for the trapezoidal-based variant). RombergMatrixIntegrator provides the same algorithm for matrix-valued integrands, running one PolynomialInterpolator per matrix element and requiring all of them to converge together.

Class Quadrature rule

RombergTrapezoidalQuadratureIntegrator

Trapezoidal Rule Integration — finite, non-singular intervals. The library’s default integrator + quadrature combination.

RombergMidPointQuadratureIntegrator

Midpoint Rule Integration for Improper Integrals — open formula for improper integrals.

RombergInfinityMidPointQuadratureIntegrator

Midpoint Rule Integration for Improper Integrals — semi-infinite or infinite ranges.

RombergLowerSquareRootMidPointQuadratureIntegrator

Midpoint Rule Integration for Improper Integrals — inverse square-root singularity at the lower bound.

RombergUpperSquareRootMidPointQuadratureIntegrator

Midpoint Rule Integration for Improper Integrals — inverse square-root singularity at the upper bound.

RombergExponentialMidPointQuadratureIntegrator

Midpoint Rule Integration for Improper Integrals — exponentially decaying integrand up to infinity. Its constructor takes only the lower limit; the upper limit is always treated as infinite.

RombergDoubleExponentialRuleQuadratureIntegrator

Double Exponential Rule Integration — arbitrary endpoint singularities.

Each has a Matrix-valued counterpart. IntegratorType.ROMBERG paired with QuadratureType.TRAPEZOIDAL is this package’s default integrator/quadrature combination (Integrator.DEFAULT_INTEGRATOR_TYPE / Integrator.DEFAULT_QUADRATURE_TYPE), used whenever a caller does not specify one explicitly:

final var integrator = Integrator.create(a, b, listener); // Romberg + trapezoidal, default eps
final var result = integrator.integrate();
Romberg integration is the most efficient method in this package for sufficiently smooth (e.g. analytic) integrands over intervals with no singularities, including at the endpoints — in such cases it can take many, many fewer function evaluations than Simpson’s Rule Integration or the plain quadrature integrator. It is, correspondingly, the least forgiving option: it tends to converge poorly, or not at all, when the integrand is not smooth or contains a singularity that the chosen quadrature rule was not designed to handle. Pick a quadrature rule from Midpoint Rule Integration for Improper Integrals or Double Exponential Rule Integration that matches the integrand’s actual singularities before reaching for Romberg’s method.

Example

final var listener = new SingleDimensionFunctionEvaluatorListener() {
    @Override
    public double evaluate(final double point) {
        final var x2 = point * point;
        return Math.pow(point, 4.0) * Math.log(point + Math.sqrt(x2 + 1.0));
    }
};

final var integrator = new RombergTrapezoidalQuadratureIntegrator(0.0, 2.0, listener);
final var result = integrator.integrate();