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.
--
Waldek Hebisch
--
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/20210408161321.GA18690%40math.uni.wroc.pl.
