The volume of the unit sphere\index{unit sphere}, $4\pi/3$, can be expressed as a triple integral in spherical polar co-ordinates with constant limits: $$ \int_0^{2\pi}\int_0^\pi\int_0^1 r^2\sin\theta\;\mathrm{d}r\mathrm{d}\theta\mathrm{d}\phi. $$
In [x]: 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)
Or in cartesian coordinates with limits as functions: $$ 8 \int_0^1 \int_0^{\sqrt{1-x^2}} \int_0^{\sqrt{1-x^2-y^2}} \; \mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x, $$ where the integral is in the positive octant of the three-dimensonal cartesian axes.
In [x]: A, _ = tplquad(lambda z, y, x: 1,
0, 1,
lambda x: 0, lambda x: np.sqrt(1 - x**2),
lambda x,y: 0, lambda x,y: np.sqrt(1 - x**2 - y**2))
In [x]: 8*A
Out[x]: 4.188790204786391