2012-12-05 5 views
7

1 차 및 2 차 (iv 및 kv)의 수정 된 Bessel 함수를 사용하는 코드가 있습니다. Annoyingly 그들에는 한계가 있고, 그것들은 iv (0,713)과 kv (0,697)이고, 각각 하나를 더하면 무한과 0을 각각 얻습니다. 이 문제는 저보다 더 높은 값 (보통 2000 이상)을 사용해야하기 때문에 문제가됩니다. 나는 이것들에 의해 나눌 때 나는 0이나 무한대로 다이빙을 끝내고, 이는 내가 에러 또는 0을 얻는다는 것을 의미한다.큰 지수로 작동하는 Python의 Bessel 함수

나는 더 작고 훨씬 더 큰 숫자에 대처할 수있는 더 나은 함수 나이 큰 숫자로 작업 할 파이썬을 수정하는 방법이 있습니까? scipy bessel functions을 사용하고 있습니다. 파이썬이 700을 훨씬 뛰어 넘는이 기능을 수행 할 수없는 이유에 관해서 진짜 문제가 무엇인지 확신 할 수 없습니다. 파이썬은 기능입니까 아니면 파이썬입니까?

파이썬이 이미 그것을하고 있는지는 모르지만 예를 들어 처음 5-10 자리 * 10^x 만 필요합니다. 말하자면, 1000 자리 숫자가 모두 필요하지는 않습니다. 아마도 이것이 Wolfram Alpha가 어떻게 작동하는지에 비해 파이썬이 어떻게 작동하는지에 대한 문제일까요?

+2

을 scipy.special.kve. scipy는 실제로 베셀 함수를 구현하는 C 코드 주위에 래퍼를 제공합니다. 따라서 두 배의 범위로 제한됩니다. – sizzzzlerz

+0

예, 방금 sys.float_info를 수행하고 무엇을 추측합니까? max_10_exp = 308 이는 한계에서 베셀 함수에 대한 해답입니다. 이것은 나에게 나쁜 소식이다. Wolfram Alpha는 어떻게 그것을 해결할 수 있습니까? – Rapid

+1

Magic? 몰라요.하지만 꽤 정교한 툴인 Wolfram의 Mathematica에서 Alpha가 코드 기반을 얻었습니다. 그들은 베셀과 같은 초월 함수에 대해 본질적으로 무제한의 정밀도와 범위를 반환 할 수있는 일종의 알고리즘을 구현했습니다. – sizzzzlerz

답변

7

Scipy의 ivkv 함수는 배정도 부동 소수점을 사용하는 경우 얻을 수있는만큼 훌륭합니다. 위의 주석에서 언급했듯이 결과가 부동 소수점 범위에서 오버플로되는 범위에서 작업하고 있습니다.

조정 가능한 정밀도 (소프트웨어) 부동 소수점을 사용하는 mpmath 라이브러리를 사용하여이 문제를 해결할 수 있습니다. (그것은하지만, 파이썬, MPFR 비슷) :

In [1]: import mpmath 

In [2]: mpmath.besseli(0, 1714) 
mpf('2.3156788070459683e+742') 

In [3]: mpmath.besselk(0, 1714) 
mpf('1.2597398974570405e-746') 
+0

이것은 정답을 얻은 것 같아서 대단히 감사합니다. 내가 가진 유일한 문제는 내 코드의 나머지 부분에 mpf 파일을 통합하려고하는 것입니다. float가 아니기 때문에 많은 오류가 발생하고 mpf는 배열을 허용하지 않습니다. 당신은 내 질문에 대답했지만 mpf와 mpfr 등등에 대한 상당한 연구를해야 할 것입니다. 1.2597e-746을 플로트로 쓰는 방법이 없으며 단지 11 가지 밖에 없다는 것을 의미합니다. 저장할 수 있습니다. 흠. – Rapid

+1

모든 것을'mpf'로 유지해야합니다. 그러나 좋은 생각은 로그 자체가 아닌 작은 수의 로그로 작업하는 것입니다. 숫자를 추가하려면'logaddexp'를 사용하십시오. 또는 문제를 다시 생각하고 거대한 숫자와 작은 숫자를 없애기 위해 수학을 할 수 있습니다. –

+0

@Rapid : 원한다면 소수점과 지수를 따로 저장하는 클래스를 직접 작성할 수는 있지만 그만한 가치가 있습니까? – endolith

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] 

이렇게하면 기본 데이터는 높은 정밀도를 가지고 있지 않는 찾고있는 큰 값을 실패합니다.

+0

큰 값에 대해서도 실패한다면, 어떻게 대답할까요? 그것들은 헌신적 인 기능보다 정확합니까? – endolith

+0

@endolith 마지막 줄이 분명하지 않으면 죄송합니다. 필자가 시사하는 바는 큰 값은 정확하기 위해서 상응하는 큰 정밀도 ('mpmath'에서 설정)가 필요하다는 것입니다. mpmath를 사용하면이 값을 설정할 수 있다는 장점이 있습니다. – Hooked

1

문제가있을 수 있습니다. 큰 양의 x에 대해 임의의 뉴에 대한 점근선 kv (nu, x) ~ e^{- x}/\ sqrt {x}가 있습니다. 따라서 큰 x의 경우 매우 작은 값으로 끝납니다. 대신 Bessel 함수의 로그로 작업 할 수 있다면 문제는 사라질 것입니다. Scilab은이 점을 이용합니다 : 그것의 디폴트는 0이지만, 1로 설정하면 exp (x) * kv (nu, x)가 리턴됩니다. 이것은 적당한 크기의 모든 것을 유지합니다.

는 사실, 같은 scipy에서 사용할 수 있습니다 - 나는 그것에게 이중 부동 소수점 범위 문제로 너무 많은 파이썬 문제를 생각하지 않는다