from numpy import *
    
def trapezoid(xi, f, F=[]):
    """
    Use trapezoid quadrature based on the subdivision xi
    Input:
    xi: location of nodes, defining the subdivision of the interval
    (including the endpoints)
    f:  function to integrate
    F:  antiderivative (for computing errors)
    
    
    Output:
    
    I = the numerical approximation to the integral on each subinterval
    E = error (absolute value) on each subinterval (only if F is known)
    
    """
    m = len(xi)-1
    E = zeros(m)
    I = zeros(m)
    for i in range(m):
        hi = xi[i+1]-xi[i]
        Ii_trapezoid = hi/2*(f(xi[i])+f(xi[i+1]))
        if F:
            Ii_exact = F(xi[i+1])-F(xi[i]) # in reality we of course do not know this
            E[i] = abs(Ii_trapezoid-Ii_exact)
        I[i] = Ii_trapezoid
    return I, E

## Test Trapezoid
if __name__ == "__main__":
    N = 1000
    a = 0
    b = 1
    x = linspace(a,b,N)
    f = lambda x: sqrt(x)
    F = lambda x: x*sqrt(x)*2./3.
    I, E = trapezoid(x, f, F)
    Iq= sum(I)
    print("N = %6d, numerical = %e, error = %e" % (N,Iq,F(b)-F(a)-Iq))