我想编写一个自定义函数来对具有特定公差的表达式(python 或 lambda 函数)进行数值积分。我知道使用 scipy.integrate.quad 可以简单地更改 epsabs 但我想使用 numpy 自己编写函数。
来自 this blogpost我知道函数:
def integrate(f, a, b, N):
x = np.linspace(a+(b-a)/(2*N), b-(b-a)/(2*N), N)
fx = f(x)
area = np.sum(fx)*(b-a)/N
return area
给出了 N 段的数值积分。我如何编写另一个函数或扩展此函数以获取 tol 输入并增加 N，直到两个后续计算之间的差异小于给定的公差？
最佳答案
使用您拥有的函数，可以从合理的 N(例如 5)开始，并不断将数字加倍，直到达到所需的精度。
def integrate_tol(f, a, b, tol):
N = 5
old_integral = integrate(f, a, b, N)
while True:
N *= 2
new_integral = integrate(f, a, b, N)
if np.abs(old_integral - new_integral) < tol:
return (4*new_integral - old_integral)/3
old_integral = new_integral
一个简单的测试:
f = lambda x: np.exp(x)
print(integrate_tol(f, -1, 1, 1e-9))
print(np.exp(1)-np.exp(-1)) # exact value for comparison
打印
2.3504023872876028
2.3504023872876028
不能保证错误确实小于 tol(但话又说回来，scipy.quad 也不能保证)。在实践中，误差会比 tol 小得多，因为我使用了技巧，称为 Richardson extrapolation :返回值 (4*new_integral - old_integral)/3 通常比新旧近似值本身准确得多。 (说明:由于 integrate 使用中点规则，N 的每次加倍都会将误差减少大约 4 倍。因此，采用组合 4*new_integral - old_integral 几乎抵消了排除这两个结果中的残余误差。)
(注意:在while循环中，不建议从N=1开始；它可能还不够，并且由于某些数字巧合而过早停止的风险更高，例如函数在一堆中为零地方。)
关于python - 使用 NumPy 编写一个函数来计算具有特定公差的积分，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/50033275/