Simpson’s Rule Integration

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

Simpson’s rule is Richardson extrapolation applied to a quadrature rule with an error series that contains only even powers of the step size — obtained here for the cost of essentially one extra combination step per refinement, with no additional function evaluations. It is described in Press et al., 2007, Sections 4.1 and 4.2 (pages 158, 160, 164–165), and implemented as the generic SimpsonIntegrator base class (and SimpsonMatrixIntegrator), specialized to each quadrature rule in the package.

Problem statement

The classical three-point Simpson’s rule integrates over a pair of intervals:

and, applied repeatedly, gives the extended Simpson’s rule with its characteristic alternating , interior weights:

The "deep fact" that makes this cheap

Trapezoidal Rule Integration noted that the extended trapezoidal rule’s error series contains only even powers of . Suppose the trapezoidal rule is evaluated with steps, giving , and again with steps, giving . Since doubling the steps quarters the leading ( ) error term, the combination:

cancels that leading term exactly — and because there is no odd-order term to leave behind, the next surviving error term is . This combination is exactly Simpson’s rule; it is also the numerically preferred way to evaluate it, since it reuses two successive trapezoidal-style refinements instead of recomputing the alternating weights from scratch.

This argument only depends on the underlying rule’s error series being even in the step size — a property shared by every quadrature rule in this package (each one is derived from either the extended trapezoidal rule or the extended midpoint rule, both of which have this property). That is why SimpsonIntegrator<T extends Quadrature> in the library is generic over any Quadrature, not hardwired to the trapezoidal rule.

Algorithm

flowchart TD A["st_n = quadrature.next() (stage n estimate from the underlying rule)"] --> B["s_n = (4 * st_n - st_(n-1)) / 3"] B --> C{"|s_n - s_(n-1)| < eps * |s_(n-1)| ?\n(checked after 5 stages, up to 20)"} C -->|No| A C -->|Yes| D["Return s_n"]

Library implementation

SimpsonIntegrator<T extends Quadrature> wraps any of the package’s quadrature rules and applies the combination every stage, checking convergence after a minimum of 5 stages and up to a maximum of 20, to a relative tolerance eps (default 1e-10). SimpsonMatrixIntegrator does the same for matrix-valued integrands. Six concrete classes pair this method with each applicable quadrature rule:

Class Quadrature rule

SimpsonTrapezoidalQuadratureIntegrator

Trapezoidal Rule Integration — finite, non-singular intervals.

SimpsonMidPointQuadratureIntegrator

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

SimpsonInfinityMidPointQuadratureIntegrator

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

SimpsonLowerSquareRootMidPointQuadratureIntegrator

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

SimpsonUpperSquareRootMidPointQuadratureIntegrator

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

SimpsonDoubleExponentialRuleQuadratureIntegrator

Double Exponential Rule Integration — arbitrary endpoint singularities.

Each has a Matrix-valued counterpart. Any of them can also be obtained through the package’s factory methods:

final var integrator = Integrator.create(
        a, b, listener, IntegratorType.SIMPSON, QuadratureType.TRAPEZOIDAL);
final var result = integrator.integrate();
IntegratorType.SIMPSON is generally more efficient than plain IntegratorType.QUADRATURE (fewer function evaluations for the same accuracy) whenever the integrand has a finite fourth derivative — i.e. a continuous third derivative. For light-duty work on reasonably smooth functions, Simpson’s method paired with the trapezoidal rule is a good default; for very smooth (e.g. analytic) integrands, Romberg Integration typically does even better.

Example

final var listener = new SingleDimensionFunctionEvaluatorListener() {
    @Override
    public double evaluate(final double point) {
        return Math.sin(point);
    }
};

final var integrator = new SimpsonTrapezoidalQuadratureIntegrator(0.0, Math.PI, listener);
final var result = integrator.integrate(); // approximately 2.0, in fewer stages than the plain quadrature integrator