Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.graphics.apps.gnuplot > #1324

Re: Airy function wrong

Followups directed to: comp.graphics.apps.gnuplot

Show all headers | View raw

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



> 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" <sfeam@users.sourceforge.net> wrote in message
> news:k01e9m$sr4$1@dont-email.me...
>> Alex van der Spek wrote:
>>
>>> 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
>>
>> 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

Back to comp.graphics.apps.gnuplot | Previous | Next — Previous in thread | Find similar

Thread

Airy function wrong "Alex van der Spek" <zdoor@xs4all.nl> - 2012-08-09 16:25 +0200
  Re: Airy function wrong sfeam <sfeam@users.sourceforge.net> - 2012-08-09 15:31 -0700
    Re: Airy function wrong "Alex van der Spek" <zdoor@xs4all.nl> - 2012-08-11 20:41 +0200
      Re: Airy function wrong sfeam <sfeam@users.sourceforge.net> - 2012-08-14 11:29 -0700

csiph-web