JuliaMath / Cubature.jl

One- and multi-dimensional adaptive integration routines for the Julia language
Other
126 stars 21 forks source link

hcubature hangs on some integrand #25

Closed mzaffalon closed 3 years ago

mzaffalon commented 7 years ago

hcubature seems to hang for some integrands when the integrand is symmetric in the domain of integration and the tolerance is sufficiently high.

julia> hcubature(x -> (1.0+(x[1]*x[3]*cos(x[2]))^2), [0,0,-0.2], [0.2,2π,0.2], abstol=1e-6)
(0.5026995050032178,3.3306690738754696e-16)
julia> hcubature(x -> (1.0+(x[1]*x[3]*sin(x[2]))^2), [0,0,-0.20001], [0.2,2π,0.2], abstol=1e-6) # note the domain of integration
(0.502712074725032,4.440892098500626e-16)
julia> hcubature(x -> (1.0+(x[1]*x[3]*sin(x[2]))^2), [0,0,-0.2], [0.2,2π,0.2], abstol=1e-6) # hangs forever
julia> hcubature(x -> (1.0+(x[1]*x[3]*sin(x[2]))^2), [0,0,-0.2], [0.2,2π,0.2], abstol=1e-3) # note the reduced tolerance
(0.5026887210432646,3.325688687261241e-5)

(Sorry for the contrived example: I could not find one with lower dimensionality.)

mzaffalon commented 3 years ago

The same integral does not hang with HCubature:

using HCubature
julia> @time hcubature(x -> (1.0+(x[1]*x[3]*sin(x[2]))^2), [0,0,-0.2], [0.2,2π,0.2], atol=1e-6)
  0.051035 seconds (200.74 k allocations: 10.392 MiB)
(0.5026995050032179, 0.0)

EDIT: ... because it was fixed in https://github.com/JuliaMath/HCubature.jl/issues/4#issuecomment-319518345

stevengj commented 3 years ago

Yes, same as JuliaMath/HCubature.jl#4 and stevengj/cubature#7. This needs to be fixed in the underlying C library.

stevengj commented 3 years ago

Should be fixed by https://github.com/JuliaRegistries/General/pull/25426.

stevengj commented 3 years ago

Yup, it works now:

julia> Cubature.hcubature(x -> (1.0+(x[1]*x[3]*sin(x[2]))^2), [0,0,-0.2], [0.2,2π,0.2], abstol=1e-6)
(0.5026995050032179, 1.1102230246251565e-16)