How to define limits of integration for a triple integral in SciPy using tplquad?

Question

I cannot understand how to approach defining the limits of integration for a triple integral in Scipy with the tplquad function. I am using the book Learning Scientific Programming With Python By Christian Hill and trying to solve problem E8.14:

The volume of the unit sphere, 4π/3, can be expressed as a triple integral in spherical polar co-ordinates with constant limits:

My attempt to evaluate this integral:

import scipy.integrate as integrate
#define the integrand
f = lambda r, theta, phi: r**2 * np.sin(theta)
#integrate with tplquad
V_unit_sphere, _ = integrate.tplquad(f, 0,1,0,np.pi,0,2*np.pi)
V_unit_sphere

My output was 165... clearly not equal to 4pi/3.

Hill's Solution:

from scipy.integrate import tplquad
In [x]: tplquad(lambda phi, theta, r: r**2 * np.sin(theta),
                0, 1,
                lambda theta: 0, lambda theta: np.pi,
                lambda theta, phi: 0, lambda theta, phi: 2*np.pi)
Out[x]: (4.18879020478639, 4.650491330678174e-14)

I do not understand why Hill defined the integrand in the order of phi, theta, r when the integral is taken the other way around, and I do not understand why even though all the limits of integration are constants, they must be defined as depending on some variable. I also do not understand why the limits of integration for theta are dependent on theta, and why the limits for phi are dependent upon both phi and theta.

jared · Answer 1 · 2023-06-12T18:43:50.503

Although the solution in the textbook works, I don't think it is the correct way to solve this problem.

According to the documentation of scipy.integrate.tplquad, the integrated function should be defined as f(z,y,x), and it will be integrated in that order (z, y, then x). In your example, the order of integration is r, theta, phi, so the function should be f(r, theta, phi), as you defined it.

The limits are defined in the opposite order, so the first two limits are the limits of integration for x (phi in your case), followed by y (theta), and finally z (r). Because the order of integration is z, y, then x, we can have the limits for z be functions of x and y, and the limits for y can be functions of x. In this case, the limits are all constants, so we can just provide constants (Hill used functions that returned constants, which is more general but has the same result).

That said, I believe the code should be:

import numpy as np
from scipy.integrate import tplquad

volume = tplquad(lambda r, theta, phi: r**2*np.sin(theta),
                 0, 2*np.pi,            # limits of integration for phi
                 0, np.pi,              # limits of integration for theta
                 0, 1)                  # limits of integration for r

print(volume)           # (4.18879020478639, 4.650491330678174e-14)

Moral of the story: Just because it works and was written in a book doesn't mean it is the right way to do it. Always test things yourself and read the documentation.

How to define limits of integration for a triple integral in SciPy using tplquad?

1 Answers1