# Error in NIntegrate command

I am trying to evaluate (numerically) the integral $$\int_0^1\int_0^1\int_0^1\int_0^1\int_0^1\int_0^1\frac{dx\,dy\,dz\,ds\,dt\,du}{(x - s)^2 + (y - t)^2 + (z - u)^2}$$ with Mathematica using the following

NIntegrate[1/((x - s)^2 + (y - t)^2 + (z - u)^2),
{x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {s, 0, 1}, {t, 0, 1}, {u, 0, 1}]


but surprisingly, it keeps giving me a negative answer! It would be highly appreciated, if someone could give an alternative method (command) to evaluate this integral.

• Why are your limits in your code not the same as in the problem? Commented Nov 5, 2015 at 0:15
• Not that this should affect the sign, but which set of limits are you trying to impose? Commented Nov 5, 2015 at 0:15
– BigM
Commented Nov 5, 2015 at 0:16
• Isn't the issue that this integral diverges very badly around the xyzstu origin? Commented Nov 5, 2015 at 0:26
• I'm kind of concerned by the fact that Integrate[1/((x - s)^2 + (y - t)^2 + (z - u)^2), x] evaluates quickly but Integrate[1/((x - s)^2 + (y - t)^2 + (z - u)^2), {x,0,1}] can be done by hand before Mathematica finishes. Commented Nov 5, 2015 at 0:35

I'm not confident of this answer, but it seems worth presenting:

NIntegrate[2^3/((x - s)^2 + (y - t)^2 + (z - u)^2), {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {s, 0, x}, {t, 0, y}, {u, 0, z}, AccuracyGoal -> 16] // Quiet

(* 5.63378 *)

The basic insight is that each "paired" integrations (e.g., $x$ and $s$; $y$ and $t$; $z$ and $u$) is over the unit square. Each such pair can be replaced by a symmetric integration over half that square (e.g., $0 < x < 1$ with $0 < s < x$). Thus there will be an undercounting of $2^3 = 8$.

One can extend this core idea to integrating each "pair" of variables over a smaller portion of their respective unit square and multiply the overall integral accordingly.

Anyway... perhaps this will point others in the proper direction.

I now believe my answer is correct.

Here is why: Make a change of variables: $a = x - s$; $b = y - t$ and $c = z - u$. Then the integral over the unit square defined by $x$ and $s$ can be recast as a one-dimensional integral over the new variable $0 \le a \le 1$. However, there is a scaling factor that is needed which you can see if you draw the unit square and consider symmetry. There is hence a factor of $2 (1 - a)$ and likewise for the $b$ and $c$ variables.

NIntegrate[(8 (1 - a) (1 - b) (1 - c))/(a^2 + b^2 + c^2), {a, 0, 1}, {b, 0, 1}, {c, 0, 1}, AccuracyGoal -> 16]

(* 5.63371 *)

If someone has the time, he or she could go one step further and replace the integration over the cube defined by $a$, $b$ and $c$ to a one-dimensional integration over a variable $q$ that goes on the diagonal of that cube, from $(a,b,c) = (0,0,0)$ to $(1,1,1)$. Be careful with scaling and normalization. If we're lucky, we might get a closed-form answer (presumably involving inverse trigonometric functions). Here's a figure illustrating the new variable $q$ in red:

Wouldn't that be cool? A six-dimensional numerical integration replaced by a one-dimensional symbolic integration?!

I'm 95% confident that David G. Stork's answer is in the ballpark:

n = 100000000;
x = RandomReal[1, n]; y = RandomReal[1, n]; z = RandomReal[1, n];
s = RandomReal[1, n]; t = RandomReal[1, n]; u = RandomReal[1, n];
r = 1/((x - s)^2 + (y - t)^2 + (z - u)^2);
xbar = Mean[r];
xsd = StandardDeviation[r];
lowerCL = xbar - 1.96 xsd/n^0.5;
upperCL = xbar + 1.96 xsd/n^0.5;

{lowerCL, upperCL}

(* {5.60976, 5.66793} *)