3

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]
DJElectric
  • 349
  • 1
  • 4
  • 18

1 Answers1

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):

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

In [13]: intensity(0.001)
Out[13]: (0.00012034021042092513, 1.3360447239754727e-18)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214