Các hàm Bessel trong Python hoạt động với số mũ lớn

Question

Tôi có một số mã sử dụng các hàm Bessel đã sửa đổi của cả thứ tự thứ nhất và thứ 2 (iv và kv). Dường như họ dường như có giới hạn, đó là iv (0,713) và kv (0,697), thêm một cho mỗi và bạn nhận được vô cùng và 0 tương ứng. Đây là một vấn đề đối với tôi bởi vì tôi cần phải sử dụng các giá trị cao hơn điều này, thường lên đến 2000 hoặc hơn. Khi tôi cố gắng phân chia chúng, tôi sẽ lặn bằng 0 hoặc vô cùng có nghĩa là tôi bị lỗi hoặc số không, không phải thứ tôi muốn.

Tôi đang sử dụng scipy bessel functions, có chức năng nào tốt hơn có thể đối phó với nhiều số nhỏ hơn và lớn hơn nhiều hay cách sửa đổi Python để làm việc với những số lớn này. Tôi không chắc vấn đề thực sự ở đây là tại sao Python không thể làm việc này vượt quá 700, có phải là chức năng hay là Python?

Tôi không biết liệu Python có đang thực hiện hay không nhưng tôi chỉ cần 5-10 chữ số đầu tiên * 10^x ví dụ; đó là để nói rằng tôi sẽ không cần tất cả 1000 chữ số, có lẽ đây là vấn đề với cách Python đang làm việc nó ra so với cách Wolfram Alpha đang làm việc nó ra?

Answer 1

Các chức năng iv và kv trong Scipy ít hoặc nhiều như bạn có thể nhận được nếu sử dụng dấu chấm động máy chính xác kép. Như đã lưu ý trong các chú thích ở trên, bạn đang làm việc trong phạm vi mà các kết quả tràn từ phạm vi dấu chấm động.

Bạn có thể sử dụng thư viện mpmath, điều chỉnh độ nổi chính xác (phần mềm) có thể điều chỉnh được để thực hiện việc này. (Nó tương tự như MPFR, nhưng trong 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

mpmath là một thư viện tuyệt vời và là con đường để đi cho các tính toán chính xác cao. Cần lưu ý rằng các hàm này có thể được tính toán từ các thành phần cơ bản của chúng. Vì vậy, bạn không bị buộc phải tuân theo hạn chế của sciper và bạn có thể sử dụng một thư viện có độ chính xác cao khác. Chẳng hạn tối thiểu dưới đây:

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] 

Điều này cho phép:

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] 

này sẽ thất bại cho các giá trị lớn bạn đang tìm kiếm trừ khi dữ liệu cơ bản có độ chính xác cao.

Answer 3

Có thể là vấn đề với chức năng. Đối với x dương lớn, có kv tiệm cận (nu, x) ~ e^{- x}/\ sqrt {x} đối với bất kỳ nu nào. Vì vậy, đối với x lớn bạn kết thúc với các giá trị rất nhỏ. Nếu bạn có thể làm việc với nhật ký hàm Bessel, các vấn đề sẽ biến mất. Scilab khai thác điều này tiệm cận: nó có một băng tham số mặc định là 0, nhưng khi đặt thành 1 sẽ trả về điểm kinh nghiệm (x) * kv (nu, x), và điều này giữ mọi thứ có kích thước hợp lý.

Thực tế, cũng có sẵn trong scipy - scipy.special.kve