title: python quad integration seems inaccurate
body: <p>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?</p>
  
  <p>Here's the code:</p>
  
  <pre><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
  </code></pre>
  
  <p>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.</p>
  
  <p>Any help would be appreciated.</p>
  
  <p>Here is the Mathematica code:</p>
  
  <pre><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}]
  </code></pre>
  
  <p>and then</p>
  
  <pre><code>NumInt[0.0006, 0.000150, 0.000164, 0.8, 1.06, 0.0000006328, 0.001]
  </code></pre>