Here's a paper that discusses Bessel function computation in the context of arbitrary precision. Miller's algorithm (backward recurrence) is one method, but there are other methods for large argument.
https://people.eecs.berkeley.edu/~fateman/papers/hermite.pdf Have fun. On Thursday, April 8, 2021 at 10:17:40 AM UTC-7 Kurt Pagani wrote: > 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/215bd762-68e9-464a-8fdf-8109a280645fn%40googlegroups.com.
