Python 中的辛普森法则

对于数值方法类，我需要编写一个程序来计算使用辛普森复合规则的定积分。我已经做到了这一点（见下文），但我的答案不正确。我正在使用 f(x)=x 测试该程序，在 0 到 1 之间积分，结果应为 0.5。我得到 0.78746...等等。 我知道 Scipy 中有辛普森规则，但我真的需要自己编写它。

我怀疑这两个循环有问题。我之前尝试过“for i in range(1, n, 2)”和“for i in range(2, n-1, 2)”，这给了我 0.41668333... 等的结果。 我还尝试了“x += h”，我尝试了“x += i*h”。第一个选项给了我 0.3954，第二个选项给了我 7.9218。

# Write a program to evaluate a definite integral using Simpson's rule with
# n subdivisions

from math import *
from pylab import *

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a
    for i in range(1,n/2):
        x += 2*h
        k += 4*f(x)
    for i in range(2,(n/2)-1):
        x += 2*h
        k += 2*f(x)
    return (h/3)*(f(a)+f(b)+k)

def function(x): return x

print simpson(function, 0.0, 1.0, 100)

您可能忘记初始化x在第二次循环之前，起始条件和迭代次数也已关闭。这是正确的方法：

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a + h
    for i in range(1,n/2 + 1):
        k += 4*f(x)
        x += 2*h

    x = a + 2*h
    for i in range(1,n/2):
        k += 2*f(x)
        x += 2*h
    return (h/3)*(f(a)+f(b)+k)

您的错误与循环不变式的概念有关。不要过多讨论细节，通常更容易理解和调试在周期结束时推进的周期，而不是在开始时，这里我移动了x += 2 * h行到末尾，这使得很容易验证求和的开始位置。在您的实现中，有必要分配一个奇怪的x = a - h对于第一个循环仅添加2 * h将其作为循环中的第一行。