python quad integration seems inaccurate

python quad integration seems inaccurate



I'm a bit of a newbie to python and am trying to numerically integrate a function. Everything seems to work, but the results I'm getting vary significantly from what I get in Mathematica (which I know to be correct). Can someone help me figure out what's going on?



Here's the code:


def integrand(x, d, a, b, l, s, wavelength, y):
return b*(np.sinc((np.pi*a/(wavelength*s))*(y + s*b*x/l))**2)*np.cos((np.pi*d/(wavelength*s))*(y + s*b*x/l))**2


def intensity(y):
result, error = si.quad(integrand, -1/2, 1/2, epsrel = 1e-16, epsabs = 1e-16,
args=(0.0006, 0.000150, 0.000164, 0.8, 1.06, 0.0000006328, y))
return result



If I calculate, for example, intensity(0), I get 0.0001580120220796804 in python and 0.000158898 in Mathematica. Within 0.5%, so that seems okay. But if I calculate intensity(0.001) I get 1.8729902318383768e-05 in python and 0.00012034 in Mathematica, which are different by nearly an order of magnitude. Note that I have tried reducing the absolute and relative errors, but this doesn't have any effect.



Any help would be appreciated.



Here is the Mathematica code:


NumInt[d_, a_, b_, l_, s_, lambda_, y_] := NIntegrate[b Sinc[(a Pi/(s lambda)) (y - (s*b*
x/l))]^2 Cos[(d Pi/(s lambda)) (y - (s*b*x/l))]^2, x, -1/2,
1/2]



and then


NumInt[0.0006, 0.000150, 0.000164, 0.8, 1.06, 0.0000006328, 0.001]




1 Answer
1



numpy.sinc is defined as sin(pi*x)/(pi*x). Mathematica's Sinc function does not include the factor of pi. To fix the discrepancy, remove np.pi from the sinc() argument in the Python code. When I make that change, I get results that agree with Mathematica (I modified intensity() to also return the error that is returned by quad):


numpy.sinc


sin(pi*x)/(pi*x)


Sinc


pi


np.pi


sinc()


intensity()


quad


In [12]: intensity(0)
Out[12]: (0.00015889773970382816, 1.764119291800849e-18)

In [13]: intensity(0.001)
Out[13]: (0.00012034021042092513, 1.3360447239754727e-18)





Brilliant! I should have thought to check the sinc definition. Thanks so much!!
– DJElectric
Aug 22 at 18:43






By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.