python - scipy.special.ellipkinc 中的错误 - 不完全椭圆积分

Question

在使用 SciPy 的 scipy.special.ellipeinc 和 ellipkinc 时，似乎存在一些数值不稳定的孤岛。例如，

>>> from scipy.special import ellipkinc
>>> ellipkinc(0.9272952180016123, 0.68359375000000011)
nan
>>> ellipkinc(0.9272952180016123, 0.6835937500000002)
2.0518660200390668
>>> ellipkinc(0.9272952180016123, 0.68359375)
1.0259330100195332
>>> ellipkinc(0.9272952180016123, 0.68359374)
1.0259330081166262

这发生在 k^2.sin^2(phi) 接近 0.3 的地方，但这里的椭圆积分本身并没有什么异常，所以大概是数值的东西。我不太了解这个算法的内部工作原理，无法说出问题所在，那么我最好的选择是什么？我想到了:圆(0.68359375000000011,8)说，但这肯定会减慢我的代码速度？

Answer 1

(这更多是对问题的扩展评论，而不是对如何处理的回答。)

这看起来像是 ellipkinc 中的错误。我只得到一个浮点值，其中函数返回 nan，以及四个相邻的浮点值，其中函数返回“正确”值的两倍:

In [91]: phi = 0.9272952180016123

In [92]: mbad = 0.68359375000000011

In [93]: m = np.nextafter(mbad, 0)

In [94]: mvals = []

In [95]: for j in range(10):
   ....:     mvals.append(m)
   ....:     m = np.nextafter(m, 1)
   ....:     

In [96]: mvals
Out[96]: 
[0.68359375,
 0.68359375000000011,
 0.68359375000000022,
 0.68359375000000033,
 0.68359375000000044,
 0.68359375000000056,
 0.68359375000000067,
 0.68359375000000078,
 0.68359375000000089,
 0.683593750000001]

In [97]: ellipkinc(phi, mvals)
Out[97]: 
array([ 1.02593301,         nan,  2.05186602,  2.05186602,  2.05186602,
        2.05186602,  1.02593301,  1.02593301,  1.02593301,  1.02593301])