Doubt in solution for evaluating $\int_0^1\int_0^1\int_0^1(1+u^2+v^2+w^2)^{-2}du~dv~dw$. I have two doubts in the answer for evaluating the following integral:
$$\int_0^1\int_0^1\int_0^1(1+u^2+v^2+w^2)^{-2}du~dv~dw$$

Solution: call this integral as $I$. By symmetry we may compute it over the domain $\{(u,v,w):0\leq v\leq u\leq 1\}$ and then double the result.
Substitute $u=r\cos(\theta)$, $v=r\sin(\theta)$, $w=\tan(\phi)$. Now the limits of integration become $0\leq \theta,\phi\leq \pi/4$ and $0\leq r\leq  \sec(\theta)$. Then finally we have:
$$I=2\int_0^{\pi/4}\int_0^{\pi/4}\int_0^{\sec(\theta)}\frac{r\sec^2(\phi)}{(r^2+\sec^2(\phi))^2}dr~d\theta ~d\phi.$$

My doubts are:
1.) I don't get how by symmetry are we computing $I$ on $\{(u,v,w):0\leq v\leq u\leq 1\}$ .
2.) I don't understand why we substitute $u=r\cos(\theta)$, $v=r\sin(\theta)$, $w=\tan(\phi)$?
Kindly help me with above doubts. Thanks in advance.
 A: This integral is a real beaut. I'll try to answer your questions, explain the substitutions, and then finally reduce the integral to a (very difficult) known definite integral.


*

*"I don't get how by symmetry are we computing $I$ on $\{(u,v,w) : 0\le v\le u \le 1\}$."


A picture or two is worth a thousand words: setting $w=0$ (this doesn't change things), compare the whole and half ranges.





*"I don't understand why we substitute $u=r\cos(\theta),$ $v=r\sin(\theta)$, $w=\tan(\phi)$."


Basically, we want to exploit the radial symmetry of the integrand. After using the symmetry trick, the region of integration is a triangular prism. You have correctly computed the new limits of integration and the Jacobian determinant, so we're ready to go.

Start with the new integral
$$
I =  \int _{0}^{\pi/4}\int_{0}^{\pi/4} \int _0^{\sec(\theta)} \frac{2r \sec^2(\phi)}{(r^2+\sec^2(\phi))^2}\,\mathrm{d}r \mathrm{d}\theta \mathrm{d}\phi
$$Integrating with respect to $r$ doesn't pose too much trouble:
$$
I = \left. \int _{0}^{\pi/4}\int_{0}^{\pi/4} \frac{-\sec ^2(\phi )}{r^2+\sec ^2(\phi )}\right|_{0}^{\sec(\theta)}\,\mathrm{d}\theta \mathrm{d}\phi
$$
$$
 = \int _{0}^{\pi/4}\int_{0}^{\pi/4} \frac{1}{1+\cos ^2(\theta ) \sec ^2(\phi )}\,\mathrm{d}\theta \mathrm{d}\phi
$$This is basically an arctan integral, we just have to scale it appropriately. Using $\displaystyle{\frac{\mathrm{d}}{\mathrm{d}x} \frac{\arctan\left(\frac{\tan (x)}{\sqrt{1+M^2}}\right)}{\sqrt{1+M^2}}=\frac{1}{1+M^2\cos^2(x)}}$, we have
$$
I = \left. \int _{0}^{\pi/4}\frac{\arctan\left(\frac{\tan (\theta )}{\sqrt{1+\sec ^2(\phi )}}\right)}{\sqrt{1+\sec ^2(\phi
   )}} \right|_{0}^{\pi/4}\,\mathrm{d}\phi
$$
$$
=  \int _{0}^{\pi/4}\frac{\arctan\left(\frac{1}{\sqrt{1+\sec ^2(\phi )}}\right)}{\sqrt{1+\sec ^2(\phi
   )}} \,\mathrm{d}\phi
$$Let $\sec(\phi)=\sqrt{y^2+1}$. Then $\sec(\phi)\tan(\phi)\mathrm{d}\phi = \frac{y}{\sqrt{y^2+1}}\mathrm{d}y$, or $\mathrm{d}\phi=\frac{1}{1+y^2} \mathrm{d}y$:
$$
=  \int _{0}^{1}\frac{\arctan\left(\frac{1}{\sqrt{2+y^2}}\right)}{\sqrt{2+y^2}(1+y^2)} \,\mathrm{d}y
$$Use the reciprocal relation of arctan:
$$
=  \int _{0}^{1}\frac{\pi/2-\arctan\left({\sqrt{2+y^2}}\right)}{\sqrt{2+y^2}(1+y^2)} \,\mathrm{d}y
$$
$$
=  \underbrace{\frac{\pi}{2}\int _{0}^{1}\frac{1}{\sqrt{2+y^2}(1+y^2)} \,\mathrm{d}y}_{I_1}-\underbrace{\int _{0}^{1}\frac{\arctan\left({\sqrt{2+y^2}}\right)}{\sqrt{2+y^2}(1+y^2)} \,\mathrm{d}y}_{I_2}
$$For the first integral, we have
$$
I_1=\frac{\pi}{2}\int _{0}^{1}\frac{(2+y^2)}{(2+y^2)^{3/2}(1+y^2)} \,\mathrm{d}y
$$
$$
=\frac{\pi}{4}\int _{0}^{1}\frac{(2+y^2)}{(1+y^2)}\cdot \frac{2}{(2+y^2)^{3/2}} \,\mathrm{d}y
$$
$$
=\frac{\pi}{2}\int _{0}^{1}\frac{1}{2}\cdot\frac{1}{1+\left(\frac{y}{\sqrt{2+y^2}}\right)^2}\cdot \frac{2}{(2+y^2)^{3/2}} \,\mathrm{d}y
$$
$$
=\left.\frac{\pi}{2}\arctan\left(\frac{y}{\sqrt{2+y^2}}\right)\right|_0^1 = \frac{\pi^2}{12}
$$$I_2$ is much more difficult to evaluate; fortunately, it is known as Ahmed's Integral and a very nice calculation giving the value of $\displaystyle{\frac{5\pi^2}{96}}$ can be found here. Putting it together, we have $\displaystyle{I = \frac{\pi^2}{32}}$.
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
$\ds{\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
\pars{1 + u^{2} + v^{2} + w^{2}}^{-2}\,\,\dd u\,\dd v\,\dd w}:\
{\Large ?}}$.

With $\ds{\large\vec{r} \equiv u\,\hat{u} + v\,\hat{v} +w\,\hat{w}}$:
\begin{align}
&\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
\pars{1 + u^{2} + v^{2} + w^{2}}^{-2}\,\,\dd u\,\dd v\,\dd w} =
\iiint_{\large\pars{0,1}^{3}}{\dd^{3}\vec{r} \over
\pars{1 + r^{2}}^{2}}
\\[5mm] = &\
\iiint_{\large\pars{0,1}^{3}}\bracks{%
-\,{1 \over 4\pi}\,\nabla^{2}\Phi\pars{r}}{\dd^{3}\vec{r}}
\\[2mm] &\ \mbox{where}\
\Phi\pars{r} = 4\pi\int_{0}^{\infty}
{1 \over \pars{1 + r\,'^{2}}^{2}}\,{r\, '^{2}\,\dd r\, ' \over \max\braces{r,r'}} = 2\pi\,{\arctan\pars{r} \over r}
\end{align}
It turns out -via
Gauss-Ostrogradsky's Divergence Theorem- that
\begin{align}
&\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
\pars{1 + u^{2} + v^{2} + w^{2}}^{-2}\,\,\dd u\,\dd v\,\dd w} 
\\[5mm] = &
-\,{1 \over 4\pi}\iint_{S}
\nabla\Phi\pars{r}\cdot\dd\vec{S}\quad
\pars{\begin{array}{l}
\ds{S}\ \mbox{is the}\ \ds{\pars{0,1}^{3}}\
\mbox{box surface}
\end{array}}
\end{align}
By symmetry considerations:

*

*The contribution, to the surface integration, from the sides that intersects the origin of coordinates $\underline{vanish\ out}$.

*Each of the $\underline{\color{red}{three}}$ remaining sides yields the same contribution to the surface integration.

Then,
\begin{align}
&\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
\pars{1 + u^{2} + v^{2} + w^{2}}^{-2}\,\,\dd u\,\dd v\,\dd w} 
\\[5mm] = &\
\color{red}{3}\braces{-\,{1 \over 4\pi}\iint_{\pars{0,1}^{2}}
\bracks{\partiald{\Phi\pars{r}}{z}}_{\ z\ = 1}\dd x\,\dd y}
\\[5mm] = &\
-\,{3 \over 2}\iint_{\pars{0,1}^{2}}
\braces{\partiald{}{r}\bracks{{\arctan\pars{r} \over r}}
\, {1 \over r}}_{\ z\ =\ 1}\dd x\,\dd y
\\[5mm] = &\
3\int_{1}^{\infty}\bracks{\int_{0}^{1}\int_{0}^{1}
{\dd x\,\dd y \over \pars{x^{2} + y^{2} + t^{2} + 1}^{2}}}\dd t
\\[5mm] = &\
3\int_{1}^{\infty}{\mrm{arccot}\pars{\root{2 + t^{2}}} \over
\pars{1 + t^{2}}\root{2 + t^{2}}}\,\dd t
\\[5mm] = &\
{3\pi \over 2}\
\underbrace{\int_{1}^{\infty}
{\dd t \over \pars{1 + t^{2}}\root{2 + t^{2}}}}
_{\ds{\pi \over 12}}\ -\
3\
\underbrace{\int_{1}^{\infty}{\arctan\pars{\root{2 + t^{2}}} \over
\pars{1 + t^{2}}\root{2 + t^{2}}}\,\dd t}
_{\ds{\pi^{2} \over 32}}
\\[5mm] = &\
\bbx{\large{\pi^{2} \over 32}} \approx 0.3084 \\ &
\end{align}
