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)
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:
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).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.)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
.