slatec bessel_k(n, %i*x): lack of precision for large x
Subject: slatec bessel_k(n, %i*x): lack of precision for large x
From: Edwin Woollett
Date: Fri, 2 Mar 2012 11:47:05 -0800
A test of the slatec fortran - lisp translated
bessel_k(n,%i*x) for large (real) x shows
either a lack of precision or refusal to compute
for large enough values of x. This may be
the reason quad_qagi has problems with
getting convergence for semi-infinite integrals
involving bessel_k(2,%i*x) which I have been
looking at.
First for the large x test: the ierr=3 means
"lack of precision warning" in evaluation of
the bessel function, while ierr=4 means
" no computation because abs(z) is too
large".
(See the slatec fortran code zbesk.f link on the
page
http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/src/numerical/slatec/fortran/
I load nint.mac so I can use my
functions fchop and cfloat here:
-------------------------------------------------------
(%i1) load(nint);
(%o1) "c:/work2/nint.mac"
(%i2) for n:2 thru 10 do
print(n," ",fchop(cfloat(bessel_k(2,%i*10^n))))$
2 0.120695-0.0338173*%i
3 -0.03892*%i-0.00748561
4 0.0111478*%i-0.00572773
zbesk ierr = 3
zbesk ierr = 3
5 0.00270057*%i-0.00290084
zbesk ierr = 3
zbesk ierr = 3
6 0.00114035-5.20003431E-4*%i
zbesk ierr = 3
zbesk ierr = 3
7 1.36403863E-4*%i-3.72120436E-4
zbesk ierr = 3
zbesk ierr = 3
8 -5.03601919E-5*%i-1.14768525E-4
zbesk ierr = 3
zbesk ierr = 3
9 8.18451269E-6-3.87789902E-5*%i
zbesk ierr = 4
zbesk ierr = 4
10 0.0
-------------------------------------
Already at arg = %i*10^4 bessel_k(2,arg)
has precision problems.
Finally, ar arg = %i*10^10 bessel_k(2,arg)
shuts down with the return of 0.0.
------------------------
How this affects ability of quad_qagi is less clear.
I use nint (which will call quad_qagi here) for a series
of less and less convergent (or easy) integrals.
The real and imaginary parts of the integrand
are separately submitted to quad_qagi by nint.
--------------------------------------------------
/* the denominator x^6 heavily damps the
small wandering values of the numerator,
so quag_qagi doesn't ask for large x values
from bessel_k here: */
(%i3) nint(bessel_k(2,%i*x)/x^6,x,1,inf);
(%o3) 0.0535655*%i-0.401416
(%i4) noutL;
(%o4) [[[qagi,-0.401416,1.06713183E-9,195,0],
[qagi,0.0535655,1.62460657E-10,255,0]]]
(%i5) nint(bessel_k(2,%i*x)/x^5,x,1,inf);
(%o5) 0.0728441*%i-0.474075
(%i6) noutL;
(%o6) [[[qagi,-0.474075,2.49004643E-9,285,0],
[qagi,0.0728441,5.43186958E-10,405,0]]]
/* here zbesk starts to complain about loss of
precision although the final answer agrees with
Mma = wolfram alpha */
(%i7) nint(bessel_k(2,%i*x)/x^4,x,1,inf);
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
(%o7) 0.107823*%i-0.576049
/* this calculation also agrees with Mma, but the zbesk
error messages clutter up the output, and
show why it is desireable to have a way to
surpress the printing of these messages to
the screen! */
(%i8) nint(bessel_k(2,%i*x)/x^3,x,1,inf);
(%i8) nint(bessel_k(2,%i*x)/x^3,x,1,inf);
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
(%o8) 0.17808*%i-0.724001
(%i9) noutL;
(%o9) [[[qagi,-0.724001,1.95810235E-9,3135,0],
[qagi,0.17808,2.86307006E-10,5925,0]]]
running
(%10) nint(bessel_k(2,%i*x)/x^2,x,1,inf);
generates _many_ pages of the same error message
and "too many subintervals done".
.....................
.......................
zbesk ierr = 3
zbesk ierr = 3
zbesk ierr = 3
quadpack error code = too many subintervals done
real part returns false
(%o10) false
note the imaginary part gave correct answer:
(%i11) noutL;
(%o11) [[qagi,-0.939769,3.75759298E-7,23985,1]]
Mma's answer = 0.332602 i - 0.939769
-----------------------------------
The bottom line here is: are there better bessel function
routines available which are reasonably accurate for large
abs(z) args??
Ted Woollett