Python中使用大指数的Bessel函数

Question

我有一些代码使用了修改后的贝塞尔函数的一阶和二阶（iv和kv）。令人恼火的是，它们似乎有极限，它们是iv（0713）和kv（0697），每加一个，分别得到无穷和0。这对我来说是个问题，因为我需要使用高于此值的值，通常高达2000或更多。当我试图除以这些值时，我最终得到了0或无穷大，这意味着我要么得到了错误，要么得到了0，这都不是我想要的。

我使用的是scipy bessel functions，有没有更好的函数可以处理更小、更大的数字，或者是修改Python以处理这些大数字的方法。我不确定这里真正的问题是什么，为什么Python不能在700之后解决这些问题，是函数还是Python？

我不知道Python是否已经这样做了，但我只需要前5-10个数字*10^x；也就是说，我不需要所有1000个数字，也许这就是Python是如何工作的问题，而Wolfram Alpha是如何工作的？

Answer 1

mpmath是一个非常棒的库，是进行高精度计算的途径。值得注意的是，这些函数可以从它们更基本的组成部分计算出来。因此，您不必遵守scipy的限制，您可以使用不同的高精度库。下面是最简单的例子：

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]

这就提供了：

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]

除非底层数据具有高精度，否则对于要查找的大值，这将失败。

Answer 2

可能是函数有问题。对于大正x，任何nu都有渐近的kv（nu，x）~e^{-x}/\sqrt{x}。所以对于大x，最终得到的是非常小的值。如果您能够使用贝塞尔函数的日志，问题将消失。Scilab利用了这种渐近性：它的参数ice默认为0，但当设置为1时，将返回exp（x）*kv（nu，x），这将保持所有内容的大小合理。

实际上，在scipy-scipy.special.kve中也有相同的功能

Answer 3

Scipy中的iv和kv函数或多或少都是使用双精度机器浮点时可以得到的。如上所述，您正在结果从浮点范围溢出的范围内工作。

您可以使用mpmath库（它执行可调精度（软件）浮点）来解决这个问题。（类似于MPFR，但在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')