辛普森在Python中的规则

Question

对于数值方法类,我需要编写一个程序来评估Simpson复合规则的定积分.我已经走到这一步了(见下文),但我的回答并不正确.我用f(x)= x测试程序,积分在0到1之间,结果应为0.5.我得到了0.78746 ...等我知道Scipy有一个Simpson规则,但我真的需要自己编写.

我怀疑这两个循环有问题.我之前试过"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)

Answer 1

您可能忘记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其添加为循环中的第一行。