Re: Bessel Functions and Thumbscrews

[PhilipNienhuis]

Robert T. Short wrote
> 
> On 03/11/2012 10:29 AM, PhilipNienhuis wrote:
>> Robert T. Short wrote
>>> Following up on a post and bug report several weeks ago.
> 
<snip>


>> (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.
>>
>>
> This is good information.  The octave algorithms are fast, but don't 
> even have the same precision as the A&S tables.  For my application the 
> octave routines are plenty good enough but it is a very good idea to be 
> sure you have good results before you start!
> 
> Are the A&S tables correct as far as they go?
> 
> Bob
> 

The A&S tables.... you mean on pages 416-421 (at the time I only programmed
the modified I0 I1 K0 K1 functions so that's as far as I go).
I assume that those tables are simply accurate. But mind you, they may not
be (probably aren't) based on the polynomial approximations given in A&S -
A&S themselves state that those polynomial approximations for I and K have
an error in the order of 10^-7 to 10^-8.

Inspecting the NUMAL source code reveals that many intermediary results are
computed to 14 digits precision using only simple arithmetic operations
(adding and multiplying) to retain as much precision as possible. I think it
was a very carefully thought out and optimally tuned numerical library.
While the NUMAL library was probably optimized for the CDC Cyber mainframes
of those times that used 60 bit words (and a mere 100K or so RAM), a
colleague successfully translated the I0/I1/K0/K1 Bessel functions to Turbo
Pascal which used 10-byte (80 bits) floats internally in the 387
coprocessors. His results proved more stable than mine; I suspect that
MS-Fortran only used 8 byte rather than 10 byte floats.

FYI, Octave-Forge has the gsl package which also contains Bessel functions.
I perused the relevant website but I couldn't find information on the
precision of the GSL special functions. 
The GSL refers to A&S, so I wonder if they are more accurate than Octave's
built-in versions.
So I tried (in Octave)
   exp(X)*bessel_Kn(0,X)   ## gsl pkg
and
   besselk (X, 0, "opt")      ## built-in
for X = 0.1, 1 and 10.0.
Both commands gave an answer that differed from A&S table 9.8 only in the
rightmost (11th) digit. 

>From this very limited sample I conclude that Octave's built-in and gsl's
Bessel functions (at least K0) are fairly accurate. But again, there may be
more demanding situations like I had at hand where 10 digits precision
simply isn't sufficient.

When I have time (not soon I'm afraid), I'll build a test program using my
old FORTRAN versions of the NUMAL procedures to do some more comparisons.

Philip


--
View this message in context: 
http://octave.1599824.n4.nabble.com/Bessel-Functions-and-Thumbscrews-tp4439290p4470467.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.