Hcubature_count?

Specific Domains Numerics
,
1

Is there a _count equivalent for hcubature in the same way that quadgk_count is for quadgk? Thx.

2

Nothing is currently built in. Would be easy to add since quadgk_count is just a 6-line wrapper around quadgk, i.e.

function hcubature_count(f, a, b; kws...)
    count = 0
    i = hcubature(a, b; kws...) do x
        count += 1
        f(x)
    end
    return (i..., count)
end

If someone wants to submit a PR to HCubature.jl with something like this (and a test, and docs), that would be welcome. hcubature_count · Issue #68 · JuliaMath/HCubature.jl · GitHub

3

Thx!

Not sure how to give an interpretation to the large number I obtain

# non-singular example 
hcubature(x -> x[1]/norm(x)^3, (1,1,1), (2,2,2))

(0.08562994005985379, 1.2758662568828662e-9)

# non-singular example with count 
hcubature_count(x -> x[1]/norm(x)^3, (1,1,1), (2,2,2))[3]

10659

# singular example 

hcubature(x -> x[1]/norm(x)^3, (0,0,0), (1,1,1))

(0.9693880526621552, 1.4444861385160079e-8)

# singular example with count 
hcubature_count(x -> x[1]/norm(x)^3, (0,0,0), (1,1,1))[3]

1098471

4

This is only about 100 points per dimension:

julia> cbrt(1098471)
103.18016038906488

which is quite reasonable for integrating a singular integrand to 8 digits (the default rtol).

1 Like
5

Thx again.

What do you consider to be a good measure of the degree of singularity of the integrand? Does information from truncated Taylor series map on-to-one to number of function evaluations?

I would expect 1/norm(x)^3 to be more singular than x[1]/norm(x)^3, but hcubature shows a lower function evaluation count.

The count for (x[1]*x[1])/norm(x)^3 is lower than for x[1]/norm(x)^3.

Is my understanding correct that one could output and plot the history x. Do you consider these plots valuable to include in a tutorial?

6

If it’s singular, you don’t even have a Taylor series.

If the function is analytic (has a convergent Taylor series in a neighborhood of the integration domain), then convergence is determined by things like the closest singularity/pole in the complex plane and/or how oscillatory the function is.

If it’s not analytic, then convergence is usually determined by the biggest “non-smoothness”. For example, the lowest derivative that is discontinuous. Or, if it has an integrable singularity (a point where the integrand blows up in the domain, but slowly enough that the integral is still finite), then by the fastest rate of blowing up.

That’s because for 1/norm(x)^3 it simply gives up without giving a finite answer: this is not an integrable singularity in 3d (the integral is \infty).

julia> hcubature_count(x -> 1/norm(x)^3, (0,0,0), (1,1,1)) # non-integrable singularity
(Inf, NaN, 67287)

julia> hcubature_count(x -> x[1]/norm(x)^3, (0,0,0), (1,1,1)) # integrable singularity
(0.9693880526621552, 1.4444861385160079e-8, 1098471)
7

Many thanks for your input. Much appreciated.