I have a segment of code meant to be used in a script to find the error in a numerical integration using the trapezoidal rule by fitting the error to a polynomial. This segment of code is throwing a float divide by zero error and I don't see why or how to fix it.
Could someone help lead me to the answer?
def trap(f,a,b,dx,exact):
N = int(numpy.round(float(b-a)/dx))
w=(b-a)/N
sum = f(a)/2.0 + f(b)/2.0
for i in range(1,N):
sum += f(a+i*w)
area = sum * w
errorf = exact-area
# If the error crosses 0, a polynomial approximation
# to the absolute value will go crazy.
return errorf
This alternate method throws the same error
# alternate way to handle dx not a divisor of b-a
def alt_trap(f,a,b,dx,exact):
N = int(numpy.floor(float(b-a)/dx))
sum = f(a)/2.0 + f(a+N*dx)/2.0
for i in range(1,N):
sum+= f(a+i*dx)
area = sum*dx
# now add one trapezoid between a+Ndx and b
area += 1/2*(b-(a+N*dx))*(f(b)+f(a+N*dx))
errorf = exact-area
return errorf
if dx == 0:dx = 1e-8... might do what you want ...basically you cant divide by zero so you must define what you expect to happen in that case ... - Joran Beasley