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

It might help to compare the lambda functions provided to [the requirements](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.tplquad.html). I think you've got them in the wrong order. — tadman, May 30 '23 at 18:37
@tadman I think so... I wouldn't know what is wrong with them... — TheHunter, May 30 '23 at 18:39
It might help to match them up to the arguments using labels, like `gfun=...`, or even better, specify those inline, like `gfun=lambda ...` — tadman, May 30 '23 at 18:39

score 0 · Answer 1 · answered May 30 '23 at 19:38

0

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)

answered May 30 '23 at 19:38

TheHunter

111
2

For your first bound (which are the bounds for `x`), I think you want `lambda y, t: ...` instead of `lambda t, x: ...`. Using `x` in place of `y` for the name is harmless, of course, but using `t, y` in place of `y, t` will end up doing the wrong integration. – Mark Dickinson May 31 '23 at 18:26
Will it do a wrongful integration because it integrates ``t`` instead of ``y`` or because it integrates ``x`` and ``y`` is the wrong order? – TheHunter Jun 01 '23 at 12:45
SciPy will pass the value of `y` as the first argument to that bound, so your bounds for `x` will effectively become `a + t0 - y` and `b + t0 - y` instead of `a + t0 - t` and `b + t0 - t`. – Mark Dickinson Jun 01 '23 at 17:17

Mark Dickinson · Answer 2 · 2023-05-31T18:34:29.067

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)

Numerically evaluate triple integral with non-constant limits

2 Answers2