Finding the volume of a torus

Consider a torus of average radius $R$ and cross sectional radius $r$. The volume of this shape may be evaluated analytically in cartesian coordinates as a volume of revolution: $$ V = 2\int_{R-r}^{R+r} 2\pi x z\;\mathrm{d}x, \quad \mathrm{where}\;z = \sqrt{r^2 - (x-R)^2}. $$ The centre of the torus is at the origin and the $z$ axis is taken to be its symmetry axis.

The integral is tedious but yields to standard methods: $V = 2\pi^2Rr^2$. Here we take a numerical approach with the values $R=4$, $r=1$:

In [x]: import numpy as np
In [x]: from scipy.integrate import quad
In [x]: R, r = 4, 1
In [x]: f = lambda x, R, r: x * np.sqrt(r**2 - (x-R)**2)
In [x]: V, _ = quad(f, R-r, R+r, args=(R, r))
In [x]: V *= 4 * np.pi
In [x]: Vexact = 2 * np.pi**2 * R * r**2
In [x]: print('V = {} (exact: {})'.format(V, Vexact))
Out[x]: V = 78.95683520871499 (exact: 78.95683520871486)