Evaluate $\int_{0}^1 \int_{0}^1 \int_{0}^1 \frac{1}{x^2+y^2+z^2}dxdydz$ I am trying to evaluate the following integral:
$$I =\int_{0}^1 \int_{0}^1 \int_{0}^1 \frac{1}{x^2+y^2+z^2}dxdydz.$$
Here is my trial. By using spherical coordinate, let 
$x = r\sin(\theta)\cos(\phi), y = r\sin(\theta)\sin(\phi), z = r\cos(\phi).$
Then one can write $I$ as
$$ I = \int_{\Omega} \frac{1}{r^2} r^2\sin(\theta)drd\theta d\phi
= \int_{\Omega} \sin(\theta)drd\theta d\phi$$
where $\Omega$ is something I can't figure out...
Here $\Omega$ is
$$ \Omega = \{(r,\theta,\phi) : 0\le x,y,z \le 1\}.$$
Since $0\le z=r\cos(\phi) \le1$, we have
$$ 0 \le \arccos(1/r) \le \phi \le \pi/2. $$
Since $0\le y=r\sin(\theta)\cos(\phi) \le 1$, we have
$$ 0 \le \theta \le \arcsin(1/(r\cos(\phi)).$$
But I don't see how such classification help in evaluating $I$.
Any suggestions/comments/answers will be very appreciated.
Thanks in advance.
 A: I might have a better idea than parametrizing a cube in spherical coordinates (ugh), which is to further exploit symmetry:
$$ \mathcal{J}=\iiint_{(-1,1)^3}\frac{d\mu}{x^2+y^2+z^2}=8\iiint_{(0,1)^3}\frac{d\mu}{x^2+y^2+z^2}=48\iiint_{0\leq z\leq y\leq x\leq 1}\frac{d\mu}{x^2+y^2+z^2} $$
Now the substitution $z=yc, dz=y\,dc$ turns the RHS into
$$ 48\int_{0}^{1}\int_{0}^{x}\int_{0}^{1}\frac{y}{x^2+y^2+c^2 y^2}\,dc  \,dy\,dx $$
and the substitution $y=xb, dy=x\,db$ turns the last integral into
$$ 48\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\frac{x^2 b}{x^2+b^2 x^2+c^2 b^2 x^2}\,dc  \,db\,dx = 48\int_{0}^{1}\int_{0}^{1}\frac{ b}{1+b^2+c^2 b^2}\,dc  \,db$$
or
$$ 24\int_{0}^{1}\frac{\log(2+c^2)}{1+c^2}\,dc=24\int_{0}^{\pi/4}\log(2+\tan^2\theta)d\theta  $$
which is an integral in a single-variable, requiring the dilogarithms machinery to be computed in closed form. Numerically, $\mathcal{J}\approx\frac{353}{23}$. In explicit terms
$$ \mathcal{J}=24\left[\text{Im}\,\text{Li}_2((3-2\sqrt{2})i)-G+\frac{\pi}{4}\log(3+2\sqrt{2})\right]$$
where $G$ is Catalan's constant.
A: Here is another way: Let $\mathbf{F}$ be defined by
$$ \mathbf{F} = \left( \frac{x}{x^2+y^2+z^2}, \frac{y}{x^2+y^2+z^2}, \frac{z}{x^2+y^2+z^2} \right). $$
Then it is easy to check that $\operatorname{div}\mathbf{F} = \frac{1}{x^2+y^2+z^2}$. So if $I$ denotes the integral in question, then by the divergence theorem1) we get
$$ 8I
= \int_{[-1, 1]^3} \operatorname{div}\mathbf{F} \, dV
= \int_{\partial [-1, 1]^3} \mathbf{F} \cdot \mathbf{n} \, dA. $$
Again by symmetry, 
$$ 8I
= 6 \int_{\{1\}\times[-1,1]^2} \mathbf{F} \cdot \mathbf{e}_1 \, dydz
= 24 \int_{[0,1]^2} \frac{dydz}{1+y^2+z^2}.$$
Applying the polar coordinate change, we have
$$ \int_{[0,1]^2} \frac{dydz}{1+y^2+z^2}
= \int_{0}^{\frac{\pi}{4}} \int_{0}^{\sec\theta} \frac{2r}{1+r^2} \, drd\theta
= \int_{0}^{\frac{\pi}{4}} \log(1 + \sec^2\theta) \, d\theta,$$
which leads to the same integral as in Jack D'Aurizio's answer
$$ I = 3\int_{0}^{\frac{\pi}{4}} \log(2 + \tan^2\theta) \, d\theta. $$

1) It can be easily shown, by excising a small ball around $0$ and letting the radius zero, that the singularity of $\mathbf{F}$ at $0$ poses no issue.
