# HCubature convergence

**URL:** <https://discourse.julialang.org/t/hcubature-convergence/111840>\
**Category:** General Usage\
**Tags:** hcubature\
**Created:** [March 19, 2024, 8:51pm UTC](https://discourse.julialang.org/t/hcubature-convergence/111840 "2024-03-19T20:51:29Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [March 19, 2024, 8:51pm UTC](https://discourse.julialang.org/t/hcubature-convergence/111840/1 "2024-03-19T20:51:29Z")

</div>

I do not understand why [`HCubature.jl`](https://github.com/JuliaMath/HCubature.jl) `v1.6.0` makes two different estimates of the error `E` ([HCubature.jl/src/HCubature.jl at master · JuliaMath/HCubature.jl · GitHub](https://github.com/JuliaMath/HCubature.jl/blob/master/src/HCubature.jl#L158)) for the following two versions of the integrand:

```julia
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.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 19, 2024, 9:37pm UTC](https://discourse.julialang.org/t/hcubature-convergence/111840/2 "2024-03-19T21:37:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [March 19, 2024, 10:28pm UTC](https://discourse.julialang.org/t/hcubature-convergence/111840/3 "2024-03-19T22:28:14Z")

</div>

@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](https://github.com/SciML/Integrals.jl/issues/169)

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [March 20, 2024, 7:22am UTC](https://discourse.julialang.org/t/hcubature-convergence/111840/4 "2024-03-20T07:22:29Z")

</div>

> [@sgaure](#):
>
> And, the global variables should be declared as `const`, otherwise `integr` must do some type checking on them every time.

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