1
votes

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
1
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
Oh wow, that was obvious now that you point it out. Thank you! - Errata

1 Answers

1
votes

dx is the only possible parameter that can be 0 and throw ZeroDivisionError. You can catch the exception and decide what to do or fix the input - depends on your logic.

In order to provide a default value just add:

dx = dx or 0.0000001 # 0 is false and python is awesome to support this syntax

if you want to Try and catch the exception (Better to ask for forgive then permission)

import sys # At the top
...
try:
    N = int(numpy.floor(float(b-a)/dx))
except ZeroDivisionError, e:
    print a, b, dx 
    N = sys.maxint