2
$\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
$\endgroup$
  • 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 Sep 28 '18 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 Sep 28 '18 at 10:39
  • 4
    $\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 Sep 28 '18 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 Oct 2 '18 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 Oct 7 '18 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.