10

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?

4
  • 2
    I don't think its a python problem so much as a double floating point range issue. scipy is providing a wrapper around C code which is actually implementing the bessel function. As such, its limited to the range a double can accommodate.
    – sizzzzlerz
    Dec 5, 2012 at 15:31
  • yes, I just did sys.float_info and guess what? max_10_exp=308 which is more or less exactly the answer to the bessel functions at the limit. This is pretty bad news for me. How is Wolfram Alpha able to work it out though?
    – Rapid
    Dec 5, 2012 at 15:34
  • 1
    Magic? I don't know but I'm pretty sure Alpha gets its code base from Wolfram's Mathematica which is a pretty sophisticated tool. They've implemented some sort of algorithm that allows them to return essentially unlimited precision and range for transcendental functions like the bessel.
    – sizzzzlerz
    Dec 5, 2012 at 16:09
  • @sizzzzlerz: It's actually written in Fortran :) Wolfram Alpha is probably using arbitrary-precision math, like the mpmath library for Python.
    – endolith
    Sep 16, 2015 at 0:00

4 Answers 4

8

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')
3
  • This seems to get the correct answers so thank you very much. The only issue I've got is trying to incorporate the mpf stuff into the rest of my code. As it's not a float it's giving me a lot of errors and mpf doesn't seem to accept arrays? You've answered my question though so I'm just going to have to do a fair amount of research into mpf and mpfr etc. Is there no way to have 1.2597e-746 stroed as a float, I mean there is only 11 things to store. Hmmm.
    – Rapid
    Dec 5, 2012 at 17:12
  • 1
    You need to keep everything as mpf. However, a good idea can also be to work with the logarithms of small numbers rather than the logarithms themselves. If you need to add numbers, use logaddexp. Or, you can rethink the problem and do some mathematics to try to get rid of the huge and tiny numbers.
    – pv.
    Dec 5, 2012 at 21:06
  • @Rapid: You could write your own class to store the decimal points and the exponent separately if you wanted, but is that really worth the effort?
    – endolith
    Sep 16, 2015 at 0:02
4

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.

3
  • That's exactly what I was looking for :) Dec 10, 2019 at 10:52
  • I might need to ask this in another function but—any tips of what happens if kve is overflowing? Scipy's kve(480.49, 26.44) overflows because the order (first argument) is so high 😕. Is my only choice to use Mpmath (which handles this fine)? Feb 7, 2023 at 2:46
  • For large-order kve, I managed to use one of the asymptotic forms on dlmf.nist.gov/10.41 Feb 7, 2023 at 3:36
3

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.

2
  • If those will also fail for large values, then how is this an answer? Are those more precise than the dedicated functions?
    – endolith
    Sep 16, 2015 at 0:03
  • @endolith sorry if that last line isn't clear. What I'm implying is that large values needed a correspondingly large precision (set in mpmath) to be accurate. The advantage of using mpmath comes from being able to set this value.
    – Hooked
    Sep 16, 2015 at 13:58
1

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

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.