CODEWITHSUNDEEP
pythonscipybessel-functions
Rapid
Rapid

Reputation: 1502

Bessel functions in Python that work with large exponents

I've got some code that uses the modified Bessel functions of both 1st and 2nd order (iv and kv). Annoyingly they seem to have limits, those are iv(0,713) and kv(0,697), add one to each and you get infinity and 0 respectively. This is a problem for me because I need to use values higher than this, often up to 2000 or more. When I try to divide by these I end up dividing by 0 or by infinity which means I either get errors or zeros, neither of which I want.

I'm using the scipy bessel functions, are there any better functions that can cope with much smaller and much larger numbers, or a way of modifying Python to work with these big numbers. I'm unsure what the real issue here is as to why Python can't work these out much beyond 700, is it the function or is it Python?

I don't know if Python is already doing it but I'd only need the first 5-10 digits *10^x for example; that is to say I wouldn't need all 1000 digits, perhaps this is the problem with how Python is working it out compared to how Wolfram Alpha is working it out?

Upvotes: 10

Views: 11309

Answers (4)

Timothy Brandt
Timothy Brandt

Reputation: 41

You can do this straightforwardly using the exponentially scaled modified Bessel functions, which will not overflow. These are implemented as special.ive and special.kve. For example, the modified Bessel function of the first kind, special.iv(0, 1714), will overflow. However, its logarithm will be perfectly well-defined, as long as you aren't taking the log of something that has already overflowed:

In [1]: import numpy as np

In [2]: from scipy import special

In [3]: np.log(special.iv(0, 1714))
Out[3]: inf

In [4]: np.log(special.kv(0, 1714))
Out[4]: -inf

In [5]: np.log(special.ive(0, 1714)) + 1714
Out[5]: 1709.3578418673253

In [6]: np.log(special.kve(0, 1714)) - 1714
Out[6]: -1717.4975741044941

Other functions that readily overflow are also available as logs or in scaled versions.

Upvotes: 4

Hooked
Hooked

Reputation: 88208

mpmath is a fantastic library and is the way to go for high-precision calculations. It is worth noting that these functions can be computed from their more basic constituents. Thus, you are not forced to abide by scipy's restriction and you can use a different high precision library. Minimal example below:

import numpy as np
from scipy.special import *

X = np.random.random(3)

v = 2.000000000

print "Bessel Function J"
print jn(v,X)

print "Modified Bessel Function, Iv"
print ((1j**(-v))*jv(v,1j*X)).real
print iv(v,X)   

print "Modified Bessel Function of the second kind, Kv"
print (iv(-v,X)-iv(v,X)) * (np.pi/(2*sin(v*pi)))
print kv(v,X)

print "Modified spherical Bessel Function, in"
print np.sqrt(np.pi/(2*X))*iv(v+0.5,X)
print [sph_in(floor(v),x)[0][-1] for x in X]   

print "Modified spherical Bessel Function, kn"
print np.sqrt(np.pi/(2*X))*kv(v+0.5,X)
print [sph_kn(floor(v),x)[0][-1] for x in X]

print "Modified spherical Bessel Function, in"
print np.sqrt(np.pi/(2*X))*iv(v+0.5,X)
print [sph_in(floor(v),x)[0][-1] for x in X]

This gives:

Bessel Function J
[ 0.01887098  0.00184202  0.08399226]

Modified Bessel Function, Iv
[ 0.01935808  0.00184656  0.09459852]
[ 0.01935808  0.00184656  0.09459852]

Modified Bessel Function of the second kind, Kv
[  12.61494864  135.05883902    2.40495388]
[  12.61494865  135.05883903    2.40495388]

Modified spherical Bessel Function, in
[ 0.0103056   0.00098466  0.05003335]
[0.010305631072943869, 0.00098466280846548084, 0.050033450286650107]

Modified spherical Bessel Function, kn
[   76.86738631  2622.98228411     6.99803515]
[76.867205587011171, 2622.9730878542782, 6.998023749439338]

Modified spherical Bessel Function, in
[ 0.0103056   0.00098466  0.05003335]
[0.010305631072943869, 0.00098466280846548084, 0.050033450286650107]

This will fail for the large values you are looking for unless the underlying data has high precision.

Upvotes: 3

ChrisR
ChrisR

Reputation: 11

Could be the problem is with the function. For large positive x, there is the asymptotic kv(nu,x) ~ e^{-x}/\sqrt{x} for any nu. So for large x you end up with very small values. If you are able to work with the log of the Bessel function instead, the problems will vanish. Scilab exploits this asymptotic: its has a parameter ice which defaults to 0, but when set to 1 will return exp(x)*kv(nu,x), and this keeps everything of reasonable size.

Actually, the same is available in scipy - scipy.special.kve

Upvotes: 1

pv.
pv.

Reputation: 35135

The iv and kv functions in Scipy are more or less as good as you can get if using double precision machine floating point. As noted in the comments above, you are working in the range where the results overflow from the floating point range.

You can use the mpmath library, which does adjustable precision (software) floating point, to get around this. (It's similar to MPFR, but in Python):

In [1]: import mpmath

In [2]: mpmath.besseli(0, 1714)
mpf('2.3156788070459683e+742')

In [3]: mpmath.besselk(0, 1714)
mpf('1.2597398974570405e-746')

Upvotes: 8

Related Questions