Double integral Simpsons rule in python where the limits are functions

Question

I have created a function where it can solve double integrals using the Simpsons rule but it only works when the limits of integration are constants

def double_integration(a,b,y1,y2,f,n):
    
    ### First, check that the number of strips is even.
    if n % 2 != 0:
        sys.exit("We require an even number of strips in Simpson's rule.")

    #Set up evaluation points in each direction
    xvals = np.linspace(a,b,n+1)
    yvals = np.linspace(y1,y2,n+1)

    hx = (b-a)/n #Mesh width in x direction
    hy = (y2-y1)/n #Mesh width in y direction

    #Set temporary weights in each coordinate direction
    weights = np.ones(n+1)
    weights[0] = 1/3
    weights[n] = 1/3
    
    for w in range (1,n):
        if w % 2 != 0:
            weights[w] = 4/3
        elif w % 2 == 0:
            weights[w] = 2/3

    #Scale with the mesh width in each direction to obtain actual weights
    x_weights = hx*weights
    y_weights = hy*weights

    integral = 0.0 #Initialise the approximation

    for i in range(n+1): #Loop over points in x-direction
        xi = xvals[i] 
        for j in range(n+1): #Loop over points in y-direction 
            yj = yvals[j]

            #Update the approximation
            integral += x_weights[i]*y_weights[j]*f(xi,yj)

    return integral

#First Test - Rectangular domain
f = lambda x,y: x**2*y
a = 1; b = 2
y1 = lambda x: np.ones(np.shape(x))
y2 = lambda x: 2*np.ones(np.shape(x))
n = 10
integral1 = double_integration(f,a,b,y1,y2,n)

#Second Test - Triangular domain
f = lambda x,y: np.ones(np.shape(x))
a = 0; b = 2
y1 = lambda x: -x
y2 = lambda x: x
n = 20
integral2 = double_integration(f,a,b,y1,y2,n)

but it doesn't pass these tests and will give the error: TypeError: unsupported operand type(s) for *: 'function' and 'float'. But i'm unsure on how to fix it

Answer 1

You do one little mistake y1 and y2 are not a number here and you deal like numbers.