9

我有一些代码使用修改后的一阶和二阶（iv 和 kv）贝塞尔函数。令人讨厌的是，它们似乎有限制，它们是 iv(0,713) 和 kv(0,697)，每加一个，你分别得到无穷大和 0。这对我来说是个问题，因为我需要使用高于此的值，通常高达 2000 或更多。当我尝试除以这些时，我最终会除以 0 或无穷大，这意味着我要么得到错误，要么得到零，这两者都不是我想要的。

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

我不知道 Python 是否已经在这样做，但我只需要前 5-10 个数字 *10^x 例如；也就是说，我不需要所有 1000 位数字，与 Wolfram Alpha 的工作方式相比，这可能是 Python 工作方式的问题？

4

4 回答 4

8

如果使用双精度机器浮点，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')
于 2012-12-05T16:56:12.603 回答
4

您可以使用指数缩放的修正贝塞尔函数直接执行此操作，该函数不会溢出。这些实现为special.ive和special.kve。例如，第一类修正贝塞尔函数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

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

于 2018-06-15T21:34:21.163 回答
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]

除非基础数据具有高精度，否则您正在寻找的大值将失败。

于 2012-12-05T17:31:53.170 回答
1

可能问题出在功能上。对于大的正 x，任何 nu 都有渐近的 kv(nu,x) ~ e^{-x}/\sqrt{x}。因此，对于大 x，您最终会得到非常小的值。如果您能够改为使用 Bessel 函数的日志，问题就会消失。Scilab 利用了这种渐近性：它有一个参数 ice，默认为 0，但当设置为 1 时，将返回 exp(x)*kv(nu,x)，这样可以保持一切合理的大小。

实际上，在 scipy - scipy.special.kve 中也可以使用

于 2015-01-20T09:58:20.753 回答