Simpson rule integration,Python

Question

I wrote this code,but I not sure if it is right.In Simpson rule there is condition that it has to has even number of intervals.I dont know how to imprint this condition into my code.

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

Answer 1

Since you already seem to be using numpy you may also consider using scipy which conveniently provides a Simpson's rule integration routine.

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

See http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.simps.html for the full documentation (where they discuss the handling of even/odd intervals)

Answer 2

Interestingly enough, you can find it in the Wikipedia entry:

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

where you use it as:

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

---

Note how it bypasses your parity problem by requiring a function and n.

In case you need it for a list of values, though, then after adapting the code (which should be easy), I suggest you also raise a ValueError in case its length is not even.