This page contains an interesting identity$$\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx=\frac{1}{15} \left(\sqrt{2}+2+5 \log \left(\sqrt{2}+1\right)\right)$$
Which calculates the average distance between two random points in the unit square. I'm wondering if there's an elementary solution? Any help will be appreciated.
Evaluate $\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx$
definite integralsintegration
Related Solutions
One can exploit a symmetry and split the integral in half along the line $y=x$ then use the following modified polar coordinates:
$$x = s^{\frac{2}{3}}\cos\theta \hspace{10 pt} y = s^{\frac{2}{3}}\sin\theta$$
$$\implies 2\cdot \frac{2}{3}\int_0^{\frac{\pi}{4}} \int_0^{\sec^{\frac{3}{2}}\theta} \frac{\log\left(s^{\frac{4}{3}}\right)}{\sqrt{\cos\theta+\sin\theta}}dsd\theta = \frac{16}{9}\int_0^{\frac{\pi}{4}} \frac{\sec^{\frac{3}{2}}\theta}{\sqrt{\cos\theta+\sin\theta}}\left(\frac{3}{2}\log(\sec\theta)-1\right)d\theta$$
then let $x=\tan\theta$
$$ \implies \frac{4}{3}\int_0^1 \frac{\log(1+x^2)}{\sqrt{1+x}}dx - \frac{16}{9}\int_0^1 \frac{1}{\sqrt{1+x}}dx$$
The integral on the right evaluates to $\frac{32}{9}(\sqrt{2}-1) \equiv \frac{32}{9}a$. The integral on the left becomes
$$ = \frac{8}{3}\sqrt{2}\log 2 - \frac{16}{3}\int_0^1 \frac{x\sqrt{1+x}}{1+x^2}dx$$
We have one and a half of the three terms the user posted, but this last integral is tricky and isn't yielding to a variety of methods. I will try to finish later, but in the meantime if anyone has any clever suggestions for this last integral I will be happy try it.
In order to prove that $S=-\frac{\pi^2}{4}$ , where \begin{equation*} S=\int_{0}^{1}\dfrac{\sqrt{x-x^3}\log(x)}{x(x^2+1)}\, dx +\int_{0}^{\frac{\pi}{2}}\dfrac{x\sqrt{\cos(x)}}{\sin(x)}\, dx, \end{equation*} we will use Cauchy's integral theorem. The integrand will be \begin{equation*} f(z)= \dfrac{2i\log(z)\sqrt{\frac{1+z^2}{2}}}{\sqrt{z}(1-z^2)} \end{equation*} where \begin{equation*} \log(z) = \ln|z| +i\arg(z) \quad \mbox{ with } -\pi<\arg[z)<\pi \end{equation*} and \begin{equation*} \sqrt{z} =e^{\frac{1}{2}\log(z)}. \end{equation*} Let $C$ be the contour of the first quadrant of the unit circle. Then \begin{equation*} \int_{C}f(z)\, dz = 0.\tag{1} \end{equation*} Some additional justifications are probably needed.
Split $C=C_1+C_2+C_3$ where $C_1$ is the path along the real axis from $0$ to $1$, $C_2$ is the arc from $1$ to $i$ and $C_3$ is the path along the imaginary axis from $i$ to $0$.
The singularity in $1$ is removable.
Close to $0$ we have to slightly modify $C_2$ and $C_3$. At the point $ir$ we leave $C_3$ and follow a small arc $C_r $ with radius $r$ to $r$ on $C_1$. According to the ML inequality \begin{equation*} \int_{C_r}f(z)\, dz \to 0 \quad \mbox{ as } r\to 0. \end{equation*}
The point $i$ can be treated analogously.
We will now study (1). \begin{equation*} \int_{C_1}f(z)\, dz = 2i\int_{0}^{1}\dfrac{\ln(x)\sqrt{\frac{1+x^2}{2}}}{\sqrt{x}(1-x^2)}\, dx \end{equation*} with real part $= 0$.
The arc $C_2$ can be described as $z=e^{it}, \quad 0 <t< \frac{\pi}{2}$. Then \begin{gather*} \int_{C_2}f(z)\, dz = \int_{0}^{\frac{\pi}{2}}\dfrac{2i\log(e^{it})\sqrt{\frac{1+e^{i2t}}{2}}}{\sqrt{e^{it}}(1-e^{i2t})}ie^{it}\, dt =\\[2ex]\int_{0}^{\frac{\pi}{2}}\dfrac{2i^{2}(\ln|e^{it}|+it)\sqrt{\cos(t)}e^{it/2}}{e^{it/2}(e^{-it}-e^{it})}\, dt= \int_{0}^{\frac{\pi}{2}}\dfrac{t\sqrt{\cos(t)}}{\sin(t)}\, dt. \end{gather*}
Now we proceed to \begin{gather*} \int_{C_3}f(z)\, dz = -\int_{0}^{1}\dfrac{2i\log(iy)\sqrt{\frac{1-y^2}{2}}}{\sqrt{(iy)}(1+y^2)}i\, dy =\\[2ex] \int_{0}^{1}\dfrac{2\left(\ln(y)+i\frac{\pi}{2}\right)\sqrt{\frac{1-y^2}{2}}}{\frac{1+i}{\sqrt{2}}\sqrt{y}(1+y^2)}\, dy = \int_{0}^{1}\dfrac{(1-i)\left(\ln(y)+i\frac{\pi}{2}\right)\sqrt{1-y^2}}{\sqrt{y}(1+y^2)}\, dy=\\[2ex] \int_{0}^{1}\dfrac{\left(\ln(y)+\frac{\pi}{2}\right)\sqrt{1-y^2}}{\sqrt{y}(1+y^2)}\, dy+i\cdot\int_{0}^{1}\dfrac{\left(\frac{\pi}{2}-\ln(y)\right)\sqrt{1-y^2}}{\sqrt{y}(1+y^2)}\, dy. \end{gather*} Now we extract the real part of (1). We get \begin{equation*} \int_{0}^{\frac{\pi}{2}}\dfrac{t\sqrt{\cos(t)}}{\sin(t)}\, dt+\int_{0}^{1}\dfrac{\left(\ln(y)+\frac{\pi}{2}\right)\sqrt{1-y^2}}{\sqrt{y}(1+y^2)}\, dy =0. \end{equation*} Thus \begin{equation*} S=-\dfrac{\pi}{2}\int_{0}^{1}\dfrac{\sqrt{1-y^2}}{\sqrt{y}(1+y^2)}\, dy =[y=s^2] =-\pi\int_{0}^{1}\dfrac{\sqrt{1-s^4}}{1+s^4}\, ds \end{equation*} However, \begin{equation*} \dfrac{\sqrt{1-s^4}}{1+s^4} =\dfrac{1}{2}\dfrac{\sqrt{\frac{1}{s^2}-s^2}}{\left(\frac{1}{s^2}-s^2\right)^2+4}\cdot 2\dfrac{1+s^4}{s^3} = \dfrac{1}{2}\dfrac{\sqrt{u}}{u^2+4}\cdot \dfrac{du}{ds}\cdot(-1). \end{equation*} where \begin{equation*} u=\dfrac{1}{s^2}-s^2. \end{equation*} Consequently, if we make the substitution $u=\dfrac{1}{s^2}-s^2$ then \begin{equation*} S=-\pi\int_{0}^{\infty}\dfrac{\sqrt{u}}{2(u^2+4)}\, du = [u=2\sqrt{v}]=-\dfrac{\pi}{4\sqrt{2}}\int_{0}^{\infty}\dfrac{v^{\frac{3}{4}-1}}{(1+v)^{\frac{3}{4}+\frac{1}{4}}}\, dv. \end{equation*} Here we recognize the beta $\mathrm{B}$ function. See
https://en.wikipedia.org/wiki/Beta_function
If we combine this with Euler's reflection formula we get \begin{gather*} S= -\dfrac{\pi}{4\sqrt{2}}B\left(\frac{3}{4},\frac{1}{4}\right)=-\dfrac{\pi}{4\sqrt{2}}\dfrac{\Gamma(\frac{3}{4})\Gamma(\frac{1}{4})}{\Gamma(1)} =\\[2ex] -\dfrac{\pi}{4\sqrt{2}}\Gamma\left(\frac{3}{4}\right)\Gamma\left(1-\frac{3}{4}\right)= -\dfrac{\pi}{4\sqrt{2}}\dfrac{\pi}{\sin(\frac{3\pi}{4})} = -\dfrac{\pi^2}{4}. \end{gather*}
Best Answer
By symmetry:
$$I=4\int _0^1\int _0^1\int _0^x\int _0^z\sqrt{(z-w)^2+(x-y)^2} \, dw \, dy \, dx dz$$
Substituting $w \to z w, \quad y \to x y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2(1-w)^2+x^2(1-y)^2}~ x z \, dw \, dy \, dx dz$$
Substituting $w \to 1-w, \quad y \to 1-y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2 w^2+x^2 y^2}~ x z \, dw \, dy \, dx dz$$
Substituting $x^2 = u, \quad z^2 = v$:
$$I=\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{v w^2+u y^2}~ \, dw \, dy \, du dv$$
Substituting $v w^2=p, \quad u y^2=q$:
$$I=\int _0^1\int _0^1\int _0^{y^2}\int _0^{w^2} \sqrt{p+q}~ \, dp dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{2}{3} \int _0^1\int _0^1\int _0^{y^2}\left((q+w^2)^{3/2}-q^{3/2} \right) dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{4}{15} \int _0^1\int _0^1\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
By symmetry:
$$I=\frac{8}{15} \int _0^1\int _0^w\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
Substituting: $y = w s$:
$$I=\frac{8}{15} \int _0^1\int _0^1 w^2 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds dw}{s^2 }$$
$$I=\frac{8}{45} \int _0^1 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds }{s^2 }$$
The last single integral can be evaluated by different methods, and the result agrees with the one from the OP.
One good method here is integration by parts. I leave this as an exercise.