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]
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):
In [12]: intensity(0)
Out[12]: (0.00015889773970382816, 1.764119291800849e-18)
In [13]: intensity(0.001)
Out[13]: (0.00012034021042092513, 1.3360447239754727e-18)