Airy function wrong

[Alex van der Spek]

Just finished compiling the NetLib FORTRAN routines from Donald Amos for the
Airy function. This triggered me to see if gnuplot had the Airy function. It
does! Unfortunately:

gnuplot> print airy(0.1)
0.32848812610163

The correct value according to Abramowitz Stegun is
0.32920313

The Amos routines return more digits
0.329203129943538

There must be a glaring error in the gnuplot implementation for such a large
discrepancy to occur. If I can be of help to sort out where this might be I
will be glad to do so.

Alex van der Spek

[sfeam]

The gnuplot implementation is documented as follows:
/* ------------------------------------------------------------
Airy Function Ai(x)

After:
"Two-Point Quasi-Fractional Approximations to the Airy Function Ai(x)"
by Pablo Martin, Ricardo Perez, Antonio L. Guerrero
Journal of Computational Physics 99, 337-340 (1992)

Beware of a misprint in equation (5) in this paper: The second term in
parentheses must be multiplied by "x", as is clear from equation (3)
and by comparison with equation (6). The implementation in this file
uses the CORRECT formula (with the "x").

This is not a very high accuracy approximation, but sufficient for
plotting and similar applications. Higher accuracy formulas are
available, but are much more complicated (typically requiring iteration).

Added: janert (PKJ) 2009-09-05
------------------------------------------------------------ */

Figure 1a of the 1992 paper shows a maximum error between 0.003 and 0.012
in the region of x=0.1 So your result (delta = 0.007) is consistent
with the accuracy claimed for the algorithm by the original authors .

If you need a more accurate implementation and are willing to code it,
I'm sure it would be welcomed. I suspect that might require first
implementing a more general routine for Bessel functions; so far
gnuplot only handles J0 and J1.

Ethan

[Alex van der Spek]

Thank you Ethan,

Not much coding left to do.

I have all of the Amos FORTRAN routines from Netlib
(http://netlib.org/amos/) compiled to a windows DLL. All in double precision
and allowing complex input.
This covers the Bessel function J, Y, I and K (of integer and fractional
order), the Airy functions Ai and Bi and the log of the gamma function.
Should be with 15 digit precision.

I have all of the Cody FORTRAN routines from Netlib
(http://www.netlib.org/specfun/) compiled to a windows DLL too. All in
double precision for real inputs only.
This covers the Bessel functions J0, J1, Y0, Y1, K0, K1, I1 and I2, plus J,
Y, K and I of integer or fractional order up to order 10 and the psi
function, gamma function, the log thereof, Dawson's integral and the Airy
functions Ai and Bi as well as the spherical Bessel functions jn and yn up
to order 10.

I would need some handholding to contribute this to gnuplot. Perhaps you can
get me started with that? I would not even know where to begin to find the
proper documentation. Besides C is not my strongest language. Compiling to a
Linux shared object is probably less of an issue and I would assume that .so
file would also work on a mac?

By the way, my comments were in no way meant as criticism. I am painfully
aware that coding errors are part of life. I am even more painfully aware
that those coding errors that do not generate a 'bug' can remain dormant for
quite a while. Thus I thought that in case there was a dormant coding error
it would be easy to find. I wasn't aware of the algorithm you mention but I
can't see why anyone would want to use it as there are far better algorithms
available.

Regards,
Alex




"sfeam" <sf...@users.sourceforge.net> wrote in message
news:k01e9m$sr4$1...@dont-email.me...

[sfeam]

Alex van der Spek wrote:

> Thank you Ethan,
>
> Not much coding left to do.
>
> I have all of the Amos FORTRAN routines from Netlib
> (http://netlib.org/amos/) compiled to a windows DLL. All in double
> precision and allowing complex input.
> This covers the Bessel function J, Y, I and K (of integer and fractional
> order), the Airy functions Ai and Bi and the log of the gamma function.
> Should be with 15 digit precision.
>
> I have all of the Cody FORTRAN routines from Netlib
> (http://www.netlib.org/specfun/) compiled to a windows DLL too. All in
> double precision for real inputs only.
> This covers the Bessel functions J0, J1, Y0, Y1, K0, K1, I1 and I2, plus
> J, Y, K and I of integer or fractional order up to order 10 and the psi
> function, gamma function, the log thereof, Dawson's integral and the Airy
> functions Ai and Bi as well as the spherical Bessel functions jn and yn up
> to order 10.
>
> I would need some handholding to contribute this to gnuplot. Perhaps you
> can get me started with that? I would not even know where to begin to find
> the proper documentation.

I have put a patchset on SourceForge
https://sourceforge.net/tracker/?func=detail&aid=3557474&group_id=2055&atid=302055

It detects and links to an external library for Bessel functions and other
special functions including the Airy function. I used libgsl
(the gnu scientific library) for the demo rather than netlib because
I couldn't find a pre-packaged version of netlib, but it should be obvious
how to modify the wrapper functions to match the entry points for a
different library.

After applying the patch and rebuilding gnuplot, you can compare the
built-in and library functions. The relative discrepancy in values
produced by the Airy function is up to 2%, which as you point out is
pretty bad. The maximum discrepancy for other functions I tested
(Bessel, LambertW) ranged from 10^(-14) to 10^(-17).

One advantage to using the netlib routines would be that they handle complex
values, correct? Neither gnuplot nor libgsl currently provide that.

Ethan