# HCubature convergence

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.

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

@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

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