Funkcje Bessela w Pythonie, które działają z dużymi wykładnikami

Question

Mam trochę kodu, który używa zmodyfikowanych funkcji Bessela zarówno pierwszego, jak i drugiego rzędu (iv i kv). Irytujące wydają się mieć ograniczenia, są to iv (0,713) i kv (0,697), dodają po jednym, a ty dostajesz odpowiednio nieskończoność i 0. Jest to dla mnie problem, ponieważ muszę użyć wartości wyższych niż ta, często do 2000 lub więcej. Kiedy próbuję przez nie podzielić, kończę nurkowanie 0 lub nieskończoność, co oznacza, że ​​albo otrzymuję błędy, albo zero, z których żadne nie chcę.

Używam scipy bessel functions, czy są jakieś lepsze funkcje, które radzą sobie z dużo mniejszymi i znacznie większymi liczbami lub sposób modyfikowania Pythona do pracy z tymi dużymi liczbami. Nie jestem pewien, jaki jest prawdziwy problem, dlaczego Python nie może ich przetworzyć poza 700, czy jest to funkcja, czy to Python?

Nie wiem, czy Python już to robi, ale potrzebuję tylko pierwszych 5-10 cyfr * 10^x na przykład; to znaczy, że nie potrzebowałbym wszystkich 1000 cyfr, być może jest to problem z tym, jak działa Python w porównaniu do tego, jak pracuje Wolfram Alpha?

Answer 1

Funkcje w Scipy są mniej więcej tak dobre, jak można uzyskać przy użyciu zmiennoprzecinkowej maszyny podwójnej precyzji. Jak zauważono w powyższych komentarzach, pracujesz w zakresie, w którym wynik przekracza zakres zmiennoprzecinkowy.

Można użyć biblioteki mpmath, która wykonuje zmiennoprzecinkową precyzję (oprogramowanie), aby obejść ten problem. (Jest to podobne do MPFR, ale w Pythonie):

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

mpmath jest fantastyczny biblioteki i jest droga dla obliczeń o wysokiej precyzji. Warto zauważyć, że funkcje te można obliczyć na podstawie ich bardziej podstawowych składników. W związku z tym nie jesteś zmuszony przestrzegać ograniczeń scipy i możesz użyć innej biblioteki o wysokiej precyzji. Minimalny przykład poniżej:

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] 

Daje:

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] 

to zawiedzie dla dużych wartościach szukasz, chyba że dane źródłowe ma wysoką precyzję.

Answer 3

Może to być problem związany z funkcją. Dla dużych dodatnich x istnieje asymptotyczny kv (nu, x) ~ e^{- x}/\ sqrt {x} dla dowolnego nu. Tak więc dla dużych x otrzymujesz bardzo małe wartości. Jeśli jesteś w stanie pracować z logiem funkcji Bessela, problemy znikną. Scilab wykorzystuje ten asymptotyczny: ma parametr ice, który domyślnie przyjmuje wartość 0, ale po ustawieniu na 1 zwróci exp (x) * kv (nu, x), a to zachowuje wszystko w rozsądnym rozmiarze.

Właściwie, to samo jest dostępny w scipy - scipy.special.kve