mpmath icon indicating copy to clipboard operation
mpmath copied to clipboard

Negative integer order of the BesselI function

Open tk-yoshimura opened this issue 1 year ago • 3 comments

mpmath version: 1.3.0

I am using mpmath to generate values for verification of my own complex Bessel function calculation routine. https://github.com/tk-yoshimura/ComplexBessel

I have discovered that mpmath cannot compute for large negative integers in Bessel I functions, even though it can use the inversion formula. https://dlmf.nist.gov/10.27

Example:

import mpmath as mpf

mpf.mp.dps = 64

x = 2

nus = [31, 32, 63, 64, 127, 128,
      -31, -31.25, -32, -32.25, 
      -63, -63.25, -64, -64.25, 
      -127, -127.25, -128, -128.25
]

for nu in nus:
    print(f"nu = {nu}, x = {x}")
    
    try:
        y = mpf.besseli(nu, x)
        print(f"{mpf.nstr(y, 40, strip_zeros=False, min_fixed=2, max_fixed=1)}")
    except ValueError as e:
        print(e)

    print("\n")

Error Message:

nu = -127, x = 2

hypercomb() failed to converge to the requested 216 bits of accuracy
using a working precision of 4814 bits. The function value may be zero or
infinite; try passing zeroprec=N or infprec=M to bound finite values between
2^(-N) and 2^M. Otherwise try a higher maxprec or maxterms.

tk-yoshimura avatar Oct 14 '24 14:10 tk-yoshimura

I have discovered that mpmath cannot compute for large negative integers in Bessel I functions

@tk-yoshimura, have you tried suggestions?

E.g.

>>> besseli(-127, 2, maxprec=8000)
mpf('3.345358761443415013354345973251886375421555647081543375756063117036e-214')
>>> nstr(_, n=20)
'3.3453587614434150134e-214'

Mathematica:

In[4]:= N[BesselI[-127, 2], 20]

                                -214
Out[4]= 3.3453587614434150134 10

skirpichev avatar Oct 14 '24 16:10 skirpichev

Thank you for your comment. I will use that code to generate the validation values.

It seems odd that at the same precision setting, the same nu, x would yield Bessel J and Y values, and nu=-127.25 would yield Bessel I values, but only nu=-127 would not. Will there be a correction by inversion formula in the future?

tk-yoshimura avatar Oct 15 '24 04:10 tk-yoshimura

Will there be a correction by inversion formula in the future?

This does make sense at least as an enhancement request.

skirpichev avatar Oct 15 '24 04:10 skirpichev

PR is ready: https://github.com/mpmath/mpmath/pull/889

skirpichev avatar Dec 05 '24 05:12 skirpichev