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)
Comments
Post a Comment