<jupyter_start><jupyter_text>This is gonna be fun!!! BD  first we import the stuff we need!<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>now we define the function we want to integrate!<jupyter_code>def fun(x):
    return np.exp(-2*x)*np.cos(10*x)<jupyter_output><empty_output><jupyter_text>now we define its integral so we know what the right answer is! XD<jupyter_code>def fun_integral(x):
    return (1/52.)*np.exp(-2*x)*(5*np.sin(10*x)- np.cos(10*x))<jupyter_output><empty_output><jupyter_text>defining the single interval value of the trapezoid method! :-)<jupyter_code>def trapezoid_interval(f, x, h):
    return 0.5*h*(f(x+h)+ f(x))<jupyter_output><empty_output><jupyter_text>executing the trapezoid method: adding the single intervals values<jupyter_code>def trapezoid_method(f, a, b, N):
    # f = function to integrate
    # a = lower limit of integration
    # b = uper limit of integration
    # N = number of intervals to use
    
    #defining the number line interval to integrate over! :D
    x = np.linspace(a, b, N)
    h = x[1] - x[0]
    
    #defining which will hold the integral value as we sum
    Fint = 0.0
    
    #performing the intergal !!! :D
    for i in range(0, len(x)-1, 1):
        Fint += trapezoid_interval(f, x[i], h)
    
    #return the answer! :D
    return Fint<jupyter_output><empty_output><jupyter_text>Lets do the integral and compare it to the correct answer! :D<jupyter_code>print("Trapezoid method result")
print(trapezoid_method(fun, 0, np.pi, 3500))
print("Analytical solution")
Answer = fun_integral(np.pi)- fun_integral(0)
print(Answer)<jupyter_output><empty_output><jupyter_text>as we can see, about 3500 intervals is enough to reach the tolerance of 10^-6  Defining the single interval element of simpson's method! :-)<jupyter_code>def simpsons_interval(f, x, h):
    return h*( f(x) + 4*f(x+h) + f(x+2*h))/3.<jupyter_output><empty_output><jupyter_text>Executing simpsons method by summing the elements ! :D<jupyter_code>def simpsons_method(f, a, b, N):
    #f = func to integrate 
    #a = lower limit of integration 
    #b = upper limit of integration
    #N = number of intervals to use
    
    #the number of elements will be N -1 
    # so if N is odd we dont need to adjust anything
    
    #define the number line interval to integrat over
    x = np.linspace(a, b, N)
    h = x[1] - x[0]
    
    #define the variable to hold the integral value as we sum intervals! :D
    
    Fint = 0.0
    
    #executing simpsons method
    for i in range(0, len(x)-2, 2):
        Fint += simpsons_interval(f, x[i], h)
        
    #applying simpsons method to the last element sepratelt in case the
    #total number of intervals is even
    
    if((N%2)==0):
        Fint += simpsons_interval(f, x[-2], 0.5*h)
        
    return Fint<jupyter_output><empty_output><jupyter_text>Lets print our answer from simpsons method and compare it to the analytical solution! :-)<jupyter_code>print("Simpson's method result")
print(simpsons_method(fun, 0, np.pi, 150))
print("Analytical solution")
Answer = fun_integral(np.pi)- fun_integral(0)
print(Answer)<jupyter_output><empty_output><jupyter_text>as we can see within 150 intervals we got a result within the tolerance !  Defining the romberg integration interval value! :D<jupyter_code>def romberg_interval(f, a, b, i):
    
    #we need the diff b-a
    
    h = b-a
    
    #increment between new function evaluations!
    dh = h/2.**(i)
    
    # cofactor
    K = h/2.**(i+1)
    
    #function evaluations!
    M = 0.0
    for j in range(2**i):
        M += f(a + 0.5*dh + j*dh)
        
    return K*M<jupyter_output><empty_output><jupyter_text>executing the romberg method ! :)<jupyter_code>def romberg_integration(f, a, b, tol):
    
    #iteraton variable! :D
    i = 0
    
    #max number of iterations! 
    
    imax = 1000
    
    #defining an error estimate - set to a large value
    
    delta = 100.0*np.fabs(tol)
    
    #setting an array of integral answers
    I = np.zeros(imax, dtype = float)
    
    #0th romberg integration
    
    I[0] = 0.5*(b-a)**(f(a) + f(b))
    
    #iterate by 1
    
    i += 1
    
    while(delta>tol):
        
        I[i] = 0.5*I[i-1] + romberg_interval(f, a, b, i)
        
        #compute the new fractional error estimate
        
        delta = np.fabs((I[i]-I[i-1])/I[i])
        
        print(i, I[i], I[i-1], delta)
        
        if(delta>tol):
            
            i+=1
            
            #if we wave reached the ma xnumber of iterations
            
            if(i>imax):
                print("Max iterations reached. :(")
                raise StopIteration('Stopping iterations after', i)
                
    return I[i]<jupyter_output><empty_output><jupyter_text>Lets compute the romberg integral and compare it to the analytical solution!<jupyter_code>print("Romberg method result")
tolerance = 1.0e-6
RI = romberg_integration(fun, 0, np.pi, tolerance)
print(RI)
print("Analytical solution")
Answer = fun_integral(np.pi)- fun_integral(0)
print(Answer)<jupyter_output><empty_output>