Bug in SciPy’s erf function

Last night I produced the plot below and was very surprised at the jagged spike. I knew the curve should be smooth and strictly increasing.

My first thought was that there must be a numerical accuracy problem in my code, but it turns out there’s a bug in SciPy version 0.8.0b1. I started to report it, but I saw there were similar bug reports and one such report was marked as closed, so presumably the fix will appear in the next release.

The problem is that SciPy’s erf function is inaccurate for arguments with imaginary part near 5.8. For example, Mathematica computes erf(1.0 + 5.7i) as -4.5717×1012 + 1.04767×1012 i. SciPy computes the same value as -4.4370×1012 + 1.3652×1012 i. The imaginary component is off by about 30%.

Here is the code that produced the plot.

from scipy.special import erf
from numpy import linspace, exp
import matplotlib.pyplot as plt

def g(y):
    z = (1 + 1j*y) /  sqrt(2)
    temp = exp(z*z)*(1 - erf(z))
    u, v = temp.real, temp.imag
    return -v / u

x = linspace(0, 10, 101)
plt.plot(x, g(x))
Tagged with: , , ,
Posted in Python, Software development
3 comments on “Bug in SciPy’s erf function
  1. Ralf says:

    In [31]: sp.__version__
    Out[31]: ’0.9.0.dev’

    In [32]: sp.special.erf(1.0 + 5.7j)*1e-12
    Out[32]: (-4.5717045780553551+1.0476748318787288j)

  2. pv says:

    The fix is also included in the current Scipy 0.8.0 release.

  3. Meziane says:

    BONJOUR,
    erf(z) =int(exp-t2dt) 0<t<x (1)
    erf(z)=1-[1+0.278393z+0.230389z2+0.000972z3+0.078108z4]puiss-4 +e(z) (2)
    jai pas compris comment trouve (2) aprtir de (1)
    marci

1 Pings/Trackbacks for "Bug in SciPy’s erf function"
  1. [...] This post was mentioned on Twitter by Planet Python, Alltop Programming. Alltop Programming said: Bug in SciPy’s erf function http://bit.ly/cb9SPk [...]