CODEWITHSUNDEEP
pythonintegralsimpsons-rule
Filip Sokołowki
Filip Sokołowki

Reputation: 61

Simpson rule integration,Python

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

Upvotes: 2

Views: 7389

Answers (2)

user1984528
user1984528

Reputation: 453

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)

Upvotes: 1

Ami Tavory
Ami Tavory

Reputation: 76297

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.

Upvotes: 4

Related Questions