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

you should trying printing Min, range, len. One of those seem to be a variable of type List rather than a function — Garrett R, Jan 27 '16 at 23:39
Did you mean `min()`? The builtin for getting the smallest integer from a collection? — kylieCatt, Jan 27 '16 at 23:41
this `a = min(range(len(data)))` gives always `a = 0` and this `b = max(range(len(data)))` gives always `b = len(data)-1` - maybe you try to do something different. — furas, Jan 27 '16 at 23:46

score 4 · Answer 1 · answered Jan 27 '16 at 23:52

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.

score 1 · Answer 2 · answered Jan 28 '16 at 00:00

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)

Simpson rule integration,Python

2 Answers2