Incorrect and variable results for smooth function depending on initdiv

[johnbcoughlin]

I'm getting inconsistent and incorrect results for integrating a very smooth function:

julia> f = v -> begin
                  c2 = norm(v)^2
                  c2 * pi^(-3/2) * exp(-c2)
              end
#3402 (generic function with 1 method)

julia> a
3-element Array{Float64,1}:
 -100.0
 -100.0
 -100.0

julia> b
3-element Array{Float64,1}:
 100.0
 100.0
 100.0

julia> hcubature(f, a, b, atol=0, initdiv=1)
(0.0, 0.0)

julia> hcubature(f, a, b, atol=0, initdiv=2)
(0.18749999994468158, 2.7930336351243728e-9)

julia> hcubature(f, a, b, atol=0, initdiv=3)
(1.4999999999365934, 2.235159620943183e-8)

julia> hcubature(f, a, b, atol=0, initdiv=4)
(0.18749999994468125, 2.7930336351226787e-9)

The correct answer is 1.5, as can be verified with Mathematica. Perhaps this has something to do with the function being radially symmetric? It also probably has to do with the integrand being zero (to within machine precision) at the boundaries. If I use boundaries at (-4, 4) instead of (-100, 100), I get the right answer.

[stevengj]

It also probably has to do with the integrand being zero (to within machine precision) at the boundaries.

Yup. Your integrand underflows to 0.0 for |x| > 30 or so. So, if you have too few initial points in a subregion then the quadrature rule won't be able to distinguish it from a function that is identically zero in that region, and won't refine.

I'm guessing you really want to simulate an infinite integration domain? It is better to do that with a coordinate transformation than by using a huge box. Otherwise, you simply need to set a large enough initdiv.