CODEWITHSUNDEEP
pythonscipyintegral
bin
bin

Reputation: 39

Calculating Dirichlet Integral in Python

I'm trying to calculate the integral

sin(x)/x , x = [0,inf]

I do the following :

import math
from scipy.integrate import quad

t = float("inf")
def integrand(x):
    return (math.sin(x))/x
ans, err = quad(integrand, 0, t)

print(ans)

However I get the wrong answer: It should be Pi/2 and the answer from Python is 2.247.

Do you know what may be the problem?

Upvotes: 1

Views: 498

Answers (1)

JDong
JDong

Reputation: 2304

You're calculating a difficult integral -- the Dirichlet integral isn't even defined unless you use very specific definitions of the integral. https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html the function you're using -- quad, is for General purpose integration. Take into account that it's using an approximation algorithm.

Did you look at the warning? My output shows:

/Users/jdong/venv2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning)

... and my error is 3.2903230524472544. So technically PI/2 is within the error bounds.

Also, you're calculating to infinity, where it's not really necessary. The sinc function tends to zero moderately quickly, so you can simply do quad(integrand,0,100) and get a decent result. Note that as you increase the limit however, you will need to increase the number of subdivisions used by quad. e.g. quad(integrand,0,100000,limit=10000)

Upvotes: 4

Related Questions