MNUM – Simpson’s Rule in Python

def simpson(f, a, b, n):
    "Approximate the definite integral of f from a to b by Simpson's rule."

    if n % 2 != 0:
        print "Ups: n must be even!"
        return -1
        
    h  = (float(b) - a)/n
    
    si = 0.0
    sp = 0.0
    
    for i in range(1, n, 2):
        xk = a + i*h
        si += f(xk)
    
    for i in range(2, n, 2):
        xk = a + i*h
        sp += f(xk)
        
        
    s = 2*sp + 4*si + f(a) + f(b)

    return (h/3)*s

def f(x):
    return x**4

ni = 50
nf = 1000000
n = ni
a = -20
b = 0
s = []
qc = []
ec = []
t = 0.0
div = 0.0
i = 0.0

while n < nf:
    t = simpson(f, a, b, n)
    s.append(t)
    n = n * 2

for i in range(0, len(s)-2):
    div = (s[i+2] - s[i+1])
    if div == 0:
        break
    
    t = (s[i+1] - s[i]) / div
    qc.append(t)
    t = div / 15
    ec.append(t)

for i in range(0, len(qc)):
    print "%.12f, %.12f, %.12f => qc=%.12f, e=%.12f" % (s[i], s[i+1], s[i+2], qc[i], ec[i])