Trapezoidal rule in Python

Question

I'm trying to implement the trapezoidal rule in Python 2.7.2. I've written the following function:

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += h * f(a)
    for i in range(1, n):
        s += 2.0 * h * f(a + i*h)
    s += h * f(b)
    return s

However, f(lambda x:x**2, 5, 10, 100) returns 583.333 (it's supposed to return 291.667), so clearly there is something wrong with my script. I can't spot it though.

Answer 1

You are off by a factor of two. Indeed, the Trapezoidal Rule as taught in math class would use an increment like

s += h * (f(a + i*h) + f(a + (i-1)*h))/2.0

(f(a + i*h) + f(a + (i-1)*h))/2.0 is averaging the height of the function at two adjacent points on the grid.

Since every two adjacent trapezoids have a common edge, the formula above requires evaluating the function twice as often as necessary.

A more efficient implementation (closer to what you posted), would combine common terms from adjacent iterations of the for-loop:

f(a + i*h)/2.0 + f(a + i*h)/2.0 =  f(a + i*h) 

to arrive at:

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

print( trapezoidal(lambda x:x**2, 5, 10, 100))

which yields

291.66875