Integration with Scipy giving incorrect results with negative lower bound

Question

I am attempting to calculate integrals between two limits using python/scipy.

I am using online calculators to double check my results (http://www.wolframalpha.com/widgets/view.jsp?id=8c7e046ce6f4d030f0b386ea5c17b16a, http://www.integral-calculator.com/), and my results disagree when I have certain limits set.

The code used is:

import scipy as sp
import numpy as np

def integrand(x):
    return np.exp(-0.5*x**2)

def int_test(a,b):
    # a and b are the lower and upper bounds of the integration
    return sp.integrate.quad(integrand,a,b)

When setting the limits (a,b) to (-np.inf,1) I get answers that agree (2.10894...) however if I set (-np.inf,300) I get an answer of zero.

On further investigation using:

for i in range(50):
    print(i,int_test(-np.inf,i))

I can see that the result goes wrong at i=36.

I was wondering if there was a way to avoid this?

Thanks,

Matt

Answer 1

I am guessing this has to do with the infinite bounds. scipy.integrate.quad is a wrapper around quadpack routines.

https://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html

In the end, these routines chose suitable intervals and try to get the value of the integral through function evaluations and then numerical integrations. This works fine for finite integrals, assuming you know roughly how fine you can make the steps of the function evaluation.

For infinite integrals it depends how well the algorithms choose respective subintervals and how accurately they are computed.

My advice: do NOT use numerical integration software AT ALL if you are interested in accurate values for infinite integrals.

If your problem can be solved analytically, try that or confine yourself to certain bounds.