Python中的梯形规则

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/21146540/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-18 22:10:34  来源:igfitidea点击:

Trapezoidal rule in Python

pythonnumerical-integration

提问by user3199345

I'm trying to implement the trapezoidal rule in Python 2.7.2. I've written the following function:

我正在尝试在 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

However, f(lambda x:x**2, 5, 10, 100) returns 583.333 (it's supposed to return 291.667), so clearly there is something wrong with my script. I can't spot it though.

但是, f(lambda x:x**2, 5, 10, 100) 返回 583.333(它应该返回 291.667),所以很明显我的脚本有问题。不过我看不出来。

回答by unutbu

You are off by a factor of two. Indeed, the Trapezoidal Ruleas taught in math class would use an increment like

你错了两倍。实际上,数学课中教授的梯形规则将使用类似的增量

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

(f(a + i*h) + f(a + (i-1)*h))/2.0is averaging the height of the function at two adjacent points on the grid.

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

Since every two adjacent trapezoids have a common edge, the formula above requires evaluating the function twice as often as necessary.

由于每两个相邻的梯形都有一个公共边,因此上面的公式需要将函数计算为必要的两倍。

A more efficient implementation (closer to what you posted), would combine common terms from adjacent iterations of the for-loop:

更有效的实现(更接近您发布的内容)将结合来自相邻迭代的公共术语for-loop

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

to arrive at:

到达:

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))

which yields

这产生

291.66875

回答by abarnert

The trapezoidal rulehas a big /2fraction (each term is (f(i) + f(i+1))/2, not f(i) + f(i+1)), which you've left out of your code.

梯形规则有一个大的/2分数(每学期是(f(i) + f(i+1))/2,不是f(i) + f(i+1)),你已经离开你的代码了。

You've used the common optimization that treats the first and last pair specially so you can use 2 * f(i)instead of calculating f(i)twice (once as f(j+1)and once as f(i)), so you have to add the / 2to the loop step and to the special first and last steps:

您已经使用了特殊处理第一对和最后一对的通用优化,因此您可以使用2 * f(i)而不是计算f(i)两次(一次 asf(j+1)和一次 as 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

You can obviously simplify the loop step by replacing the 2.0 * … / 2.0with just .

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

However, even more simply, you can just divide the whole thing by 2 at the end, changing nothing but this line:

然而,更简单的是,你可以在最后将整个事物除以 2,只改变这一行:

return s / 2.0