Path: csiph.com!usenet.pasdenom.info!weretis.net!feeder4.news.weretis.net!eternal-september.org!feeder.eternal-september.org!mx04.eternal-september.org!.POSTED!not-for-mail From: sfeam Newsgroups: comp.graphics.apps.gnuplot Subject: Re: Airy function wrong Followup-To: comp.graphics.apps.gnuplot Date: Tue, 14 Aug 2012 11:29:06 -0700 Organization: gnuplot development team Lines: 122 Message-ID: References: <5023c85e$0$6946$e4fe514c@news2.news.xs4all.nl> <5026a763$0$6930$e4fe514c@news2.news.xs4all.nl> Reply-To: sfeam@users.sourceforge.net Mime-Version: 1.0 Content-Type: text/plain; charset="ISO-8859-1" Content-Transfer-Encoding: 7Bit Injection-Date: Tue, 14 Aug 2012 18:29:09 +0000 (UTC) Injection-Info: mx04.eternal-september.org; posting-host="8e86a57dfa599721f116da4577f3d1af"; logging-data="11702"; mail-complaints-to="abuse@eternal-september.org"; posting-account="U2FsdGVkX1+jspHB2AJbhSeyzZ0FKB31" User-Agent: KNode/4.4.9 Cancel-Lock: sha1:520MXXN6vLXaD7F+vQM6uoBC/3c= Xref: csiph.com comp.graphics.apps.gnuplot:1324 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" 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