CODEWITHSUNDEEP
python
phchen
phchen

Reputation: 123

Python integrate with infinite

Sometimes, I get wrong solution when I integrate with infinite boundaries in Python. Here is a simple example to illustrate my confusion:

from scipy import *
import numpy as np
from scipy.integrate import quad

def integrand1(x):
    output = exp(-(x-1.0)**2.0)
    return output

def integrand2(x):
    output = exp(-(x-100.0)**2.0)
    return output

solution1 = quad(integrand1,-np.inf,np.inf)
print solution1

solution1 = quad(integrand2,-np.inf,np.inf)
print solution2

while the output is:

(1.7724538509055159, 3.668332157626072e-11)
(0.0, 0.0)

I don't understand why the second integral is wrong while the first one avoid to be wrong. It will be great to tell me some tricks to handle infinite in Python.

Upvotes: 2

Views: 18928

Answers (1)

themiurge
themiurge

Reputation: 1669

There is nothing inherently wrong with your code. The results you get are due to the quad algorithm being an approximate method whose accuracy, from what I gathered doing some tests, greatly depends on where the midpoint of the integration interval is located, with respect to the x-axis interval where the integrand is significantly different from 0.

The midpoint of the integration interval in case of a (-inf,+inf) interval is always 0 (see comment to relevant Fortran code here, starting at line 238), and (sadly) cannot be configured. Your integrand2 function is centered on x=100, which is too far from the quad algorithm's midpoint for it to be accurate enough.

It would be nice to be able to specify the midpoint in case of integration between -inf and +inf, but the good news is that you can implement it yourself with function decorators. First you need a wrapper for your integrand function, in order to shift it arbitrarily over the x axis:

def shift_integrand(integrand, offset):
    def dec(x):
        return integrand(x - offset)
    return dec

This generates a new function based on any integrand you like, just shifting it along the x axis according to the offset parameter. So if you do something like this (using your integrand1 and integrand2 functions):

new_integrand1 = shift_integrand(integrand1, -1.0)
print new_integrand1(0.0)
new_integrand2 = shift_integrand(integrand2, -100.0)
print new_integrand2(0.0)

you get:

1.0
1.0

Now you need another wrapper for the quad function, in order to be able to pass a new midpoint:

def my_quad(func, a, b, midpoint=0.0, **kwargs):
    if midpoint != 0.0:
        func = shift_integrand(func, -midpoint)
    return quad(func, a, b, **kwargs)

Finally, knowing where your integrands are centered, you can call it thus:

solution2 = my_quad(integrand2, -np.inf, np.inf, midpoint=100.0)
print solution2

Which yields:

(1.772453850905516, 1.4202639944499085e-08)

Upvotes: 6

Related Questions