Python中的Simpson规则

问题内容：

对于数值方法类，我需要编写一个程序以使用Simpson的合成规则来评估定积分。我已经走了这么远（见下文），但是我的答案不正确。我正在测试f（x）=
x的程序，该程序在0到1之间进行积分，结果应为0.5。我得到0.78746
…等。我知道Scipy中有可用的Simpson规则，但是我真的需要自己编写它。

我怀疑这两个循环有问题。之前，我尝试过“ for i in range（1，n，2）”和“ for 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为循环中的第一行。

Python中的Simpson规则

微信关注