Numerically evaluate triple integral with non-constant limits

Question

Numerically evaluate triple integral with non-constant limits.

I want to numerically perform a triple integral of a two-variable functionf(x,y) where both x and y have integration limits [a+t0-t,b+t0-t] and the third integral goes over t with integration limits [t0,t0+a]. I have tried following several examples using tplquad or nquad from scipy.integrate but nothing seems to work... This is my previous attempt:

from scipy.integrate import quad, dblquad, tplquad, nquad

def f(x, y):
    # Define your function f(x, y) here
    return x + y  # Example function: x + y

def integrate_function(a, b, t0):
    # Define the limits of integration
    x_lower = lambda t: a + t0 - t
    x_upper = lambda t: b + t0 - t
    y_lower = lambda x, t: a + t0 - t
    y_upper = lambda x, t: b + t0 - t
    t_lower = t0
    t_upper = t0 + a

    # Perform the triple integration
    result, _ = tplquad(f, t_lower, t_upper, y_lower, y_upper, x_lower, x_upper)

    return result

# Example usage
a = 1
b = 2
c = 3
t0 = 0
result = integrate_function(a, b, t0)

Which yields the error:

TypeError: integrate_function.<locals>.<lambda>() missing 1 required positional argument: 't'

My current attempt is able to evaluate a double integral using nquad:

def f(x, y):
    return x*y

def bounds_y(t0, a):
    return [t0, t0+a]

def bounds_x(y,a,b,t0):
    return [a+t0-y, b+t0-y]

def integrate(a,b,t0):
    return nquad(f, [lambda y: bounds_x(y,a,b,t0), bounds_y(t0,a)])

integrate(2,5,0)

Unfortunately when I try to implement a third integral:

def f(x, y, t):
    return 1

def bounds_t(t0, a):
    return [t0, t0+a]

def bounds_x(y,a,b,t0):
    return [a+t0-y, b+t0-y]

def integrate(a,b,t0):
    return nquad(f, [lambda t: bounds_x(t,a,b,t0), lambda t: bounds_x(t,a,b,t0), bounds_t(t0,a)])

integrate(2,5,0)

it gives the error TypeError: integrate.<locals>.<lambda>() takes 1 positional argument but 2 were given

Answer 1

Apparently, this does the trick:

def f(x, y, t):
    return x*y

def bounds_t(t0, a):
    return [t0, t0+a]

def bounds_x(y,a,b,t0):
    return [a+t0-y, b+t0-y]

def integrate(a,b,t0):
    return nquad(f, [lambda t,x: bounds_x(t,a,b,t0), lambda t: bounds_x(t,a,b,t0), bounds_t(t0,a)])

integrate(2,5,0)

Answer 2

It does indeed require some careful reading of the documentation to use nquad and tplquad correctly. Here are solutions to your integration problem via both nquad and tplquad:

from scipy.integrate import tplquad, nquad


def f(x, y, t):
    return x * y


def integrate_via_tplquad(a, b, t0):
    t_lower = t0
    t_upper = t0 + a
    y_lower = lambda t: a + t0 - t
    y_upper = lambda t: b + t0 - t
    # Note order of parameters here: (t, y) rather than (y, t)
    x_lower = lambda t, y: a + t0 - t
    x_upper = lambda t, y: b + t0 - t
    result, _ = tplquad(f, t_lower, t_upper, y_lower, y_upper, x_lower, x_upper)
    return result


def integrate_via_nquad(a, b, t0):
    t_bounds = (t0, t0 + a)
    y_bounds = lambda t: (a + t0 - t, b + t0 - t)
    # Order of parameters (y, t) is opposite to that used in tplquad.
    x_bounds = lambda y, t: (a + t0 - t, b + t0 - t)
    result, _ = nquad(f, [x_bounds, y_bounds, t_bounds])
    return result


print("Result via tplquad: ", integrate_via_tplquad(2, 5, 0))
print("Result via nquad: ", integrate_via_nquad(2, 5, 0))

When I run the above on my machine, I get:

Result via tplquad:  118.50000000000001
Result via nquad:  118.50000000000001

Note that it's really important to get the order of parameters right, and that to make matters worse the conventions for tplquad and nquad are opposite. If you accidentally reverse y and t in the bounds for x in either solution, you get an incorrect answer:

Result via tplquad:  25.499999999999996
Result via nquad:  25.499999999999996

If you're using nquad, you could also pass a, b and t0 directly as arguments to the integrand and the bound functions instead of having to use the closures above. Here's what that looks like:

def integrand(x, y, t, a, b, t0):
    return x * y

def x_bounds(y, t, a, b, t0):
    return a + t0 - t, b + t0 - t

def y_bounds(t, a, b, t0):
    return a + t0 - t, b + t0 - t

def t_bounds(a, b, t0):
    return t0, t0 + a

print("Result via nquad (2): ", nquad(integrand, [x_bounds, y_bounds, t_bounds], (2, 5, 0)))

The result on my machine:

Result via nquad (2):  (118.50000000000001, 1.3156142841808107e-12)