Robert T. Short wrote
Following up on a post and bug report several weeks ago.
I have attached a patch in which I added a bunch of tests to
bessel{j,y,i,k}. The existing tests were fine, but the Amos code
(libcruft/amos) uses different algorithms for different orders and
arguments and the existing tests were for small values of the order and
magnitude. I created a subset of the tables in Abramowitz and Stegun
and compared the octave results to the tables. In particular for large
order/arguments Amos used the asymptotic expansion and so one might
expect some differences (and there are).
I also added tests to cover the bug that set me off on this mission -
checking that negative integer orders give the correct answers.
Bob
(Sorry for reacting a bit late)
The A&S Bessel functions are fairly "rough" approximations.
Some years back I used Bessel functions in an numerical inverse Laplace
transform (Stehfest) and it turned out that an accuracy of at least 10-12
significant digits was needed to get some stability. The A&S approximations
were insufficiently accurate by far.
I my programs of the time (some 15+ years back) I used (FORTRAN) versions
for Ix and Kx that I based on ALGOL functions from the NUMAL library (CWI,
Amsterdam) from decades ago. I only used real arguments, but perhaps the
principle works for complex arguments as well.
AFAICR my (those) FORTRAN versions were quite a bit more involved than
earlier ones I based on A&S and thus a bit slower (containing a series of
IFs depending on argument magnitude), but MUCH more accurate. A bit of a
trade-off.
Perhaps this NUMAL library can be used for Octave as well. It seems to live
here these days:
http://www.softwarepreservation.org/projects/ALGOL/applications
and the license seems a bit BSD like.
Any license gurus who can shed a light here?
Philip
--
View this message in context:
http://octave.1599824.n4.nabble.com/Bessel-Functions-and-Thumbscrews-tp4439290p4464264.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.