python 四边形积分似乎不准确
python quad integration seems inaccurate
我是 python 的新手,正在尝试对函数进行数值积分。一切似乎都有效,但我得到的结果与我在 Mathematica 中得到的结果有很大不同(我知道这是正确的)。谁能帮我弄清楚这是怎么回事?
代码如下:
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
如果我计算,例如,强度 (0),我在 python 中得到 0.0001580120220796804,在 Mathematica 中得到 0.000158898。在 0.5% 以内,所以这似乎还可以。但是,如果我计算 intensity(0.001),我在 python 中得到 1.8729902318383768e-05,在 Mathematica 中得到 0.00012034,它们相差近一个数量级。请注意,我已尝试减少绝对和相对误差,但这没有任何效果。
如有任何帮助,我们将不胜感激。
这是 Mathematica 代码:
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}]
然后
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
函数不包含pi
的因数。要解决此差异,请从 Python 代码的 sinc()
参数中删除 np.pi
。当我进行更改时,我得到的结果与 Mathematica 一致(我将 intensity()
修改为 return return 由 quad
编辑的错误:
In [12]: intensity(0)
Out[12]: (0.00015889773970382816, 1.764119291800849e-18)
In [13]: intensity(0.001)
Out[13]: (0.00012034021042092513, 1.3360447239754727e-18)
我是 python 的新手,正在尝试对函数进行数值积分。一切似乎都有效,但我得到的结果与我在 Mathematica 中得到的结果有很大不同(我知道这是正确的)。谁能帮我弄清楚这是怎么回事?
代码如下:
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
如果我计算,例如,强度 (0),我在 python 中得到 0.0001580120220796804,在 Mathematica 中得到 0.000158898。在 0.5% 以内,所以这似乎还可以。但是,如果我计算 intensity(0.001),我在 python 中得到 1.8729902318383768e-05,在 Mathematica 中得到 0.00012034,它们相差近一个数量级。请注意,我已尝试减少绝对和相对误差,但这没有任何效果。
如有任何帮助,我们将不胜感激。
这是 Mathematica 代码:
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}]
然后
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
函数不包含pi
的因数。要解决此差异,请从 Python 代码的 sinc()
参数中删除 np.pi
。当我进行更改时,我得到的结果与 Mathematica 一致(我将 intensity()
修改为 return return 由 quad
编辑的错误:
In [12]: intensity(0)
Out[12]: (0.00015889773970382816, 1.764119291800849e-18)
In [13]: intensity(0.001)
Out[13]: (0.00012034021042092513, 1.3360447239754727e-18)