Python中的梯形规则答案

【问题标题】：Trapezoidal rule in PythonPython中的梯形规则
【发布时间】：2014-02-04 11:02:42
【问题描述】：

我正在尝试在 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），所以很明显我的脚本有问题。不过我看不出来。

【问题讨论】：

可能并非巧合，583.333 / 2 = 291.6665。
提示：您返回的正好是应返回的两倍。 trapezoidal rule 公式中有一个 / 2 分数。你的代码没有。

标签： python numerical-integration

【解决方案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 对网格上两个相邻点的函数高度进行平均。

由于每两个相邻的梯形都有一条公共边，因此上面的公式需要根据需要对函数进行两次评估。

更有效的实现（更接近您发布的内容）将结合 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

【讨论】：

没有。您在此处两次除以二（删除 2.0 * 并添加 / 2.0，但您忘记将最后一项除以二，并且您还使用了错误的规则形式。
不过，如果他这样做，他将不得不更改更多代码，而且效率会低于他现在所做的，因为它需要多个函数评估。

【解决方案2】：

trapezoidal rule 有一个很大的 /2 分数（每个术语是 (f(i) + f(i+1))/2，而不是 f(i) + f(i+1)），这是您在代码中遗漏的。

您已经使用了专门处理第一对和最后一对的通用优化，因此您可以使用 2 * f(i) 而不是计算 f(i) 两次（一次作为 f(j+1) 和一次作为 f(i)），所以你必须将/ 2 添加到循环步骤以及特殊的第一步和最后一步：

s += h * f(a) / 2.0
for i in range(1, n):
    s += 2.0 * h * f(a + i*h) / 2.0
s += h * f(b) / 2.0

您显然可以通过将2.0 * … / 2.0 替换为… 来简化循环步骤。

但是，更简单的是，您可以在最后将整个内容除以 2，只改变这一行：

return s / 2.0

【讨论】：