2

我写了这段代码，但我不确定它是否正确。在辛普森规则中，有一个条件是它必须有偶数个间隔。我不知道如何将这个条件印在我的代码中。

def simpson(data):
    data = np.array(data)
    a = min(range(len(data)))
    b = max(range(len(data)))
    n = len(data)
    h = (b-a)/n
    for i in range(1,n, 2):
        result += 4*data[i]*h
    for i in range(2,n-1, 2):
        result += 2*data[i]*h
    return result * h /3
4

2 回答 2

4

有趣的是，您可以在Wikipedia 条目中找到它：

from __future__ import division  # Python 2 compatibility

def simpson(f, a, b, n):
    """Approximates the definite integral of f from a to b by the
    composite Simpson's rule, using n subintervals (with n even)"""

    if n % 2:
        raise ValueError("n must be even (received n=%d)" % n)

    h = (b - a) / n
    s = f(a) + f(b)

    for i in range(1, n, 2):
        s += 4 * f(a + i * h)
    for i in range(2, n-1, 2):
        s += 2 * f(a + i * h)

    return s * h / 3

你在哪里使用它：

simpson(lambda x:x**4, 0.0, 10.0, 100000)

请注意它是如何通过需要一个函数和n.

但是，如果您需要它作为值列表，那么在调整代码之后（这应该很容易），我建议您也提出 aValueError以防它的长度不均匀。

于 2016-01-27T23:52:31.013 回答
1

由于您似乎已经在使用 numpy，您也可以考虑使用 scipy，它可以方便地提供 Simpson 的规则集成例程。

from scipy.integrate import simps
result=simps(data)

有关完整文档，请参阅http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.simps.html（他们讨论偶数/奇数间隔的处理）

于 2016-01-28T00:00:54.997 回答