scipy特殊椭圆函数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时，就会发生这种情况，但椭圆积分本身并没有什么异常，所以这可能是一个数值问题。我对这个算法的内部工作原理了解不多，所以我最好的选择是什么？我想到了： round(0.68359375000000011,8) 但是这肯定会让我的代码变慢？在

Answer 1

（这与其说是对问题的回答，不如说是对问题的延伸评论。）

这看起来像是ellipkinc中的一个bug。我只得到一个浮点值，其中函数返回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])