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

consider following function:

import numpy np scipy.special import erf  def my_func(x):     return np.exp(x ** 2) * (1 + erf(x)) 

when evaluate integral of function -14 -4 using scipy's quad function, following result:

in [3]: scipy import integrate  in [4]: integrate.quad(my_func, -14, -4) /usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:289: userwarning: maximum number     of subdivisions (50) has been achieved.   if increasing limit yields no improvement advised analyze    integrand in order determine difficulties.  if position of    local difficulty can determined (singularity, discontinuity) 1    gain splitting interval , calling integrator    on subranges.  perhaps special-purpose integrator should used.   warnings.warn(msg) out[4]: (0.21896647054443383, 0.00014334175850538866) 

that is, 0.22.

however, when submit integral wolfram alpha, different result:

-5.29326 x 10 ^ 69. 

what's deal? i'm guessing has warning scipy has given me. what's best way evaluate integral in python?

note: increasing limit changes warning leaves 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: occurrence of roundoff error detected, prevents    requested tolerance being achieved.  error may    underestimated.   warnings.warn(msg) out[5]: (0.21894780966717864, 1.989164129832358e-05) 

tl;dr: integrand equivalent erfcx(-x), , implementation of erfcx @ scipy.special.erfcx takes care of numerical issues:

in [10]: scipy.integrate import quad  in [11]: 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 lambda expression changing sign of integration argument , limits:

in [14]: quad(erfcx, 4, 14) out[14]: (0.6990732491815446, 1.4463494884581349e-13) 

---

the problem numerical evaluation of 1 + erf(x) negative values of x. x decreases, erf(x) approaches -1. when add 1, catastrophic loss of precision, , sufficiently negative x (specifically x < -5.87), 1 + erf(x) numerically 0.

note default behavior @ wolfram alpha suffers same problem. had click on "more digits" twice reasonable answer.

the fix reformulate function. can express 1+erf(x) 2*ndtr(x*sqrt(2)), ndtr normal cumulative distribution function, available scipy.special.ndtr (see, example, https://en.wikipedia.org/wiki/error_function). here's alternative version of function, , result of integrating 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 @ wolfram alpha after clicking on "more digits" twice 0.6990732491815446...

this plot of function looks when use numerically stable version:

---

to avoid overflow or underflow arguments large magnitudes, can part of 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 @ali_m beat me in new question: tricking numpy/python representing large , small numbers.)

---

finally, simon byrne pointed out in answer on @ tricking numpy/python representing large , small numbers, function integrated can expressed erfcx(-x), erfcx scaled complementary error function. available scipy.special.erfcx.

for example,

in [10]: scipy.integrate import quad  in [11]: 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)