Dear R-Users, besselJ(pi/2, 0) # 0.4720012 besselJ(pi/2, 1E-14) # 0.4720012 besselJ(pi/2, 1E-15) # 4.720012e+14 There seems to be an error in besselJ(pi/2, eps), where eps is close to 0; although besselJ(pi/2, 0) is well behaved. I am not an expert in the field - but it doesn't seem right. Sincerely, Leonard [[alternative HTML version deleted]]
Confirmed on
R Under development (unstable) (2025-12-22 r89219)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.3 LTS
https://github.com/r-devel/r-svn/commits/50fa23f7748e9da48a15d045b6cb23a622192725/src/nmath/bessel_j.c
It's not surprising that this fails when nu is close to
.Machine$double.eps (approx 1e-16, the smallest value for which 1+x !=
1, i.e. precision limit of 64-bit floats), but it looks like one would
have to dig through the J_bessel C code
<https://github.com/r-devel/r-svn/blob/50fa23f7748e9da48a15d045b6cb23a622192725/src/nmath/bessel_j.c#L36>
to see exactly where the catastrophic loss of precision happens, and
what to do about it ...
if desired I could post this on the R bug tracker (seems appropriate).
cheers
Ben Bolker
On 12/28/25 11:28, Leo Mada via R-help wrote:> Dear R-Users,
>
> besselJ(pi/2, 0)
> # 0.4720012
> besselJ(pi/2, 1E-14)
> # 0.4720012
> besselJ(pi/2, 1E-15)
> # 4.720012e+14
>
> There seems to be an error in besselJ(pi/2, eps), where eps is close to 0;
although besselJ(pi/2, 0) is well behaved.
>
> I am not an expert in the field - but it doesn't seem right.
>
> Sincerely,
>
> Leonard
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Associate chair (graduate), Mathematics & Statistics
Director, School of Computational Science and Engineering
* E-mail is sent at my convenience; I don't expect replies outside of
working hours.
Dear R-Users, I suspected it is a catastrophic loss of precision. Unfortunately, it popped up in an unexpected way. I was plotting some curves using the Bessel function , something like this: curve(besselJ(x, pi/2-x), 0, 2*pi, ylim = c(-2,2), n = 101) abline(v = pi/2, col = "red", lty = 2) # n = 101 is the default; Note that the error is generated for a large set of values of x: median (besselJ(seq(0, 1, length.out = 1025), 1E-15)) # 9.38233e+14 Sincerely, Leonard ________________________________ From: Leo Mada <leo.mada at syonic.eu> Sent: Sunday, December 28, 2025 6:28 PM To: Leo Mada via R-help <r-help at r-project.org> Subject: Error in besselJ(pi/2, 1E-15)? Dear R-Users, besselJ(pi/2, 0) # 0.4720012 besselJ(pi/2, 1E-14) # 0.4720012 besselJ(pi/2, 1E-15) # 4.720012e+14 There seems to be an error in besselJ(pi/2, eps), where eps is close to 0; although besselJ(pi/2, 0) is well behaved. I am not an expert in the field - but it doesn't seem right. Sincerely, Leonard [[alternative HTML version deleted]]
besselJ(pi/2, sqrt(0.1)^(24:40))
[1] 4.720012e-01 4.720012e-01 4.720012e-01 4.720012e-01 4.720012e-01
1e-12 1e-13
1e-14
[6] 4.720012e-01 4.720012e-01 1.492599e+15 1.000000e+00 1.000000e+00
1e-15 1e-16
[11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
1 e-17 1e-18
1e-19
[16] 1.000000e+00 1.000000e+00
1e-20> besselJ(pi/1, 1E-15)
[1] -3.042422e+14
For what it's worth, the besselj function in Octave 8.4 doesn't do this.
octave:4> besselj(0.1.^(12:20), pi/2)
ans
0.4720 0.4720 0.4720 0.4720 0.4720 0.4720 0.4720
0.4720 0.4720
Bessel functions are most commonly wanted at integer and half-integer orders,
so it's not surprising that a problem like this could go unnoticed for a
while.
On Mon, 29 Dec 2025 at 05:29, Leo Mada via R-help <r-help at
r-project.org> wrote:>
> Dear R-Users,
>
> besselJ(pi/2, 0)
> # 0.4720012
> besselJ(pi/2, 1E-14)
> # 0.4720012
> besselJ(pi/2, 1E-15)
> # 4.720012e+14
>
> There seems to be an error in besselJ(pi/2, eps), where eps is close to 0;
although besselJ(pi/2, 0) is well behaved.
>
> I am not an expert in the field - but it doesn't seem right.
>
> Sincerely,
>
> Leonard
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.