Examples

The DECUHR algorithm as a pluggable Integrals.jl solver, on integrands whose exact value is known — vertex singularities in two and three dimensions, a logarithmic singularity, vector-valued and parametrized integrands, forward-mode differentiation through the quadrature, and what the return codes and diagnostics actually mean.

This file is the single source of both the documentation page you may be reading and a standalone script: run it with

julia --project=. examples/basic_usage.jl

or, from a REPL, include("examples/basic_usage.jl").

using Integrals
using DECUHR
using Printf

1 — Vertex singularity 2D, known

Analytical integral:

f = (u, p) -> (u[1] * u[2])^(-0.5)
prob = IntegralProblem(f, (zeros(2), ones(2)))
sol = solve(prob, DecuhrAlgorithm(singul = 2, alpha = -0.5); abstol = 1.0e-8)

println("I  ≈ ", sol.u)
println("|err| = ", abs(sol.u - 4.0))
println("retcode : ", sol.retcode)

I  ≈ 4.000121585481047
|err| = 0.00012158548104679312
retcode : MaxIters

2 — Automatic estimation of

Same integral, but without providing : DECUHR estimates it via DECALP.

sol2 = solve(prob, DecuhrAlgorithm(singul = 2); abstol = 1.0e-7)

println("I  ≈ ", sol2.u)
println("|err| = ", abs(sol2.u - 4.0))

I  ≈ 4.000000011204702
|err| = 1.1204702055067628e-8

3 — Smooth integrand (no singularity)

f3 = (u, p) -> sin(u[1]) * cos(u[2])
prob3 = IntegralProblem(f3, (zeros(2), fill(π / 2, 2)))
sol3 = solve(prob3, DecuhrAlgorithm(); abstol = 1.0e-10)

println("I  ≈ ", sol3.u)
println("|err| = ", abs(sol3.u - 1.0))

I  ≈ 1.0000000000000038
|err| = 3.774758283725532e-15

4 — Vector-valued integrand (NUMFUN = 2)

Two components integrated simultaneously:

f4 = (u, p) -> [u[1]^2 + u[2]^2, u[1] * u[2]]
prob4 = IntegralProblem(f4, (zeros(2), ones(2)))
sol4 = solve(prob4, DecuhrAlgorithm(); abstol = 1.0e-9)

exact4 = [2 / 3, 1 / 4]
println("I  ≈ ", sol4.u)
println("|err| = ", abs.(sol4.u .- exact4))

I  ≈ [0.6666666666666667, 0.2500000000000001]
|err| = [1.1102230246251565e-16, 1.1102230246251565e-16]

5 — 3D singularity

The 3-D vertex singularity is challenging; wrksub=50000 (now the default) allows sufficient refinement:

f5 = (u, p) -> (u[1] * u[2] * u[3])^(-1 / 3)
prob5 = IntegralProblem(f5, (zeros(3), ones(3)))
sol5 = solve(
    prob5, DecuhrAlgorithm(singul = 3, alpha = -1 / 3);
    abstol = 1.0e-7, reltol = 1.0e-7, maxiters = 1_500_000
)

println("I  ≈ ", sol5.u)
println("|err| = ", abs(sol5.u - (3 / 2)^3))

I  ≈ 3.3750066940189907
|err| = 6.69401899067168e-6

6 — Logarithmic singularity

f6 = (u, p) -> -log(u[1] * u[2])
prob6 = IntegralProblem(f6, (zeros(2), ones(2)))
sol6 = solve(prob6, DecuhrAlgorithm(singul = 2, alpha = 0.0, logf = 1); abstol = 1.0e-8)

println("I  ≈ ", sol6.u)
println("|err| = ", abs(sol6.u - 2.0))

I  ≈ 2.0000000287907733
|err| = 2.879077332096358e-8

7 — Parametrized integral

Integral depending on a parameter passed via p:

For  :  .

f7 = (u, p) -> (u[1] * u[2])^(-0.5) * exp(-p[1] * (u[1] + u[2]))
prob7 = IntegralProblem(f7, (zeros(2), ones(2)), [0.0])

for λ in (0.0, 0.5, 1.0, 2.0)
    s = solve(
        remake(prob7, p = [λ]),
        DecuhrAlgorithm(singul = 2, alpha = -0.5);
        abstol = 1.0e-8
    )
    @printf "λ = %.1f  →  I ≈ %.6f\n" λ s.u
end

λ = 0.0  →  I ≈ 4.000122
λ = 0.5  →  I ≈ 2.928494
λ = 1.0  →  I ≈ 2.231107
λ = 2.0  →  I ≈ 1.431166

8 — Automatic differentiation with ForwardDiff

The integrand is parametrized by . We compute and in forward AD mode, without finite differences.

Analytical derivative:

At  :

using ForwardDiff

f8 = (u, p) -> (u[1] * u[2])^(-0.5) * exp(-p[1] * (u[1] + u[2]))
prob8 = IntegralProblem(f8, (zeros(2), ones(2)), [0.0])

# Function I(λ). Here λ does not change the singularity structure, so we may
# either supply alpha explicitly or let it be auto-estimated — both differentiate.
I(λ) = solve(
    remake(prob8, p = [λ]),
    DecuhrAlgorithm(singul = 2, alpha = -0.5);
    abstol = 1.0e-7
).u

# First derivative
dI = ForwardDiff.derivative(I, 0.0)

# Second derivative (second-order nesting)
d2I = ForwardDiff.derivative(λ -> ForwardDiff.derivative(I, λ), 0.0)

exact_dI = -8 / 3
println("dI/dλ  ≈ ", dI, "  (exact = ", exact_dI, ")")
println("|err|    = ", abs(dI - exact_dI))
println("d²I/dλ² ≈ ", d2I)

dI/dλ  ≈ -2.666666669704931  (exact = -2.6666666666666665)
|err|    = 3.038264306809424e-9
d²I/dλ² ≈ 2.4888888906440045

The same derivative is obtained without supplying alpha: it is auto-estimated on the primal integrand (a structural property of the singularity, independent of the differentiation seed) and the integration then runs in the dual-number type.

Iauto(λ) = solve(
    remake(prob8, p = [λ]),
    DecuhrAlgorithm(singul = 2);   # alpha auto-estimated
    abstol = 1.0e-7, maxiters = 300_000
).u

println(
    "dI/dλ (auto-α) ≈ ", ForwardDiff.derivative(Iauto, 0.0),
    "   (exact = ", -8 / 3, ")"
)

dI/dλ (auto-α) ≈ -2.6666666707943176   (exact = -2.6666666666666665)

Multi-parameter gradient

We add a second parameter controlling the singularity exponent:

Analytical values at    :

f9 = (u, p) -> (u[1] * u[2])^p[2] * exp(-p[1] * (u[1] + u[2]))
prob9 = IntegralProblem(f9, (zeros(2), ones(2)), [0.0, -0.3])

# I as a function of p = [λ, μ] — singularity exponent = p[2], held fixed
Ivec(p) = solve(
    remake(prob9, p = p),
    DecuhrAlgorithm(singul = 2, alpha = -0.3);   # alpha fixed at μ₀
    abstol = 1.0e-7
).u

p0 = [0.0, -0.3]
grad = ForwardDiff.gradient(Ivec, p0)

μ = p0[2]
exact_I = 1 / (1 + μ)^2
# dI/dλ|λ=0 = -∫∫ (x₁+x₂)(x₁x₂)^μ dx = -2/((μ+1)(μ+2))
exact_dIdλ = -2 / ((1 + μ) * (2 + μ))
# dI/dμ = d/dμ [1/(1+μ)²] = -2/(1+μ)³
exact_dIdμ = -2 / (1 + μ)^3

println("I  ≈ ", exact_I)
println("∇I ≈ ", grad)
println("exact ∇I = [", exact_dIdλ, ", ", exact_dIdμ, "]")
println("|err|    = ", abs.(grad .- [exact_dIdλ, exact_dIdμ]))

I  ≈ 2.0408163265306127
∇I ≈ [-1.680672272293823, -5.830991598243292]
exact ∇I = [-1.680672268907563, -5.830903790087465]
|err|    = [3.3862599391198955e-9, 8.780815582731805e-5]

Alpha held fixed during gradient computation

When differentiating with respect to (the singularity exponent), alpha in DecuhrAlgorithm must remain a constant Float64 — it controls the extrapolation rule, not the value of the integral. For an exact gradient in , one must evaluate at and accept that the quadrature error introduces a bias of order .

9 — Budget control

What happens when the budget cannot meet the tolerance depends on how short it is, and the two cases are not the same. With a budget too small to complete even the first subdivision of a 3-D vertex singularity, the solve fails outright and there is no estimate to salvage — u comes back as zero, not as a poor approximation:

sol9 = solve(
    prob5,
    DecuhrAlgorithm(singul = 3, alpha = -1 / 3);
    abstol = 1.0e-14,   # very tight tolerance
    maxiters = 500      # far too small a budget
)

println("retcode : ", sol9.retcode)
println("u       : ", sol9.u)
println("resid   : ", sol9.resid)
println("ifail   : ", sol9.stats.ifail)

retcode : Failure
u       : 0.0
resid   : 0.0
ifail   : 13

Always test the return code before using u. With a budget large enough to refine, but still short of the requested tolerance, the return code becomes MaxIters and the value is usable — accurate to about 1 %, which is not the requested 1e-14 but is a genuine estimate:

sol9b = solve(
    prob5,
    DecuhrAlgorithm(singul = 3, alpha = -1 / 3);
    abstol = 1.0e-14,
    maxiters = 50_000
)

println("retcode : ", sol9b.retcode)
println("best estimate : ", sol9b.u)
println("|err| ≈ ", abs(sol9b.u - (3 / 2)^3))

retcode : MaxIters
best estimate : 3.401666926876933
|err| ≈ 0.02666692687693306

10 — Diagnostics: sol.stats and sol.resid

sol.stats reports the number of integrand evaluations and the raw DECUHR code, and sol.resid carries the algorithm's own error estimate. Asking the 2-D vertex singularity for 1e-12 exhausts the budget: the return code is MaxIters and the value is accurate to about 1e-4.

Note what resid says next to the true error below. The estimate is not a guaranteed bound: here it is smaller than the actual error by an order of magnitude. Treat it as an indication of how far the extrapolation has converged, not as a certificate — and where the value matters, confirm it by tightening maxiters until u stops moving.

sol10 = solve(prob, DecuhrAlgorithm(singul = 2, alpha = -0.5); abstol = 1.0e-12)

println("retcode   : ", sol10.retcode)
println("numevals  : ", sol10.stats.numevals) # integrand evaluations
println("ifail     : ", sol10.stats.ifail)    # raw DECUHR code
println("resid     : ", sol10.resid)          # the algorithm's own estimate
println("I ≈ ", sol10.u, "   true |err| = ", abs(sol10.u - 4.0))

retcode   : MaxIters
numevals  : 99905
ifail     : 1
resid     : 1.0262244253922551e-5
I ≈ 4.000121585481047   true |err| = 0.00012158548104679312

---

This page was generated using Literate.jl.