from numpy import *

def trapezeEasy(f,a,b,n):
    length = (b-a)/n
    lst = arange(a,b+length,length)
    U = zeros(len(lst))
    U[0],U[-1] = f(lst[0]),f(lst[-1])
    U[1:len(U)-1] = 2*f(lst[1:len(lst)-1])
    I = (length/2)*sum(U)
    return I

def trapezeFun(f,a,b,n,nmax,tol):
    errorEst = tol+10
    length = (b-a)/n
    I = trapezeEasy(f,a,b,n)
    while errorEst > tol and 2*n<=nmax:
        I2 = I
        old_length = length
        n = n*2
        length = length/2
        I = I2/old_length*length + (length)*sum(f(arange(a+length,b,old_length)))
        errorEst = abs(I2-I)
    return I,n,errorEst



from matplotlib import pyplot as plt
from numpy import *
def u(x):
    return cos(x)
a = 0
b = pi/2
n = 10
I = trapezeEasy(u,a,b,n)
errorExact = abs(1.0 - I)
print(" ======= Integral of sinus between 0 and pi/2 = %21.14e " % I)
print(" True error = %14.7e" % errorExact)
print(" Number of intervals = %d" % n)
I,n,errorEst = trapezeFun(u,a,b,n,200000,1e-12)
errorExact = abs(1.0 - I)
print(" ======= Integral of sinus between 0 and pi/2 = %21.14e " % I)
print(" True error = %14.7e" % errorExact)
print(" Est. error = %14.7e" % errorEst)
print(" Number of intervals = %d" % n)
plt.figure("Discovering the trapezoids integration rule :-)")
x = [a,b]
plt.plot(x,u(x),'.k',markersize=5)
x = linspace(a,b,200)
plt.fill(append(x,[0]),append(u(x),[0]),'xkcd:sky blue')
x = linspace(-pi/2,pi,300)
plt.plot(x,u(x),'-k')
plt.title('Integrating sinus between 0 and pi/2')
plt.gca().axhline(y=0,color='k',linewidth=1.0)
plt.show()
  

