[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Bessel Functions and Thumbscrews
From: |
PhilipNienhuis |
Subject: |
Re: Bessel Functions and Thumbscrews |
Date: |
Sun, 11 Mar 2012 10:29:20 -0700 (PDT) |
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.