Evaluating a complicated triple integral I was working on a problem of trying to determine synchronization of neurons and wanted to compare an analytical result to some numerical results I achieved. The analytical result depends upon evaluating this triple integral:
$$\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi} \sqrt{3+2\cos(\theta_1-\theta_2)+2\cos(\theta_2-\theta_3)+2\cos(\theta_3-\theta_1)} d\theta_1 d\theta_2d\theta_3$$
The triple integral comes from trying to evaluate $E\left[\sqrt{|e^{i\theta_1}+e^{i\theta_2}+e^{i\theta_3}|^2} \right]$ where the $\theta$'s are i.i.d uniform random variables on the unit circle.
Wolfram alpha wasn't able to produce a value so I question if it's even possible to evaluate, whether exactly or a decimal approximation.
 A: It shouldn't be too hard to approximate this numerically.  I did a quick Monte Carlo integration in R, taking $n = 10^8$ samples:
set.seed(42)
n = 10^8
theta1 = runif(n, 0, 2*pi)
theta2 = runif(n, 0, 2*pi)
theta3 = runif(n, 0, 2*pi)
integrand = sqrt(3+2*cos(theta1-theta2)+2*cos(theta2-theta3)+2*cos(theta3-theta1))
mu = mean(integrand)
sigma = sd(integrand)

gives that the estimated mean of the integrand is $\mu = 1.574591$ and the standard deviation of the integrand is $\sigma = 0.7215917$. This gives a standard error for that estimate of $\sigma/\sqrt{n} = 0.0000722$.  
Thus we get an estimated value for the integrand of 
$$(8\pi^3) \left( \mu \pm {\sigma \over \sqrt{n}} \right)$$
since the region of integration has volume $8\pi^3$.    This is $390.578 \pm 0.018$.  
As an analytical note, your integrand is unchanged when $\theta_1, \theta_2, \theta_3$ are rotated by the same angle around the unit circle.  So you may as well set $\theta_3 = 0$.  Let $\theta_3 = 0$ and you get
$$ \int_0^{2\pi} \int_0^{2\pi} \int_0^{2\pi} \sqrt{3 + 2 \cos (\theta_1 - \theta_2) + 2 \cos \theta_2 + 2 \cos \theta_1} \: d\theta_3 d\theta_2 d\theta_1$$
and since $\theta_3$ isn't in the integrand and you're working over a rectangular region, this is
$$ 2\pi \int_0^{2\pi} \int_0^{2\pi} \sqrt{3 + 2 \cos (\theta_1 - \theta_2) + 2 \cos \theta_2 + 2 \cos \theta_1} \: d\theta_2 d\theta_1$$
This probably isn't easier to do analytically but I'd think a double integral would be easier than a triple integral numerically.
