search

scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
12.95k stars 5.15k forks source link

scipy.special.hankel1(0, 10**10) return `nan`? #3737

Open ghost opened 10 years ago

ghost commented 10 years ago 
>>> import scipy.special
>>> scipy.special.hankel1(0, 10**8)
(3.2060295340412074e-05+7.3063911655217044e-05j)
>>> scipy.special.hankel1(0, 10**9)
(2.4687471886269185e-05-5.2104226538976152e-06j)
>>> scipy.special.hankel1(0, 10**10)
(nan+nan*j)

I'm not sure but i think this should be done by either Mathematica or MATLAB but i got nan when use scipy.

WarrenWeckesser commented 10 years ago

For anyone looking into this: mpmath has an implementation of hankel1. For example,

In [4]: from sympy import mpmath

In [5]: mpmath.mp.dps = 40

In [6]: mpmath.hankel1(0, 1e10)
Out[6]: mpc(real='0.000002175591750246891726859055283638209197067145', imag='-0.000007676508175792936690488745267426925572737321')
pv commented 10 years ago

The function is rapidly oscillating at large arguments, and AMOS refuses to give results --- presumably because it is not possible to guarantee relative accuracy in floating point.

The real-valued implementations of jv, yv do continue giving correct results (as far as possible in floating point) also for larger values, so fixing the issue on real axis is fairly easy.

In general (also in Matlab and whatever environment) values of oscillating functions at large floating point arguments lose relative precision, but it may still be better to return best-efforts results rather than nans. This is what libc does with sin/cos.