Maybe you are already aware of COMPUTATION OF BESSEL AND AIRY FUNCTIONS AND OF RELATED GAUSSIAN QUADRATURE FORMULAE by WALTER GAUTSCHI https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/169.pdf
The code is in Fortran though. https://www.cs.purdue.edu/archives/2001/wxg/codes/ IMO it would be easier, however, if we took it from GSL: https://www.gnu.org/software/gsl/doc/html/specfunc.html Bill wrote an interface: https://github.com/billpage/gsla/blob/master/gsl.lisp On 08.04.2021 18:13, Waldek Hebisch wrote: > This is mainly to let you know why bug fixes take so much > time. There is old bug, namely our numeric Airy functions > fails for positive real arguments (and gives quite wrong > results for some complex arguments). Problem is that > our current formula for Airy Ai uses Bessel functions and > square root. For some arguments (including positive reals) > we get wrong branch of square root and conseqently > completely wrong result. Mathematically, fix > looks easy. Namely there is formula for Airy Ai: > > Ai(x) = (3)^(-2/3)*(1/Gamma(2/3))*H01(2/3, (1/9)*x^3) - > (3)^(-1/3)*(1/Gamma(1/3)* x*H01(4/3, (1/9)*x^3) > > where H01 is hypergeometric function. H01 is analytic > for all complex x, so no problem with branches. But > it has two problem, one fundamental. First, for > positve reals H01 grows quite fast while Ai decays fast. > So we are subracting two large numbers to get tiny > difference, which leads to significant loss of > accuracy (this is fundamental problem). This one can > be corrected by expressing Ai in terms of Bessel K > function (this only for limited set of arguments, as > outherwise we would get back problems with branches). > However, our Bessel K uses formula that is similar > to expression for Ai above, and in particular has > the sane problem with loss of accuracy that we are > trying to avoid... > > Other problem is that for some arguments our H01 is > quite inaccurate. Now, H01 is closely related to Bessel > J function: > > H01(p, x) = (sqrt(-x))^(-(p-1))*Gamma(p)J_{p-1}(2 sqrt(-x)) > > For efficiency it makes some sense to have separate > special method for Bessel J and for H01, but one > can use common method just doing little adjustment > implied by equation above. Certainly, there is no > reson to have significant difference in accuracy. > For positive reals formula for H01 means that we take > square root of negative number, that is evaluate > Bessel J at imaginary arguments. Alternativly, one > can use Bessel I function. Again, for complex arguments > Bessel J and Bessel I are closely related and one > can use common method. But our (Boot) code for > Bessel J and Bessel K differ quite a lot. In particular > Bessel I raises error for most arguments. This means > we have not only problem with Airy functions, but > also we have problem with Bessel functions (and H01). > > I looked at methods used to compute Bessel functions. > There are some: > - direct use of power series (good only for small x) > - asymptotic expansions (needs large x or p) > - Chebyshev series (moderate arguments) > - recurences > > There are several asymptotic expansions, each with its > validity conditions and with different errors. It seems > that most expansions are given without error bounds, so > it is not clear where they are really useful. Recurences > have problem with stability, some are claimed stable, > but this is probably valid only for real p and real x. > My calculations suggest that in "interesting" ranges > of p and x recurences are unstable. What looked > like easy bugfix expanded to a research program... > > Anyway, it took some time and I probably will postopone > fixing this one past current release. But there are > few other "easy looking" potential fixes... > > BTW: Obvious thing to do is to look at other implementations. > ATM I have not found "full range" implementation, that > is one giving good accuracy for all arguments where > function may be computed (there good looking implementations > for real p). According to various papers Mathematica > is using multiple procision which gives good accuracy, > but is so slow for some arguments that range is > limited by compute time. AFAICS the same it true of > mpmath. > -- You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/07d9bb7d-d412-0d06-a619-7f46cd65a2f9%40gmail.com.
