Numerically evaluating $\iiint_{C} \frac{ {\rm d}x \, {\rm d}y \, {\rm d}z}{x^2 +y^2 + z^2}$

Question

Numerically evaluate $$\iiint_{C} \frac{ {\rm d}x \, {\rm d}y \, {\rm d}z}{x^2 +y^2 + z^2}$$ where $$C = [-1,1] \times [-1,1] \times [-1,1] = [-1,1]^3$$

---

I do not expect an analytic solution. I would be satisfied with a numeric value. But even this is not so easy.

How to deal with the singularity at the origin? One idea is to first integrate over the unit ball, which can be done analytically, then use monte carlo to integrate over the rest irregular region. But I am not sure whether this will converge fast enough.

Answer 1

We can use divergence theorem to say

$$\iiint_{[-1,1]^3} \frac{dV}{r^2} = \iint_{\partial [-1,1]^3} \frac{\hat{r}}{r}\cdot\vec{dS}$$

By symmetry we can choose to only integrate over the square in the $z=1$ plane and evaluate

$$\iint_{[-1,1]^2}\frac{6}{\sqrt{x^2+y^2+1}}\cos\theta\:dA = \iint_{[-1,1]^2}\frac{6}{x^2+y^2+1}\:dA$$

By even more symmetry we can reduce this to $8$ times an integral over a triangle in the first quadrant

$$\int_0^1\int_0^x \frac{48}{x^2+y^2+1}dydx$$

You can show this is equal to

$$\int_0^{\sinh^{-1}(1)}48\tan^{-1}(\tanh t)\:dt$$

by integrating directly or

$$\int_0^{\frac{\pi}{4}}24\log(\sec^2\theta+1)\:d\theta$$

by integrating in polar coordinates. Numerically these integrals both evaluate to $\approx 15.3482$

Answer 2

Here is a triple integral attempt in spherical coordinates. It does not lead to an analytic answer and I had to evaluate the integral with the help of WolframAlpha.

We consider the cube is made of $8$ unit cubes in $8$ octants. The unit cube in the first octant is bound by,

$x = 0, x = 1, y = 0, y = 1, z = 0, z= 1$. We again divide this into two parts by plane $x = y$ and integral in the sub-region bound by $y = 0, x = 1, x = y, z = 0, z = 1$ can be evaluated as below using spherical coordinates -

$x = \rho \cos \theta \sin \phi, y = \rho \sin \theta \sin \phi, z = \rho \cos \phi$. $x^2+y^2+z^2 = \rho^2$.

$0 \leq \theta \leq \frac{\pi}{4}$.

i) when we are bound by plane $z = 1$

$z = 1 \implies \rho = \sec \phi$

To find the upper bound of $\phi$, let's look at the intersection of plane $z = 1, x = 1$ which is given by,

$\rho = \sec\phi = \sec\theta \csc\phi \implies \phi = \arctan(\sec\theta)$

So, $0 \leq \phi \leq \arctan(\sec\theta)$

$I_1 = \displaystyle \int_{0}^{\pi/4} \int_0^{\arctan(\sec\theta)} \int_0^{\sec\phi} \sin\phi \ d\rho \ d\phi \ d\theta \ \approx 0.319755$

ii) when we are bound by plane $x = 1$

$x = 1 \implies \rho = \sec \theta \csc \phi$

$\arctan(\sec\theta) \leq \phi \leq \frac{\pi}{2}$

$I_2 = \displaystyle \int_0^{\pi/4} \int_{\arctan(\sec\theta)}^{\pi/2} \int_0^{\sec\theta \csc\phi} \sin\phi \ d\rho \ d\phi \ d\theta \approx 0.639510$

So considering symmetry, the integral over the entire cubic region evaluates to,

$I = 16(I_1 + I_2) \approx 15.348$.