3
$\begingroup$

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.

|improve this question
$\endgroup$
  • 1
    $\begingroup$ Why are your limits in your code not the same as in the problem? $\endgroup$ – David G. Stork Nov 5 '15 at 0:15
  • 1
    $\begingroup$ Not that this should affect the sign, but which set of limits are you trying to impose? $\endgroup$ – IPoiler Nov 5 '15 at 0:15
  • $\begingroup$ made corrections.thanks. $\endgroup$ – BigM Nov 5 '15 at 0:16
  • $\begingroup$ Isn't the issue that this integral diverges very badly around the xyzstu origin? $\endgroup$ – evanb Nov 5 '15 at 0:26
  • 2
    $\begingroup$ 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. $\endgroup$ – IPoiler Nov 5 '15 at 0:35
8
$\begingroup$

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:

enter image description here

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

|improve this answer
$\endgroup$
7
$\begingroup$

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} *)
|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.