Simpson's Rule using for loops (numerical integration)

Question

I am trying to code Simpson's Rule in python using for loops and I keep getting an assertion error and cant find out why.

def integrate_numeric(xmin, xmax, N):
    xsum = 0
    msum = 0
    h = (xmax-xmin)//N

    for i in range(0, N):
        xsum += f(xmin + i*h)
        print (xsum)

    for i in range(0,N-1):
        msum += f(xmin + (h/2) + i*h)    
        print (msum)

    I = (h/6) * (f(xmin) + 4*(msum) + 2*(xsum) + f(xmax))
    return I

f:

def f(x):
    return (x**2) * numpy.sin(x)

sample:

assert numpy.isclose(integrate_numeric(xmin=0, xmax=4, N=50), 1.096591)

Answer 1

Here's a fixed version of your code:

import numpy

def integrate_numeric(xmin, xmax, N):
    '''                                                                                                                               
    Numerical integral of f from xmin to xmax using Simpson's rule with                                                               
        N panels.                                                                                                                     
    '''
    xsum = 0
    msum = 0
    h = (xmax-xmin)/N

    for i in range(1, N):
        xsum += f(xmin + i*h)
        print(xsum)

    for i in range(0, N):
        msum += f(xmin + (h/2) + i*h)
        print(msum)

    I = (h/6) * (f(xmin) + 4*msum + 2*xsum + f(xmax))
    return I


def f(x):
    '''Function equivalent to x^2 sin(x).'''
    return (x**2) * numpy.sin(x)


assert numpy.isclose(integrate_numeric(xmin=0, xmax=4, N=50), 1.096591)

Notes:

• The ranges in the two for loops have been changed: we want the first for loop to go from xmin + h to xmin + (N-1)*h in steps of h (so N-1 total values), and the second for loop to go from xmin + h/2 to xmin + (N-1)*h + h/2 in steps of h (N total values).
• In the final computation, there's no need to apply f to msum and xsum: those values are already sums of f values. The only places we still need to evaluate f are at xmin and xmax. (Note: this was already fixed in an edit to the question.)
• The line h = (xmax-xmin)//N needs to be h = (xmax-xmin)/N. You just want a regular division here, not a floor division. This is likely the cause of the zeros you were getting originally: h would have been 0.