Algorithm

Problem class

DECUHR targets integrals of the form

\[I = \int_{\mathbf{a}}^{\mathbf{b}} f(\mathbf{x})\, d\mathbf{x}\]

where $f$ may have a homogeneous vertex singularity at the lower-left corner $(a_1, \ldots, a_s)$ of the domain:

\[f(\mathbf{x}) \;\sim\; g(\mathbf{x})\, \prod_{i=1}^{s}(x_i - a_i)^{\alpha} \quad \text{as } \mathbf{x} \to \mathbf{a}_{1:s},\]

where $\alpha > -s$ and $g$ is smooth near the corner. A logarithmic factor $\log\!\bigl(\prod_i(x_i-a_i)\bigr)$ is also supported.

Strategy

The algorithm combines three ideas:

1. Adaptive subdivision (DEADHR / DESBHR). The domain is maintained as a pool of subregions stored in a binary max-heap keyed on the local error estimate. At each step the worst-error region is subdivided:

   • Singular region (touching the singularity corner): cut into singul + 1 children by halving each singular dimension in turn.
   • Regular region: bisect along the axis with the largest fourth difference.

2. Richardson extrapolation (DEXTHR). Each time the singular region is subdivided, a new term $U_k$ of a geometric series is appended. The tableau is extrapolated to accelerate convergence from $O(h^\alpha)$ to machine precision.

3. Automatic singularity estimation (DECALP). When alpha ≤ -singul (the default), the exponent is estimated by evaluating the integrand along a ray toward the singular vertex and extrapolating the log-ratio sequence with Aitken / Bjorstad–Grosse–Dahlquist acceleration.

Integration rules

Four Genz–Malik-style fully-symmetric rules are available, selected via key:

Result, error estimate and diagnostics

DECUHR's error estimator is deliberately conservative: following the original Fortran DEXTHR, the pure extrapolation error is inflated by a heuristic constant before being compared with the tolerance. A returned retcode = MaxIters therefore frequently accompanies a result that is in fact accurate to better than the requested tolerance — the estimate simply never dropped below it within the evaluation budget.

Practical guidance:

• Trust sol.u; read sol.resid for the (conservative) error estimate.
• Raise maxiters and/or wrksub to let the estimate catch up to the true accuracy if a certified Success is required.
• sol.stats.numevals reports the number of integrand evaluations and sol.stats.ifail the raw DECUHR code — handy to confirm you hit the budget rather than a genuine problem.

Generic value type and automatic differentiation

All internal arrays that accumulate integrand values are parameterized on a value type TV, inferred at runtime from a test evaluation:

TV = typeof(f(midpoint, p))   # e.g. Float64 or ForwardDiff.Dual{...}

This means the algorithm works out of the box with dual numbers when differentiating through the integral with respect to parameters p:

using ForwardDiff, Integrals, DECUHR

f = (u, p) -> exp(-p[1] * (u[1]^2 + u[2]^2))
prob = IntegralProblem(f, (zeros(2), ones(2)))

# Derivative of the integral w.r.t. p[1] at p₀ = [1.0]
dI_dp = ForwardDiff.gradient(p -> solve(prob, DecuhrAlgorithm(alpha=0.0);
                                         abstol=1e-8).u, [1.0])

DecuhrAlgorithm(singul=2, alpha=-0.5)   # explicit α
DecuhrAlgorithm(singul=2)               # α auto-estimated — also differentiable

Error codes (ifail)

The low-level driver returns an integer ifail:

These are mapped to SciMLBase.ReturnCode values by the Integrals.jl integration hook: ifail=0 → Success, ifail=1 → MaxIters, otherwise → Failure.