# How to improve the implementation of a function involving a numerical integration?

**URL:** <https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444>\
**Category:** Performance\
**Tags:** quadgk, integral\
**Created:** [October 26, 2023, 3:37pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444 "2023-10-26T15:37:24Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [October 27, 2023, 12:25am UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/6 "2023-10-27T00:25:44Z")

</div>

> [@stevengj](#):
>
> ```julia
> julia> quadgk_count(θ -> k(0.1,cos(θ),1e1) * sin(θ)^2/π, 0,pi, rtol=1e-4)
> (0.4993508746129765, 6.152897258360213e-8, 15)
> 
> ```

The change of variables trick works about as well as the specialized quadrature for Jacobi weight functions with QuadJGK.jl. I get:

```julia
julia> using QuadJGK
[Info: Precompiling QuadJGK [f27e79a0-bbaa-47b5-b2d3-8db3cef8b089]

julia> k(x,y,l) = exp(-(x-y)^2/(2*l^2))
k (generic function with 1 method)

julia> p(x) = sqrt(1-x^2)/π
p (generic function with 1 method)

julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)))
(0.49935087461297667, 1.1102230246251565e-16, 15)

```

For this example, machine precision is still possible with 7 integrand evaluations

```julia
julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=3)
(0.4993508746129767, 1.6239870559431324e-10, 7)

```

And the same is true for the second integral

```julia
julia> ES(x,l) = quadjgk(y -> k(x,y,l) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)))[1]
ES (generic function with 1 method)

julia> quadjgk_count(x -> p(x) * ES(x,1e1), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=3)
(0.24937694744330427, 8.052006283953972e-11, 7)

```

So `7^2 = 49` function evaluations should be sufficient to compute the normalization for this large value of `l=1e1`. For \ell \ll 1, HCubature.jl is probably the better option for the 2d integral since the integrand evaluations scale as \mathcal{O}(\log(\ell)) as \ell \to 0 for cubatures of a Gaussian localized at the origin, whereas for iterated integration it will be \mathcal{O}(\log^2(\ell)) (i.e. one factor for each integral).

Update: by loosening the error tolerance, it seems like only 5 integrand evaluations are needed for each 1d integral, and the result is still accurate to machine precision

```julia
julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=2, rtol=1e-4)
(0.4993508746129768, 3.8987719608307714e-7, 5)

julia> quadjgk_count(x -> p(x) * ES(x,1e1), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=2, rtol=1e-4)
(0.24937694744330427, 1.9377299406708737e-7, 5)

```

Edit: I wasn’t thinking: k(x,y;l) is localized on the line x=y, not at a point, so the scaling of cubature is like \mathcal{O}(\log(l)/l), which makes iterated integration more efficient as l \to 0. It would be hard to get rid of the logarithmic scaling without doing something special

---

_[View the full topic](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444)._