In [1]:

    

import numpy as np


    

In [ ]:

    

#write a function that computes a 2d integral numerically
#needs to take a function, the x limit, the y limit, and assigns the number of points for x and y limit


    

In [4]:

    

#simple function designed to test
def f1(x,y):
    return (x**2)*y


    

In [17]:

    

def integral(function, x1, x2, y1, y2, xpoints = 1000., ypoints = 1000.):
    x_box = np.linspace(x1, x2, xpoints)
    y_box = np.linspace(y1, y2, ypoints)
    delta_x = (x2 - x1)/xpoints
    delta_y = (y2 - y1)/ypoints
    #the volume of the boxes
    volume = delta_x*delta_y
    total_sum = []
    for i,x in enumerate(x_box):
        a_sum = np.sum(function(x, y_box))
        total_sum.append(a_sum)
    
    function_sum = np.sum(total_sum)
    return volume*function_sum


    

In [20]:

    

integral(f1, 0, 1, 0, 1)


    

    
Out[20]:





0.16675008341675007

In [21]:

    

#percentage difference between theoretical and actual
((integral(f1, 0, 1, 0, 1) - 1./6.)/(1./6.))*100


    

    
Out[21]:





0.050050050050048922

In [ ]: