Why is my Simpson's modelling not giving the correct answer?

Question

I've been trying to do some Simpson's modelling to approximate integrals and although my code runs fine, it doesn't seem to give the correct answer. I'm new-ish to coding so I understand it's probably a dumb mistake but I can't figure out what's wrong. f(x) has already been defined and seems to work, it's just this integral approximation that's a bit dodgey.

def integrate_numeric(xmin, xmax, N):
    ''' 
    Numerical integral of f from xmin to xmax using Simpson's rule with 
        N panels.
    '''
    h = (xmax - xmin)/2*N
   
    # Odd sum
    oddtotal = 0
    for i in range(1, N, 2):
        oddtotal += f(xmin + (i*h))
    f_odd = oddtotal
    
    # Even sum
    eventotal = 0
    for i in range(2, N, 2):
        eventotal += f(xmin + (i*h))
    f_even = eventotal
    
    return (h/3)*(f(xmin)+(4*f_odd)+(2*f_even)+f(xmax))

Which answer does it give and which answer would be correct? — mkrieger1, Oct 23 '22 at 15:36
For which f, xmin, xmax and N is this? Can you please show a [mre]? — mkrieger1, Oct 23 '22 at 15:39
Possibly not, it's my first time trying Simpson's modeeling. I'm reasonably confident in h though and I intended to do it so that f_odd is f(xmin+(h * odd integers)) up to N and similarly for f_even (but this time with even integers) — Adam P, Oct 23 '22 at 15:40
`def f(x): return numpy.cos(2*x)*x**2` and between 0 and 4 for N=1 is what I've been testing — Adam P, Oct 23 '22 at 15:43

Tom McLean · Answer 1 · 2022-10-23T15:56:15.760

The ranges for odd numbers need to be from x=1 to x=N-1
The ranges for even numbers need to be from x=2 to x=N-2
h is needs to be (xmax - xmin)/N

Here is the Simpsons rule, with those corrections:

def f(x):
    return x**2


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

    odd_sum = sum(f(xmin + i*h) for i in range(1, N, 2))
    even_sum = sum(f(xmin + i*h) for i in range(2, N-1, 2))

    return h/3 * (f(xmin) + 4*odd_sum + 2*even_sum + f(xmax))


print(integrate_numeric(0, 10, 50)) -> 333.3333333333333

Alternatively, just use scipy's simpsons rule

Why is my Simpson's modelling not giving the correct answer?

1 Answers1