8
$\begingroup$

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
Improve this question
$\endgroup$
8
  • 3
    $\begingroup$ 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. 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)?$ $\endgroup$
    – gammatester
    Commented Sep 28, 2018 at 10:18
  • 2
    $\begingroup$ The starting point for me for computing special functions has always been the DLMF. According to that (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. $\endgroup$
    – Kirill
    Commented Sep 28, 2018 at 10:39
  • 6
    $\begingroup$ Out of curiosity, can you try 7.999999523162841796875 and 8.00000095367431640625 (printing the output to full precision), the two Float32 numbers next to 8? $\endgroup$
    – Kirill
    Commented Sep 28, 2018 at 18:02
  • 1
    $\begingroup$ @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 $\endgroup$
    – Chris H
    Commented Oct 2, 2018 at 14:53
  • 1
    $\begingroup$ @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 $\endgroup$
    – Chris H
    Commented Oct 7, 2018 at 14:05

1 Answer 1

Reset to default
4
$\begingroup$

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=x*x;
   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+y*1.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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.