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),所以很明显我的脚本有问题.我不能发现它.
你离开了两倍.实际上,数学课中教授的梯形法则会使用增量
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)