Convergence of hcubature

Specific Domains Numerics
1

I am trying to integrate a Lorentzian offset by a smooth function which is defined on a 3D polar domain. Here is the code:

using Cubature

lorentzian_re(freq, ampl, freq0, fwhm) = ampl ./ (1 + ((freq - freq0) / fwhm).^2)
singlet_0_1Hz(freq) = lorentzian_re(freq, 100.0, 0.0, 0.1)
freq_range = linspace(-2.5, +2.5, 101)

function lineshape(reltol, abstol, f=freq_range)
    ρ0, z0 = 2e-1, 8e-1
    field = p -> 5 * p[1] * p[3] * sin(p[2])  # hcubature chokes on this function
    #field = p -> 5 * p[1] * p[3] * cos(p[2]) # hcubature integrates this without any problem
    #field = p -> 5 * p[1] * sin(p[2]) # no problem with this one either
    inv_vol = 1.0 / (2π * ρ0^2 * z0)

    integr(polar_point, v) = begin
        # integral in polar coordinates is rho * f(rho, phi, z), normalized by the integration volume
        ρ_norm = polar_point[1] * inv_vol
        v[:] = ρ_norm * singlet_0_1Hz(f - field(polar_point))
    end

    hcubature(length(f), integr, [0.0, 0.0, -z0], [ρ0, 2π, +z0], reltol=reltol, abstol=abstol)
    #pcubature(length(f), integr, [0.0, 0.0, -z0], [ρ0, 2π, +z0], reltol=reltol, abstol=abstol)
end

The integrand is a smooth function and the convergence should not be a problem. It is not a problem when the function field is p[1]*p[3]*cos(2p[2]), but hcubature runs forever on the very similar function p[1]*p[3]*sin(2p[2]), with the cosine replaced by the sine. pcubature has no problem with either functions.

Does anybody have any idea why?

michele

2

If the integral is zero (e.g. because you are integrating an odd function), then you will need to set a nonzero abstol as explained in the Cubature.jl README.

3

Sorry for not being clearer: the integrand is not an odd function.

The integrand is the nonnegative lorentzian function 100 / (1 + 100*f^2) where f is either
f = linspace(-2.5, 2.5, 101) - ρ*z*cos(ϕ) or
f = linspace(-2.5, 2.5, 101) - ρ*z*sin(ϕ);
and because the integration is done in polar coordinates, there is the additional volume factor ρ≧0.

(You may recognize ρ*z*cos(ϕ) as the real solid harmonic n=2, m=1 written in polar coordinates; the sine term is for m=-1).

What puzzles me is that the integral containing sin(ϕ) and that containing cos(ϕ) should give the same result (the integrand being a periodic function integrated over the period 2π), but even pcubature gives two different results using reltol=1e-5, abstol=1e-3. (Left is the cosine term integral).

I tried the integrand f -> f^2 for which there is an analytical solution, and the numerical results for the sine and the cosine terms are the same.

4

I am at loss here: I am integrating the functions ρ*cos(ϕ) and ρ*sin(ϕ) on a disk using pcubature and hcubature and obtain two different results that are not within the error (top is integrals, bottom is errors)

Here is the code: the integrand is the square function.

using Cubature
using Plots
pyplot()

function test_cuba_method(field, reltol, abstol, f, method)
    ρ0 = 2e-1
    inv_vol = 1.0 / (π * ρ0^2)

    integr(polar_point, v) = begin
        ρ_norm = polar_point[1] * inv_vol
        v[:] = ρ_norm * (f - field(polar_point)).^2
    end

    method == :hcuba && return hcubature(length(f), integr, [0.0, 0.0], [ρ0, 2π], reltol=reltol, abstol=abstol)
    method == :pcuba && return pcubature(length(f), integr, [0.0, 0.0], [ρ0, 2π], reltol=reltol, abstol=abstol)
end


function test()
    freq_range = linspace(-0.5, +0.5, 101)
    field_sin = p -> 5 * p[1] * sin(p[2])
    field_cos = p -> 5 * p[1] * cos(p[2])
    @time line_cm1, err_cm1 = test_cuba_method(field_cos, 1e-6, 1e-4, freq_range, :hcuba)
    @time line_cm2, err_cm2 = test_cuba_method(field_sin, 1e-6, 1e-4, freq_range, :hcuba)
    @time line_cm3, err_cm3 = test_cuba_method(field_cos, 1e-6, 1e-4, freq_range, :pcuba)
    @time line_cm4, err_cm4 = test_cuba_method(field_sin, 1e-6, 1e-4, freq_range, :pcuba)

    plot(plot(freq_range, hcat(line_cm1, line_cm2, line_cm3, line_cm4),
              label=["cos, hcuba", "cos, hcuba", "cos, hcuba", "cos, hcuba"]'),
         plot(freq_range, hcat(err_cm1, err_cm2, err_cm3, err_cm4),
              label=["cos, hcuba", "cos, hcuba", "cos, hcuba", "cos, hcuba"]'),
              reuse=false, layout=(2,1))
end
5

I can’t reproduce your problem (Julia 0.5 on MacOS).

julia> hcubature(x -> x[1] * (0.3 - 5*x[1]*sin(x[2])).^2, [0,0], [0.2,2π], reltol=1e-6, abstol=1e-4)
(0.042725660088821185,1.348693715190561e-5)

julia> hcubature(x -> x[1] * (0.3 - 5*x[1]*cos(x[2])).^2, [0,0], [0.2,2π], reltol=1e-6, abstol=1e-4)
(0.042725660088821185,1.3486937151908646e-5)

julia> freq_range = linspace(-0.5, +0.5, 4);
julia> field_sin = p -> 5 * p[1] * sin(p[2]);
julia> field_cos = p -> 5 * p[1] * cos(p[2]);

julia> test_cuba_method(field_cos, 1e-6, 1e-4, freq_range, :hcuba)
([0.5,0.277778,0.277778,0.5],[8.94878e-5,8.83078e-5,8.71277e-5,8.59476e-5])

julia> test_cuba_method(field_sin, 1e-6, 1e-4, freq_range, :hcuba)
([0.5,0.277778,0.277778,0.5],[8.94878e-5,8.83078e-5,8.71277e-5,8.59476e-5])

i.e. I appear to be getting the same results for cos and sin, both with your posted test_cuba_method code and with a simplified 1-output case.

6

Could you please run test_cuba_method with field_sin and the :pcuba switch? That is where I see the problem with pcubature.

This is Julia 0.5 on Windows 7.

7

Okay, I can reproduce your problem:

julia> pcubature(x -> x[1] * (0.3 - 5*x[1]*cos(x[2])).^2, [0,0], [0.2,2π], reltol=1e-6, abstol=1e-4)
(0.04272566009868226,9.74283563556777e-6)

julia> pcubature(x -> x[1] * (0.3 - 5*x[1]*sin(x[2])).^2, [0,0], [0.2,2π], reltol=1e-6, abstol=1e-4)
(0.01130973355292325,3.469446951953614e-18)

pcubature starts with a low-order cubature rule and then progressively doubles the order until convergence. The problem is that for an interval [0,2π], the lowest-order rule that it starts with evaluates the function at only 3 points: 0, π, and 2π. At all three of these points, the sin function is zero, so it thinks it is just integrating the linear function x[1]*0.3, which is integrated exactly by the 3-point rule, and hence it stops.

(You can see this just by putting a println statement in the integrand to see what points it is evaluating.)

At some point, I thought about having an option to force pcubature to start with a higher-order rule, but never implemented it. In your case, however, the workaround is easy: just use cosine, or use the sine function with some arbitrary phase shift.

8

Thank you. I also forgot to report that the execution time for the integral containing the sin function was suspiciouly low.

My original question was about hcubature iterating forever on the following code:
hcubature(p -> 100 * p[1] / (1 + 20 * (0.2 - 5 * p[1] * p[3] * sin(p[2])) .^ 2), [0.0, 0.0, -0.8], [0.2, 2π, 0.8], reltol=1e-6, abstol=1e-4)
but converging quickly on
hcubature(p -> 100 * p[1] / (1 + 20 * (0.2 - 5 * p[1] * p[3] * cos(p[2])) .^ 2), [0.0, 0.0, -0.8], [0.2, 2π, 0.8], reltol=1e-6, abstol=1e-4)
A println statement does not reveal anything particular (at least to me).

They do converge to numbers differing by more than the error estimate when the tolerances are reduced to reltol=1e-3, abstol=1e-1.

9

Is there a way of specifying that the integrand is periodic in one of the integration variables so that the integration points are equally spaced? In the example above, pcubature evaluates the integrand at two independent points, [0,π]; the 2π point is wasted computation.

Or do you have any plans to add that?

10

No. If you know what order of quadrature rule you need in the periodic direction, you can just implement that quadrature yourself with equally-spaced points (so that you only do a 1d adaptive quadrature with pquadrature or hquadrature.

11

If you know what order of quadrature rule you need…

Isn’t this dependent on what the function does in the other variables? Do I figure that out by doing what Cubature does, computing in all variables, checking for convergence and refining?

You assume I know the order as if this is straightforward, which implies that I am missing the obvious.