Slow cubature

Ok, I found the old (from 2015) code where I used xgiven, it was in Fortran, not even in C.

Consider the integral between 1 and 0 of the Gaussian function

f(x) = \exp\left(-\left(\frac{x - x_0}{\mu}\right) ^ 2\right)

The value of the integral is

\frac{\sqrt{\pi} \mu}{2} \left(\mathrm{erf}\left(\frac{1 - x_0}{\mu}\right) + \mathrm{erf}\left(\frac{x_0}{\mu}\right)\right)

Let’s compute the expected result with SpecialFunctions and the parameters (x_0 = 0.5, \mu = 0.0001):

julia> using SpecialFunctions

julia> expected(x₀, μ) = 1/2 * sqrt(pi) * μ * (erf((1 - x₀) / μ) + erf(x₀ / μ))
expected (generic function with 1 method)

julia> expected(0.5, 0.0001)
0.0001772453850905516

Now let’s compute the integral numerically with Cuba.jl:

julia> using Cuba

julia> f(x, x₀=0.5, μ=0.0001) = exp(-((x - x₀) / μ) ^ 2)
f (generic function with 3 methods)

julia> integrand(x, y) = y[1] = f(x[1])
integrand (generic function with 1 method)

julia> divonne(integrand)
Component:
 1: 8.862365018700048e-5 ± 1.579627496207042e-10 (prob.: 0.0)
Integrand evaluations: 5257
Number of subregions:  40
Note: The desired accuracy was reached

julia> divonne(integrand; ngiven=1, ldxgiven=1, xgiven=[0.5;;])
Component:
 1: 0.0001772444738865438 ± 3.5121869189594667e-10 (prob.: 0.0)
Integrand evaluations: 9217
Number of subregions:  70
Note: The desired accuracy was reached

julia> cuhre(integrand)
Component:
 1: 2.608731354364103e-6 ± 2.135531666691244e-10 (prob.: 1.0)
Integrand evaluations: 35165
Number of subregions:  271
Note: The desired accuracy was reached

julia> cuhre(integrand; rtol=1e-8)
Component:
 1: 0.00017723214024825766 ± 1.7700547231374744e-12 (prob.: 1.0)
Integrand evaluations: 150865
Number of subregions:  1161
Note: The desired accuracy was reached

julia> cuhre(integrand; rtol=1e-12)
Component:
 1: 0.0001772412672726941 ± 9.869539325496199e-13 (prob.: 1.0)
Integrand evaluations: 154375
Number of subregions:  1188
Note: The desired accuracy was reached

The most accurate result is with divonne when passing the position of the peaks, which takes more evaluations than when not passing the peaks (but the result is also off from the correct result, although within the default, but quite low, tolerance), but it takes orders of magnitude fewer evaluations than cuhre, which I believe uses a similar algorithm to HCubature.