使用NumPy编写一个函数来计算具有特定容差的积分

Question

使用NumPy编写一个函数来计算具有特定容差的积分

3

我希望编写一个自定义函数，用特定的容差数值对表达式（Python或Lambda函数）进行数值积分。我知道使用scipy.integrate.quad可以简单地更改epsabs，但我想使用numpy自己编写该函数。

从这篇博客文章中我了解到这个函数：

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，直到两次计算之间的差异小于给定的容差？

- Foad S. Farimani

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

使用您拥有的函数，可以从一个合理的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也无法保证）。实际上，由于我使用了称为Richardson外推的技巧，因此误差将远小于tol：返回值(4*new_integral - old_integral)/3通常比新旧逼近本身更准确。（解释：由于integrate使用中点规则，N加倍会将误差减少约4倍。因此，取组合4*new_integral - old_integral几乎会消除这些结果中的残余误差。）

（备注：在while循环中，从N = 1开始不可取；它可能不够用，并且由于某些数值巧合，例如函数在一堆地方为零，停止太早的风险更高。）