Python中的梯形规则

use*_*345 3 python numerical-integration

我正在尝试在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
Run Code Online (Sandbox Code Playgroud)

但是,f(lambda x:x**2,5,10,100)返回583.333(它应该返回291.667),所以很明显我的脚本有问题.我不能发现它.

unu*_*tbu 5

你离开了两倍.实际上,数学课中教授的梯形法则会使用增量

s += h * (f(a + i*h) + f(a + (i-1)*h))/2.0
Run Code Online (Sandbox Code Playgroud)

(f(a + i*h) + f(a + (i-1)*h))/2.0 平均网格上两个相邻点处的函数高度.

由于每两个相邻的梯形具有共同的边缘,因此上述公式需要根据需要两次评估该函数.

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

f(a + i*h)/2.0 + f(a + i*h)/2.0 =  f(a + i*h) 
Run Code Online (Sandbox Code Playgroud)

到达:

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))
Run Code Online (Sandbox Code Playgroud)

产量

291.66875
Run Code Online (Sandbox Code Playgroud)