Numerical Analysis

import numpy as np

def midpoint(f, x0, xn, n):
    X=[0.0]*(4*n)
    Y=[0.0]*(4*n)

    dx = float(xn-x0) / float(n)
    VLines = [x0+i*dx for i in xrange(n)]
    for i in xrange(n):
        X[4*i] = X[4*i+1] = float(x0+i*dx)
        X[4*i+2] = X[4*i+3] = x0+(i+1)*dx

        Y[4*i] = Y[4*i+3] = 0.0
        Y[4*i+1] = Y[4*i+2] = f( (X[4*i]+X[4*i+3]) / 2.0)
    return X,Y,VLines

def trapezoid(f, x0, xn, n):
    X=[0.0]*(4*n)
    Y=[0.0]*(4*n)

    dx = float(xn-x0) / float(n)
    VLines = [x0+i*dx for i in xrange(n)]
    for i in xrange(n):
        X[4*i] = X[4*i+1] = float(x0+i*dx)
        X[4*i+2] = X[4*i+3] = x0+(i+1)*dx

        Y[4*i] = Y[4*i+3] = 0.0
        Y[4*i+1] = f(X[4*i])
        Y[4*i+2] = f(X[4*i+3])
    return X,Y,VLines

def simpson(f, x0, xn, n):
    X,Ym,V = midpoint(f,x0,xn,n)
    Xt,Yt,Vt = trapezoid(f,x0,xn,n)
    Y = [(2.0*M+T)/3.0 for M,T in zip(Ym,Yt)]
    return X,Y,V