8

Writing this comment reminded me of something I noticed years ago about evaluating Bessel functions of the first kind $J_n(x)$ in Excel. (BESSELJ)

I don't use Excel now but at the time I'd checked MacOS and Windows computers with several different released versions of Excel and they all had the same problem.

This was circa 2015.

When the argument passed through $x=8$, there was a step discontinuity of order $10^{-5}$ for $n<9$.

Plotting the first difference will show that the Excel values are step-wise discontinuous at $8$.

Question: Is it possible to know which algorithm had this behavior?

Excel Bessel Function Glitch

Python source file is available.

Example from Excel:

   x         J1             J3            J5
7.9997,  0.234593634,  -0.291131046,  0.185841208,
7.9998,  0.234607873,  -0.291131434,  0.185819063,
7.9999,  0.234622108,  -0.291131819,  0.185796918,
8.0000,  0.234628387,  -0.291125242,  0.185775021,
8.0001,  0.234642617,  -0.291125622,  0.185752874,
8.0002,  0.234656846,  -0.291126,     0.185730726,
8.0003,  0.234671072,  -0.291126376,  0.185708577

edit: per request:

n    x=7.99999952316284            x=8             x=8.00000095367431
--   ------------------    ------------------      ------------------
0     0.17165091838403      0.171650807317694       0.171650583550973
1     0.2346362739686       0.234628386507074       0.234628522235498
2    -0.112991846395527    -0.112993710690925      -0.112993459984572
3    -0.291132200533783    -0.291125241852537      -0.112993459984572
uhoh
  • 1,048
  • 1
  • 10
  • 17
  • 3
    Without the actual implementation, you can only guess. My guess is, that the integer order Bessel functions are computed from $J_0, J_1$ with the recurrence relation. And (and least) in the $J_1$ implementation there is a switch between two approximations. (E.g. http://www.netlib.org/fdlibm/ has a switch at $x=8$). Do the even-order functions are 'discontinuous' too? Can you check with real order Bessel functions, e.g. $J_{3\pm0.00001}(8)?$ – gammatester Sep 28 '18 at 10:18
  • 2
    The starting point for me for computing special functions has always been the DLMF. According to that (https://dlmf.nist.gov/10.74#i), it could be a switch between a power series and an asymptotic expansion, but obviously that doesn't prove anything. – Kirill Sep 28 '18 at 10:39
  • @gammatester I'm quite rusty at excel now, but here is what I can do. I've just plotted change since 7.999997 as a rough check. Both even and odd but it stops at J9! (editing question) See https://i.stack.imgur.com/0YlZA.jpg and also https://i.stack.imgur.com/nFnzV.jpg Also, I don't see non integral orders available. – uhoh Sep 28 '18 at 11:01
  • 6
    Out of curiosity, can you try 7.999999523162841796875 and 8.00000095367431640625 (printing the output to full precision), the two Float32 numbers next to 8? – Kirill Sep 28 '18 at 18:02
  • @Kirill Thanks, I'm struggling to coax digits out of Excel. At the bottom of the question I've added some numbers that I've had to hand copy/paste one at a time to reveal digits. I don't remember how to get more. At some point we may need the help of someone who knows how to use Excel. This is Excel for Mac 2011 and the help is quite awkward. signing off for the night... – uhoh Sep 28 '18 at 18:22
  • 1
    @Kirill the jump is definitely where you thought: BESSELJ(7.99999952316284,3)=-0.291132200533783000 but BESSELJ(8.00000095367431,3)=-0.291125245492850000 forcing digits until they're all zero in excel 2013 for windows. The latter matches `BESSELJ(8,3) to 8 dp – Chris H Oct 02 '18 at 14:53
  • @ChrisH I just saw this in a news feed, I don't subscribe but the title is funny https://www.wsj.com/articles/the-first-rule-of-microsoft-exceldont-tell-anyone-youre-good-at-it-1538754380 – uhoh Oct 07 '18 at 14:01
  • 1
    @uhoh, luckily I can truthfully say that my knowledge is so outdated as to be of little use to anyone - which is good as I provide similar assistance for a lot of other software – Chris H Oct 07 '18 at 14:05

1 Answers1

4

It is clear that Microsoft Excel's implementation of the BesselJ function for $N=1$ is identical to this algorithm for arguments $X < 8$:

double ax,z;
double xx,y,ans,ans1,ans2;

if ((ax=fabs(x)) < 8.0) {
   y=xx;
   ans1=x(72362614232.0+y(-7895059235.0+y(242396853.1
         +y(-2972611.439+y(15704.48260+y(-30.16036606))))));
   ans2=144725228442.0+y(2300535178.0+y(18583304.74
         +y(99447.43394+y(376.9991397+y1.0))));
   ans=ans1/ans2;

You can implement that equation in Excel and compare it to Excel's result for BESSELJ(X,1). They are identical -- and I would note, identically erroneous for all values $0 < X < 8$, which proves Excel used this equation. For instance, compared to the exact result, the errors of BESSELJ(7,1) and $X=7$ plugged into the above are both approximately $2.24457 \times 10^{-9}$.

However, Excel encoded something else for values of $X\geq8$, and the error is considerably worse (e.g. 7.96012E-6 for $X=8$). The following code gives errors that are four orders of magnitude more accurate:

If ax=fabs(x) >= 8:

z=8.0/ax;
y=z*z;
xx=ax-2.356194491;
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
ans2=0.04687499995+y*(-0.2002690873e-3
   +y*(0.8449199096e-5+y*(-0.88228987e-6
   +y*0.105787412e-6)));
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
if (x < 0.0) ans = -ans;

I thought perhaps Excel attempted to implement the above and made a typo in one of the values, but I tried varying each of the terms in turn and was unable to replicate the values BESSELJ(X,1) computes.

Anton Menshov
  • 8,672
  • 7
  • 38
  • 94
Rob Matson
  • 41
  • 2