python之sin(x)/x 的数值积分_编程设计

我想使用 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

提前致谢!

请您参考如下方法:

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

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

这是你的代码修复(我在 python2 中运行它，这是我现在使用的 pc 中安装的):

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) 的 float 舍入)/1e-13 = 1!