Re: Observed large error from spherical Bessel functions at high orders

[Patrick Alken]

Hello,

  Unfortunately some of the special functions do not produce the
expected accuracy for all inputs, and it seems you've found another
example. It would be great if you could look into algorithms for
computing these functions and possibly find a way to improve the
accuracy for these inputs. Unfortunately we don't have a lot of
developers working on GSL these days, and I don't have much expertise in
computing the bessel functions.

Please at least log this issue in the bug tracker:
http://savannah.gnu.org/bugs/?group=gsl

That way it won't be lost and someone may eventually look into it.

Thanks,
Patrick

On 7/14/20 11:05 AM, Ziqi Fan wrote:
> Dear developers and maintainers of GNU Scientific Library,
>
> My name is Ziqi Fan. I am a PhD candidate from University of Florida. I
> have been using GSL for at least 3 years for my numerical solver project.
> Recently, I met a bug in my solver and the bug should not exist from a
> theoretical perspective. So I checked very carefully my own implementation
> and found that the bug originated from using the routine for spherical
> bessel functions of the second kind: gsl_sf_bessel_yl.
>
> For instance, I tested the function using a large order n = 45, and got the
> following result: y_45(6.975948) = -726330209582507479571265224704.000000,
> err = 19103761820604644.000000, where the first term is the output and the
> second term is an estimated error provided by the routine
> "gsl_sf_bessel_yl_e". The error is relatively small compared to the output,
> but it can easily cause divergence of a numerical algorithm, considering
> the magnitude of its value.
>
> I once implemented spherical bessel functions myself using the numerical
> recipe, and I understand that large error occurs when x << n, where x is a
> positive input to the function, and n is the order of the function. I am
> wondering if it is possible to resolve the large absolute error at high
> order of the spherical bessel functions of the second kind. If I have an
> opportunity to communicate with an engineer or mathematician from you, I
> may be able to help contribute a routine with a higher precision. If so, a
> new routine with a higher precision can really contribute to the invention
> of many new algorithms for powerful numerical solvers.
>
> I really appreciate your effort in providing and maintaining this great
> software, and I look forward to hearing from you!
>
> Best,
> Ziqi