[FIXED] Numeric integral of sin(x)/x

Issue

I want to calcultate the definite integral with quadratures of sin(x)/x in python using scipy. With n = 256. It doesnt seem to work well:

 from scipy import integrate

 exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
 print("Exact value of integral:", exact)

 # Approx of sin(x)/x by Trapezoidal rule
 x = np.linspace(0, 2*np.pi, 257)
 f = lambda x : (np.sin(x))/x
 approximation = np.trapz(f(x), x)
 print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
   '+', round((2/256)**2, 5))
 print ("real error:", exact - approximation)

 # Approx of sin(x)/x by Simpsons 1/3 rule
 approximation = integrate.simps(f(x), x)
 print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
   '+', round((2/256)**4, 9))
 print ("real error:", exact - approximation)

 plt.figure()
 plt.plot(x, f(x))
 plt.show()

The exact integral is calculated well, but I get an error with quadratures... What is wrong?

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): nan + 6e-05
real error: nan
Estimated value of simpsons O(h^4): nan + 4e-09
real error: nan

Thanks in advance!

Solution

There are at least some issues in your code:

1. Your linspace starts from 0, thus when you evaluate the function to integrate, at the beginning of the trapezoidal integral, you have: sin(0)/0 = nan. You should use a numeric zero instead of an exact zero (in the example below I used 1e-12)
2. When you get the first nan, nan + 1.0 = nan: this means that in your code, when you are summing up the integral over the intervals, the first nan is messing up all of your results.
3. for python 2 only: The division 2/256 is a division between 2 integers and the result is 0. Try with 2.0/256.0 instead (thanks @MaxU for pointing that out).

This is your code fixed (I'm running it in python2, this is what is installed in the pc I'm using now):

from scipy import integrate
import numpy as np

exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
print("Exact value of integral:", exact)

# Approx of sin(x)/x by Trapezoidal rule
x = np.linspace(1e-12, 2*np.pi, 257) # <- 0 has become a numeric 0
f = lambda x : (np.sin(x))/x
approximation = np.trapz(f(x), x)
print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
  '+', round((2.0/256.0)**2, 5))
print ("real error:", exact - approximation)

# Approx of sin(x)/x by Simpsons 1/3 rule
approximation = integrate.simps(f(x), x)
print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
  '+', round((2/256)**4, 9))
print ("real error:", exact - approximation)

with its output:

('Exact value of integral:', 1.4181515761326284)
('Estimated value of trapezoidal O(h^2):', 1.41816, '+', 6e-05)
('real error:', -7.9895502944626884e-06)
('Estimated value of simpsons O(h^4):', 1.418151576, '+', 0.0)
('real error:', 2.7310242955991271e-10)

Discalimer The limit of sin(x)/x -> 1 for x -> 0, but due to floating rounding for sin(1e-12)/1e-13 = 1!

Answered By - Matteo Ragni