How to avoid catastrophic cancellation Bessel Integrals

General Usage
1

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
image873×145 13.7 KB

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.

2

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.

3

@RayleighLord I know little about Bessel functions but arb knows more.
You can access those functions throught Arblib.jl:

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)

2 Likes