sin(x)/x 的数值积分

Question

我想使用 scipy 计算 python 中 sin(x)/x 的求积的定积分。 n = 256。它似乎工作不正常：

 from scipy import integrate

 exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
 print("Exact value of integral:", exact)

 # Approx of sin(x)/x by Trapezoidal rule
 x = np.linspace(0, 2*np.pi, 257)
 f = lambda x : (np.sin(x))/x
 approximation = np.trapz(f(x), x)
 print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
   '+', round((2/256)**2, 5))
 print ("real error:", exact - approximation)

 # Approx of sin(x)/x by Simpsons 1/3 rule
 approximation = integrate.simps(f(x), x)
 print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
   '+', round((2/256)**4, 9))
 print ("real error:", exact - approximation)

 plt.figure()
 plt.plot(x, f(x))
 plt.show()

精确的积分计算得很好，但求积时出现错误...有什么问题吗？

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): nan + 6e-05
real error: nan
Estimated value of simpsons O(h^4): nan + 4e-09
real error: nan

提前致谢！

Answer 1

您的代码中至少存在一些问题：

您的 linspace 从 0 开始，因此当您计算要积分的函数时，在梯形积分的开头，您有：sin(0)/0 = nan .您应该使用 数字零 而不是 精确零 （在下面的示例中我使用了 1e-12)
当你得到第一个 nan、nan + 1.0 = nan 时：这意味着在你的代码中，当你对区间求和时，第一个 nan 搞砸了你所有的结果。
仅适用于python 2：除法2/256是2个整数之间的除法，结果为0。尝试使用 2.0/256.0（感谢@MaxU 指出这一点）。

这是你的代码修复（我运行它在python2，这是我现在使用的电脑上安装的）：

from scipy import integrate
import numpy as np

exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
print("Exact value of integral:", exact)

# Approx of sin(x)/x by Trapezoidal rule
x = np.linspace(1e-12, 2*np.pi, 257) # <- 0 has become a numeric 0
f = lambda x : (np.sin(x))/x
approximation = np.trapz(f(x), x)
print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
  '+', round((2.0/256.0)**2, 5))
print ("real error:", exact - approximation)

# Approx of sin(x)/x by Simpsons 1/3 rule
approximation = integrate.simps(f(x), x)
print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
  '+', round((2/256)**4, 9))
print ("real error:", exact - approximation)

及其输出：

('Exact value of integral:', 1.4181515761326284)
('Estimated value of trapezoidal O(h^2):', 1.41816, '+', 6e-05)
('real error:', -7.9895502944626884e-06)
('Estimated value of simpsons O(h^4):', 1.418151576, '+', 0.0)
('real error:', 2.7310242955991271e-10)

Discalimer sin(x)/x -> 1 for x -> 0 的极限，但由于 sin(1e-12)/1e-13 = 1!

的浮动四舍五入

Answer 2

NaN 表示 "Not A Number"。在你的情况下基本上是无穷大。当您创建时：

x = np.linspace(0, 2*np.pi, 257)

你创建了一个数组，其中的值为 0，然后你尝试除以 x，但你不能除以 0...

一个解决方案是使用这个：

x = np.linspace(0.1, 2*np.pi, 257)

哪个给你这个:

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): 1.31822 + 6e-05
real error: 0.099935104987
Estimated value of simpsons O(h^4): 1.318207115 + 4e-09
real error: 0.0999444614012

越接近零，近似值就越好！

Answer 3

对于 x == 0，您可以使函数 return 1（sin(x)/x 在 0 中的极限）而不是 NaN。这样，您就没有欺骗并更改您积分的间隔以排除 0.

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
print("Exact value of integral:", exact)


def f(x):
    out = np.sin(x) / x
    # For x == 0, we get nan. We replace it by the 
    # limit of sin(x)/x in 0
    out[np.isnan(out)] = 1
    return out

# Approx of sin(x)/x by Trapezoidal rule

x = np.linspace(0, 2*np.pi, 257)

approximation = np.trapz(f(x), x)
print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
  '+', round((2/256)**2, 5))
print ("real error:", exact - approximation)
 # Approx of sin(x)/x by Simpsons 1/3 rule
approximation = integrate.simps(f(x), x)
print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
  '+', round((2/256)**4, 9))
print ("real error:", exact - approximation)
plt.figure()
plt.plot(x, f(x))
plt.show()

输出：

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): 1.41816 + 6e-05
real error: -7.98955129322e-06
Estimated value of simpsons O(h^4): 1.418151576 + 4e-09
real error: 2.72103006793e-10

Answer 4

scipy.special.sici 将 return 正弦和余弦积分。

sin(x)/x 的数值积分

Numeric integral of sin(x)/x

python

numpy

integral

scipy

python-3.x