Wolfram Alpha and scipy.integrate.quad give me different answers for the same integral

Consider the following function:

import numpy as np
from scipy.special import erf

def my_func(x):
    return np.exp(x ** 2) * (1 + erf(x))

When I evaluate the integral of this function from -14 to -4 using scipy's quad function, I get the following result:

In [3]: from scipy import integrate

In [4]: integrate.quad(my_func, -14, -4)
/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:289: UserWarning: The maximum number     of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg)
Out[4]: (0.21896647054443383, 0.00014334175850538866)

That is, about 0.22.

However, when I submit this integral to Wolfram Alpha, I get a very different result:

-5.29326 X 10 ^ 69.

What's the deal? I'm guessing this has to do with the warning scipy has given me. What's the best way to evaluate this integral in python?

NOTE: Increasing the limit changes the warning but leaves the scipy result unchanged:

In [5]: integrate.quad(my_func, -14, -4, limit=10000)
/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:289: UserWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  warnings.warn(msg)
Out[5]: (0.21894780966717864, 1.989164129832358e-05)

Solution

TL;DR: The integrand is equivalent to erfcx(-x), and the implementation of erfcx at scipy.special.erfcx takes care of the numerical issues:

In [10]: from scipy.integrate import quad

In [11]: from scipy.special import erfcx

In [12]: quad(lambda x: erfcx(-x), -14, -4)
Out[12]: (0.6990732491815446, 1.4463494884581349e-13)

In [13]: quad(lambda x: erfcx(-x), -150, -50)
Out[13]: (0.6197754761443759, 4.165648376274775e-14)

You can avoid the lambda expression by changing the sign of the integration argument and limits:

In [14]: quad(erfcx, 4, 14)
Out[14]: (0.6990732491815446, 1.4463494884581349e-13)

The problem is the numerical evaluation of 1 + erf(x) for negative values of x. As x decreases, erf(x) approaches -1. When you then add 1, you get catastrophic loss of precision, and for sufficiently negative x (specifically x < -5.87), 1 + erf(x) is numerically 0.

Note that the default behavior at Wolfram Alpha suffers from the same problem. I had to click on "More digits" twice to get a reasonable answer.

The fix is to reformulate your function. You can express 1+erf(x) as 2*ndtr(x*sqrt(2)), where ndtr is the normal cumulative distribution function, available from scipy.special.ndtr (see, for example, https://en.wikipedia.org/wiki/Error_function). Here's an alternative version of your function, and the result of integrating it with scipy.integrate.quad:

In [133]: def func2(x):
   .....:     return np.exp(x**2) * 2 * ndtr(x * np.sqrt(2))
   .....: 

In [134]: my_func(-5)
Out[134]: 0.1107029852258767

In [135]: func2(-5)
Out[135]: 0.11070463773306743

In [136]: integrate.quad(func2, -14, -4)
Out[136]: (0.6990732491815298, 1.4469372263470424e-13)

The answer at Wolfram Alpha after clicking on "More digits" twice is 0.6990732491815446...

This is what the plot of the function looks like when you use a numerically stable version:

To avoid overflow or underflow for arguments with very large magnitudes, you can do part of the computation in log-space:

from scipy.special import log_ndtr

def func3(x):
    t = x**2 + np.log(2) + log_ndtr(x * np.sqrt(2))
    y = np.exp(t)
    return y

E.g.

In [20]: quad(func3, -150, -50)
Out[20]: (0.6197754761435517, 4.6850379059597266e-14)

(Looks like @ali_m beat me to it in the new question: Tricking numpy/python into representing very large and very small numbers.)

Finally, as Simon Byrne pointed out in an answer over at Tricking numpy/python into representing very large and very small numbers, the function to be integrated can be expressed as erfcx(-x), where erfcx is the scaled complementary error function. It is available as scipy.special.erfcx.

For example,

In [10]: from scipy.integrate import quad

In [11]: from scipy.special import erfcx

In [12]: quad(lambda x: erfcx(-x), -14, -4)
Out[12]: (0.6990732491815446, 1.4463494884581349e-13)

In [13]: quad(lambda x: erfcx(-x), -150, -50)
Out[13]: (0.6197754761443759, 4.165648376274775e-14)