CODEWITHSUNDEEP
pythonscipynumerical-methods
xnx
xnx

Reputation: 25548

Bug in scipy.special.ellipkinc - incomplete elliptic integral

In using SciPy's scipy.special.ellipeinc and ellipkinc, there seem to be some islands of numerical instability. For example,

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

This happens where k^2.sin^2(phi) is close to 0.3, but there is nothing unusual about the elliptic integrals themselves here, so presumably it's a numerical thing. I don't know enough about the inner workings of this algorithm to say what's wrong, so what's my best option? I thought of: round(0.68359375000000011,8) say, but this is surely going to slow my code?

Upvotes: 3

Views: 598

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114921

(This is more an extended comment about the problem than an answer to what to do about it.)

This looks like a bug in ellipkinc. I get just one floating point value where the function returns nan, and four adjacent floating point values where the function returns twice the "correct" value:

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

Upvotes: 1

Related Questions