HCubature convergence

General Usage
1

I do not understand why HCubature.jl v1.6.0 makes two different estimates of the error E (HCubature.jl/src/HCubature.jl at master · JuliaMath/HCubature.jl · GitHub) for the following two versions of the integrand:

using HCubature

function singlet(f)
    Δf² = 0.5^2
    100 * Δf² / (Δf² + 4*f^2)
end

normInf(x) = norm(x, Inf)

freqs = -2:0.01:2
δf, res = zeros(length(freqs)), zeros(length(freqs))

function integr(ppt)
    ρ, z = ppt
    δf .= freqs .- (z^2  + 0.5z - 0.5ρ^2)
    res .= 2π * ρ * singlet.(δf)
    res
end

function integr2(ppt)
    ρ, z = ppt
    δf = freqs .- (z^2  + 0.5z - 0.5ρ^2)
    res = 2π * ρ * singlet.(δf)
    res
end

a, b = (0, -1.2), (1, +1.2)
r1, _ = hcubature(integr, a, b; rtol=1e-2, atol=0.5,
                   maxevals=100_000, norm=normInf)

r2, _ = hcubature(integr2, a, b; rtol=1e-2, atol=0.5,
                   maxevals=100_000, norm=normInf)

The integration with integr terminates when maxeval=100_000 is reached, integr2 requires about 1000 evaluations.

2

I do not know the internals of hcubature, but in integr you are reusing the output vector. It could be that hcubature saves the output, without copying it, and does some comparisons or something across iterations. But then the content of the saved output is overwritten by integr; confusion ensues. If you return copy(res) from integr, it works better.

And, the global variables should be declared as const, otherwise integr must do some type checking on them every time.

1 Like
3

@sgaure is correct that to get the correct answer a new array must be allocated by every integrand evaluation and this is because HCubature.jl assumes that the integrand returns an immutable type. A similar discussion occurred here: Broken inplace integrand for HCubatureJL · Issue #169 · SciML/Integrals.jl · GitHub

4

Thank you, I am aware of that. The code used to be within a function.