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))