Algorithm

Problem class

DECUHR targets integrals of the form

where may have a homogeneous vertex singularity at the lower-left corner of the domain:

where    and is smooth near the corner. A logarithmic factor   is also supported.

Domain convention

The singularity must be located at the lower-left vertex of the integration domain, i.e. at . Shift the domain beforehand if necessary.

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 of a geometric series is appended. The tableau is extrapolated to accelerate convergence from 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:

key	Rule	Dimensions	Points per region
1	Degree-13	2D only	65
2	Degree-11	3D only	127
3	Degree-9	any nD
4	Degree-7	any nD
0	Auto	—	key=1 if n=2, key=2 if n=3, else key=3

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])

Alpha auto-estimation and AD

Automatic estimation of (triggered when alpha ≤ -singul) runs in Float64 on the primal integrand — is the structural homogeneity degree of the singularity and does not depend on the differentiation seed — after which the integration proceeds in the dual-number type. Differentiating through solve therefore works with or without an explicit alpha:

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:

Code	Meaning
0	Success — all tolerances satisfied
1	Budget (maxiters) exhausted
2–13	Invalid input parameters (key, ndim, numfun, bounds, budget, tolerances, singul, logf, minpts, emax, wrksub)
14	Integral may not converge (   even after auto-estimation)
15	emax < 1
17	wrksub too small for the given budget

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