Erratic Systematic Wormhole Variance

It is known that the average distance,

$\overline{x} = \int_0^1 \int_0^1 \int_0^1 \int_0^1 \sqrt{(x-y)^2+(z-w)^2} \text{ d}x\text{ d}y\text{ d}z\text{ d}w=\frac{1}{15} (2+\sqrt{2}+5\ln(1+\sqrt{2}))$

The variance, is

$\begin{aligned} \sigma ^2 &= \int_0^1 \int_0^1 \int_0^1 \int_0^1 (\sqrt{(x-y)^2+(z-w)^2} - \overline{x} )^2 \text{ d}x\text{ d}y\text{ d}z\text{ d}w \\ &= \int_0^1 \int_0^1 \int_0^1 \int_0 ^1 (x-y)^2+(z-w)^2 -2\overline{x} \sqrt{(x-y)^2+(z-w)^2} + \overline{x} ^2 \text{ d}x\text{ d}y\text{ d}z\text{ d}w \\ &= \int_0^1 \int_0^1 \int_0^1 \int_0 ^1 (x-y)^2+(z-w)^2 \text{ d}x\text{ d}y\text{ d}z\text{ d}w -2\overline{x} \int_0^1 \int_0^1 \int_0^1 \int_0 ^1 \sqrt{(x-y)^2+(z-w)^2} \text{ d}x\text{ d}y\text{ d}z\text{ d}w + \overline{x} ^2 \\ &= \frac{1}{3} -2\overline{x} \int_0^1 \int_0^1 \int_0^1 \int_0^1 \sqrt{(x-y)^2+(z-w)^2} \text{ d}x\text{ d}y\text{ d}z\text{ d}w + \overline{x} ^2 \\ &= \frac{1}{3} -2\overline{x} ^2+\overline{x} ^2\\ &= \frac{1}{3} - \overline{x} ^2 \\ &= \frac{1}{3} - \frac{1}{225} (2+\sqrt{2}+5\ln(1+\sqrt{2}))^2 \\ &= \frac{75-(2+\sqrt{2}+5\ln(1+\sqrt{2}))^2}{225} \\ & \approx 0.061469707 \\ \end{aligned}$

Say the particle's initial position is (a,b) and the particle's final position is (x,y).

The displacement of the particle would be

$\sqrt{\left(x-a\right)^2+\left(y-b\right)^2}$

Since we want to find the variance, we need to first find the average over all possible a, b, x, y:

$\overline { x } =\frac { \sum\text{All cases} }{ \text{No. of cases} } = \frac{1}{1}\int _0^1\frac{1}{1}\int _0^1\frac{1}{1}\int _0^1\frac{1}{1}\int _0^1\sqrt{\left(x-a\right)^2+\left(y-b\right)^2}\, dy \, dx \, da \, db = 0.521405433309$

Note that numerical methods were used to calculate the integrals.

Now to find the variance:

$s^2=\frac{\sum(\text{All cases}-\overline{x})^2}{\text{No. of cases}}=\int _0^1\int _0^1\int _0^1\int _0^1\left(\sqrt{\left(x-a\right)^2+\left(y-b\right)^2}-\overline{x}\right)^2 \, dy \, dx \, da \, db=0.061469707449$

Interesting fact: If the question asks for the root mean square, the answer would be exactly $1/\sqrt{3}$ . Also this problem actually has a closed form

Monte Carlo code

import math
import random

u = 0.521 # from previous problem

for z in range(0,10):  # run experiment 10 times

    N = 10**7   # 10 million points - Monte Carlo

    Dsqav = 0.0

    for j in range(0,N):

        x1 = random.uniform(0.0,1.0)
        y1 = random.uniform(0.0,1.0)

        x2 = random.uniform(0.0,1.0)
        y2 = random.uniform(0.0,1.0)

        Dsq = (x1-x2)**2.0 + (y1-y2)**2.0

        Dsqav = Dsqav + Dsq

    EDsq = Dsqav / N

    print EDsq,(EDsq - u**2.0)  # Var = EDsq - u**2.0