slatec bessel_k(n, %i*x): lack of precision for large x



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