# How to avoid catastrophic cancellation Bessel Integrals

**URL:** <https://discourse.julialang.org/t/how-to-avoid-catastrophic-cancellation-bessel-integrals/124363>\
**Category:** General Usage\
**Tags:** question\
**Created:** [January 2, 2025, 12:41pm UTC](https://discourse.julialang.org/t/how-to-avoid-catastrophic-cancellation-bessel-integrals/124363 "2025-01-02T12:41:22Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![RayleighLord](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rayleighlord/32/33592_2.png) [@RayleighLord](https://discourse.julialang.org/u/RayleighLord)\
**Post date:** [January 2, 2025, 12:41pm UTC](https://discourse.julialang.org/t/how-to-avoid-catastrophic-cancellation-bessel-integrals/124363/1 "2025-01-02T12:41:22Z")

</div>

Recently, I needed to evaluate some integrals involving linear combinations of J\_m and Y\_m. There are well-known analytical expressions for those, which are

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/b/db4c2ee01d425bee9abd450aa44c70cacacdd0e5.png)

However, I am finding that when \alpha and \beta are very close, and for certain m, like m=7, I think I am getting catastrophic cancellation since the results are completely wrong. A comparison with `QuadGK` confirms that indeed the results are wrong.  
Also, the results for \alpha \approx \beta should get really close to the \alpha = \beta case, which should be 1 since I am working with an orthonormal basis.

What could be a good strategy to alleviate this issue? I would like to keep using the analytical expression since it is much faster.

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [January 2, 2025, 1:12pm UTC](https://discourse.julialang.org/t/how-to-avoid-catastrophic-cancellation-bessel-integrals/124363/2 "2025-01-02T13:12:07Z")

</div>

Have you tried integrating over the difference δ = α-β starting from zero? Since you know δ = 0 you can write the result as a definite integral… sometimes this works. At least you will get a taylor expansion which may be accurate at the scale when the cancellation becomes problematic.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [January 5, 2025, 11:52pm UTC](https://discourse.julialang.org/t/how-to-avoid-catastrophic-cancellation-bessel-integrals/124363/3 "2025-01-05T23:52:51Z")

</div>

@RayleighLord I know little about Bessel functions but [`arb` knows more](https://arblib.org/acb_hypgeom.html#bessel-functions).  
You can access those functions throught Arblib.jl:

```julia
julia> Arblib.hypgeom_bessel_
hypgeom_bessel_i! hypgeom_bessel_i_0f1! hypgeom_bessel_i_asymp! hypgeom_bessel_i_integration!
hypgeom_bessel_i_scaled! hypgeom_bessel_j! hypgeom_bessel_j_0f1! hypgeom_bessel_j_asymp!
hypgeom_bessel_jy! hypgeom_bessel_k! hypgeom_bessel_k_0f1! hypgeom_bessel_k_0f1_series!
hypgeom_bessel_k_asymp! hypgeom_bessel_k_integration! hypgeom_bessel_k_scaled! hypgeom_bessel_y!

```

you will get certified bounds (note that QuadGK can report tight error estimate with very wrong answer for pathological functions)