Double Exponential Rule Integration

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

The double exponential (DE) rule integrates by a change of variables chosen so that the plain trapezoidal rule, applied to the transformed integrand, converges exponentially fast — even when the original integrand has singularities at its endpoints, and without needing to tailor the substitution to the specific kind of singularity. It is described in Press et al., 2007, Section 4.5 (pages 172–178), and implemented as DoubleExponentialRuleQuadrature (and DoubleExponentialRuleMatrixQuadrature).

Quadrature by variable transformation

Consider and a change of variables mapping to :

If is chosen to go rapidly to zero at both endpoints of the range, then even a plain, uniformly-weighted trapezoidal rule on the transformed integral becomes extremely accurate — because the endpoint weights barely matter once the transformed integrand is already close to zero there. This works even for an with an integrable singularity at or : the vanishing factor can overwhelm it.

The first transformation of this kind, the TANH rule, maps to via:

Since decays like , repeating the same idea one level deeper — applying an exponential map to itself — decays even faster. This is the DE rule:

where the constant is typically or . The double-exponential decay of is what gives the rule its name and its markedly faster convergence than the single-exponential TANH rule — comparing their respective optimal errors, for TANH versus for DE, the DE rule wins for any reasonably large number of points .

Implementation details

Evaluating the formula above directly is numerically delicate: can overflow if computed as , and forming or for tiny can lose precision right where the integrand is most sensitive (near a singularity). This package follows the same approach as Press et al., 2007 to avoid both problems. Writing (for ):

the transformed integral is evaluated symmetrically about the interval’s midpoint, using directly instead of computing or and subtracting:

Because simply underflows harmlessly to zero for large positive , only needs to be evaluated directly; the symmetry of the trapezoidal rule about the interval’s midpoint supplies the side for free.

To let a singular integrand be coded robustly near or , this package’s evaluator interface passes explicitly alongside , so that (for example) can be coded as near rather than computing the cancellation-prone directly.

Algorithm

flowchart TD A["Stage 1: s = hmax * (b-a)/2 * f(midpoint, halfwidth)"] --> B["Stage n: double the number of t-samples\nwithin (-hmax, hmax)"] B --> C["For each new t: q = exp(-2 sinh t), delta = (b-a)*q/(1+q)"] C --> D["weight = q / (1+q)^2 * cosh(t)"] D --> E["sum += weight * (f(a+delta, delta) + f(b-delta, delta))"] E --> F["s = 0.5*s_prev + (b-a) * (spacing) * sum"] F --> G{"Converged?"} G -->|No| B G -->|Yes| H["Return s"]

Library implementation

DoubleExponentialRuleQuadrature transforms the range to the fixed interval , where hmax defaults to DEFAULT_HMAX = 3.7. It accepts either:

  • a DoubleExponentialSingleDimensionFunctionEvaluatorListener, whose evaluate(x, delta) receives the singularity-handling offset described above, or

  • a plain SingleDimensionFunctionEvaluatorListener, for integrands with no singularities or only "mild" ones (no worse than logarithmic), where can safely be ignored.

DoubleExponentialRuleMatrixQuadrature and DoubleExponentialMatrixSingleDimensionFunctionEvaluatorListener mirror both for matrix-valued integrands. As with every other quadrature rule in this package, the DE rule can be paired with any of the three integration methods — directly via DoubleExponentialRuleQuadratureIntegrator, or refined further by Simpson’s Rule Integration (SimpsonDoubleExponentialRuleQuadratureIntegrator) or Romberg Integration (RombergDoubleExponentialRuleQuadratureIntegrator) — or through the factory with QuadratureType.DOUBLE_EXPONENTIAL_RULE:

final var integrator = Integrator.create(
        a, b, listener, IntegratorType.ROMBERG, QuadratureType.DOUBLE_EXPONENTIAL_RULE);
final var result = integrator.integrate();
hmax trades range-of-integration coverage in the transformed variable against how aggressively the substitution suppresses the endpoints. The default of 3.7 is adequate for mild (e.g. logarithmic) singularities; stronger singularities, such as an inverse square-root, need a larger value (around 4.3) to reach full double-precision accuracy.

Example

Integrating , singular at both endpoints:

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

final var integrator = new DoubleExponentialRuleQuadratureIntegrator(0.0, 1.0, listener);
final var result = integrator.integrate(); // converges to machine precision in a few dozen evaluations