[Python] Computing an iterated integral by iteratively applying the 1D trapezoidal rule

I have a Python function, called trap_1D, which computes the (approximate) integral of a univariate function f over an interval [a,b] using the composite trapezoidal rule. Now I want to evaluate an iterated integral of the form \int_{a}^{b} \int_{c}^{d} f(x,y) dy dx by iteratively calling trap_1D: Below is trap_1D:

import numpy as np

def trap_1D(func, a, b, N, *args):

#Computes the integral of a univariate function over [a, b] using the composite trapezoidal rule.

h = (b - a) / N #The grid spacing
x = np.linspace(a, b, N + 1) #The sub-division points: x_0, x_1,...,x_N.
y = func(x, *args)

return (np.sum(y) - (y[0] + y[-1]) / 2) * h


I've tested trap_1D on a few functions and it works. Next, I constructed a function I(x) := \int_{a}^{b} f(x,y) dy.

def create_I(f, a, b, N):
def I(x):
return trap_1D(lambda y: f(x, y), a, b, N) #integrate the function y --> f(x,y)

return I


I've tested create_I and it works. Now to test out evaluating iterated integrals I'm using f(x,y) = sin(x) * sin(y) as my test function, which I want to integrate over the rectangle [0,pi/2] x [0,pi/2]:

def f(x,y):
return np.sin(x)*np.cos(y)

N = 100
I = create_I(f, 0.0, np.pi/2, N) #The function I() is correct
I2 = trap_1D(I, 0.0, np.pi/2, N) #This gives an error


The last line gives the following error:

---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_15880\410820249.py in <module>
26 N = 100
27 I = create_I(f, 0.0, np.pi/2, 100)
---> 28 trap_1D(I, 0.0, np.pi/2, N)

~\AppData\Local\Temp\ipykernel_15880\410820249.py in trap_1D(func, a, b, N, *args)
9 y = func(x, *args)
10
---> 11 return (np.sum(y) - (y[0] + y[-1]) / 2) * h
12
13

IndexError: invalid index to scalar variable.


I don't understand this error. What is causing the problem here?

Continue reading...