Simpson's Rule

from math import log, cos
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
def f(x):
    return (log(x+1)**(1+cos(x)))

print simpson(f, 3, 4, 2)
exit()

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


# Calcular S para diferentes valores de n [n(i+1) = ni * 2]
# para poder ter S (h), S' (h'=h/2), S'' (h''=h/4)
while n < nf:
    t = simpson(f, a, b, n)
    s.append(t)
    n = n * 2

# Calcular quociente (S'-S)/(S''-S') [que deve ser aprox igual a 16]
# e erro em que e = (S''-S')/15
for i in range(0, len(s)-2):
    div = (s[i+2] - s[i+1]) # S''-S'
    if div == 0:
        break
    
    t = (s[i+1] - s[i]) / div # (S'-S)/div com div = (S''-S')<=>(S'-S)/(S''-S')
    qc.append(t)
    t = div / 15
    ec.append(t)

#qc ~= 16 implica validade => podemos usar os valores do erro: e => qualidade
for i in range(0, len(qc)):
    print "S=%.12f, S'=%.12f, S''=%.12f\n\t=> qc=%.12f, e=%.12f\n" \
            % (s[i], s[i+1], s[i+2], qc[i], ec[i])