# Convergence of hcubature

**URL:** <https://discourse.julialang.org/t/convergence-of-hcubature/660>\
**Category:** Numerics\
**Created:** [November 30, 2016, 9:12pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660 "2016-11-30T21:12:16Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [November 30, 2016, 9:12pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/1 "2016-11-30T21:12:16Z")

</div>

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

```julia
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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 30, 2016, 9:33pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/2 "2016-11-30T21:33:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 1, 2016, 1:25pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/3 "2016-12-01T13:25:08Z")

</div>

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).

![](https://global.discourse-cdn.com/julialang/original/3X/0/7/07d0fc72bd1bead2d4f8c438e0855ea50f383914.png) ![](https://global.discourse-cdn.com/julialang/original/3X/2/a/2ad1484b6134a28a5f6d63e34456ebb8a237afde.png)

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.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 1, 2016, 1:45pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/4 "2016-12-01T13:45:38Z")

</div>

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)

![](https://global.discourse-cdn.com/julialang/original/3X/0/8/08c2b45d69d1560ad2d3c41f9c6ab4fa30ab67b2.png)

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 1, 2016, 2:47pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/5 "2016-12-01T14:47:17Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 1, 2016, 2:56pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/6 "2016-12-01T14:56:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 1, 2016, 11:34pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/7 "2016-12-01T23:34:43Z")

</div>

Okay, I can reproduce your problem:

```julia
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.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 2, 2016, 6:37am UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/8 "2016-12-02T06:37:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 2, 2016, 10:16am UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/9 "2016-12-02T10:16:54Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 2, 2016, 12:24pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/10 "2016-12-02T12:24:37Z")

</div>

> [@mzaffalon](#):
>
> Is there a way of specifying that the integrand is periodic in one pof the integration variables so that the integration points are equally spaced?

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`.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [December 2, 2016, 1:38pm UTC](https://discourse.julialang.org/t/convergence-of-hcubature/660/11 "2016-12-02T13:38:43Z")

</div>

> 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.