Python中使用大指数的贝塞尔函数

Question

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

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

我不知道Python是否已经在做了，但我只需要前5-10个数字*10^x；也就是说，我不需要全部1000个数字，也许这就是Python与WolframAlpha相比的问题所在？

Answer 1

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')

Answer 2

您可以使用指数缩放的修改贝塞尔函数直接执行此操作，该函数不会溢出。它们被实现为CCD_ 4和CCD_。例如，第一类修改后的贝塞尔函数special.iv(0, 1714)将溢出。然而，只要你不使用已经溢出的东西的对数，它的对数就会定义得很好：

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

其他容易溢出的函数也可用作日志或缩放版本。

Answer 3

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 4

可能是函数出了问题。对于大的正x，任何nu都有渐近的kv(nu，x)~e^｛-x｝/\sqrt｛x｝。所以对于大的x，你会得到非常小的值。如果你能够使用贝塞尔函数的对数，问题就会消失。Scilab利用了这种渐近性：它有一个默认为0的参数ice，但当设置为1时，它将返回exp(x)*kv(nu，x)，这使所有东西都保持合理的大小。

事实上，scipy-scipy.special.kve 也有同样的功能