为什么scipy.special.hankel1(0, 10**10)返回“nan”？

Question

>>> 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)

我认为这些例子可以在Mathematica或MATLAB中运行。但是，我无法从scipy获得正确的结果。不知道为什么，可能我犯了一些错误？

- fronthem

---

在计算函数时，您可能遇到了浮点限制。该代码是FORTRAN算法，在此处列出（至少我认为这是使用的代码），其中确实提到了一些限制（尽管没有一个符合您的情况）。也许应该向scipy提交问题报告。请注意，仍然可以通过以下方式进行计算:np.complex(jv(0, 1e10), yv(0, 1e10))。 - user707650

@Evert，那-3621592938.02呢？使用你的方法会返回(6.210109276541645e-06+nanj)。 - fronthem

好的，我并不是说我的建议是总体解决方案；它只是在1e10情况下有效的一种权宜之计。也许yv在使用你的新数字时会遇到与hankel1类似的问题。我建议你向scipy团队反映这个问题。Mathematica和matlab可能正在运行比scipy使用的更近期的代码版本。 - user707650

作为一种解决方法，您可以使用Sympy。尝试“import sympy”和“sympy.hankel1(0,10**10).evalf()”。 - Dietrich

Answer 1

我不认为这是一个bug，而是一项功能。

你应该通过执行来启用错误报告功能

scipy.special.errprint(1)

然后你会注意到：

In : scipy.special.hankel1(0, 10**4)
Out: (-0.0070961603533888007+0.0036478055589866053j)

In : scipy.special.hankel1(0, 10**9)
/usr/bin/ipython:1: SpecialFunctionWarning: scipy.special/hankel1:: loss of precision
  #!/usr/bin/python3
Out: (2.4687471886269192e-05-5.2104226538976127e-06j)

In : scipy.special.hankel1(0, 10**10)
/usr/bin/ipython:1: SpecialFunctionWarning: scipy.special/hankel1:: no result obtained
  #!/usr/bin/python3
Out: (nan+nan*j)

- Stefan