Bessel函數在Python中與大指數一起工作

我有一些代碼使用第一和第二順序（iv和kv）的修改的貝塞爾函數。令人討厭的是，他們似乎有限制，那些是iv（0,713）和kv（0,697），分別加1，你分別得到無窮大和0。這對我來說是一個問題，因爲我需要使用比這更高的值，通常高達2000或更多。當我試圖用這些方法進行劃分時，我最終會以0或無窮大進行跳水，這意味着我要麼獲得錯誤或零點，這兩者都不是我想要的。Bessel函數在Python中與大指數一起工作

我使用的是scipy bessel functions，是否有更好的函數可以處理更小和更大的數字，或者是一種修改Python來處理這些大數字的方法。我不確定這裏真正的問題是爲什麼Python不能在700以外完成這些工作，是它的功能還是Python？

我不知道Python是否已經在做了，但我只需要前5-10位* 10^x的，例如;也就是說我不需要全部1000位數字，這可能與Python如何解決這個問題相比，Wolfram Alpha是如何解決這個問題的？

來源

2012-12-05 Rapid

我不認爲它是一個Python的問題，這麼多雙浮點範圍問題。 scipy提供了一個實際實現bessel函數的C代碼包裝器。因此，它僅限於雙倍可以容納的範圍。 – sizzzzlerz

是的，我只是做了sys.float_info並猜測是什麼？ max_10_exp = 308，這或多或少正好是限制條件下bessel函數的答案。這對我來說是個壞消息。 Wolfram Alpha如何能夠解決這個問題？ – Rapid

魔法？我不知道，但我很確定Alpha從Wolfram的Mathematica中獲得它的代碼庫，這是一個非常複雜的工具。他們已經實現了一些算法，使他們能夠返回基本上無限的精度和範圍，像bessel這樣的超越函數。 – sizzzzlerz

在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-05 16:56:12

這似乎得到了正確的答案，所以非常感謝你。我遇到的唯一問題是試圖將mpf的東西加入到我的其他代碼中。因爲它不是浮點數，它給了我很多錯誤，mpf似乎不接受數組？你已經回答了我的問題，所以我只需要對mpf和mpfr等進行大量的研究。是否沒有辦法將1.2597e-746作爲浮點運算，我的意思是隻有11件事情儲藏。嗯。 – Rapid

您需要將所有內容保存爲「mpf」。然而，一個好主意也可以是使用小數的對數而不是對數本身。如果您需要添加數字，請使用'logaddexp'。或者，你可以重新思考這個問題，並做一些數學試圖擺脫巨大而微小的數字。 –

@Rapid：如果你想要，你可以編寫自己的類來存儲小數點和指數，但這真的值得嗎？ – endolith

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-05 17:31:53 Hooked

如果這些也會失敗的大值，那麼這是一個答案？這些比專用功能更精確嗎？ – endolith

@ndndolith抱歉，如果最後一行不清楚。我暗示的是，大數值需要一個相對較大的精度（在mpmath中設置）纔是準確的。使用mpmath的優勢來自於能夠設置此值。 – Hooked

可能是問題與功能。對於大的正x，任何nu都有漸近的kv（nu，x）〜e^{ - x}/\ sqrt {x}。所以對於大型的x，你會得到非常小的值。如果您能夠使用貝塞爾函數的日誌，那麼問題就會消失。 Scilab利用這個漸近：它有一個默認爲0的參數ice，但當設置爲1時會返回exp（x）* kv（nu，x），並且這會保持一切合理的大小。

其實，同樣是在SciPy的可用 - scipy.special.kve

來源

2015-01-20 09:58:20 ChrisR

Bessel函數在Python中與大指數一起工作

回答

相關問題