python - Python 中的梯形法则

Question

我正在尝试在 Python 2.7.2 中实现梯形规则。我编写了以下函数:

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += h * f(a)
    for i in range(1, n):
        s += 2.0 * h * f(a + i*h)
    s += h * f(b)
    return s

但是，f(lambda x:x**2, 5, 10, 100) 返回 583.333(它应该返回 291.667)，很明显我的脚本有问题。虽然我无法发现它。

Answer 1

你的差距是原来的两倍。事实上，Trapezoidal Rule正如数学课上所教的那样，会使用像这样的增量

s += h * (f(a + i*h) + f(a + (i-1)*h))/2.0

(f(a + i*h) + f(a + (i-1)*h))/2.0 是对网格上两个相邻点的函数高度进行平均。

由于每两个相邻的梯形都有一条公共(public)边，因此上面的公式需要对函数进行两次必要的计算。

一个更有效的实现(更接近您发布的)，将结合 for-loop 相邻迭代的常用术语:

f(a + i*h)/2.0 + f(a + i*h)/2.0 =  f(a + i*h) 

到达:

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

print( trapezoidal(lambda x:x**2, 5, 10, 100))

产生

291.66875