# Mathematica integral wrong

I need to evaluate the following integral:

 Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z}, Ball[{0, 0, 0}, R]] ,
Element[{x1, y1, z1}, Ball[{0, 0, 0}, R1]],
Assumptions -> R > 0 && R1 > 0 ]


and Mathematica says the integral is 0, which is obviusly wrong. Am I doing something wrong?

• Wait, what does this mean? I want a double integral in $x,y,z$ and $x1, y1, z1$, with the former variables in a ball of radius $R$ and the latter in a ball of radius $R1$. In your suggestions, $R1$ doesn't show up Commented Jul 7 at 12:03
• I think his is looking for the gravitational energy between two balls, one placed inside the other—thus the double integral. @tommy1996q, I think there will be a “ball” of singularities in that intégrantd. But, I think your zero is probably not correct. (This comment referred to an earlier but deleted comment—i think the OP’s comment above also referred to the deleted comment) Commented Jul 7 at 12:04
• @CraigCarter nope, it's an entirely different problem, but it doesn't really matter. The point is that it's a very simple integration (conceptually at least) and I do not understand how Mathematica could give an answer so wrong. Trivially, I am integrating something which is always greater than 0. Commented Jul 7 at 12:07
• Greater than zero, yes, but infinite at points in the smaller ball—i think. The 0 seems like a bug in any case. Commented Jul 7 at 12:09
• @tommy1996q "The point is that it's a very simple integration..." So what is the correct answer when you know it is very simple? Commented Jul 7 at 16:15

Not an answer, but this won't fit in comment:

    Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z},
RegionDifference[Ball[{0, 0, 0}, bigR], Ball[{0, 0, 0}, littleR]]],
Element[{x1, y1, z1}, Ball[{0, 0, 0}, littleR]],
Assumptions -> 0 < littleR < bigR]  (*gives 0*)


and

Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z},
RegionDifference[Ball[{0, 0, 0}, bigR],
Ball[{0, 0, 0}, littleR - epsilon]]],
Element[{x1, y1, z1}, Ball[{0, 0, 0}, littleR]],
Assumptions -> 0 < littleR < bigR && 0 < epsilon < littleR] (*gives 0*)


and

NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z},
RegionDifference[Ball[{0, 0, 0}, 2], Ball[{0, 0, 0}, 1 - 0.1]]],
Element[{x1, y1, z1}, Ball[{0, 0, 0}, 1 - 0.1]]
] (*gives 61.2 but with warnings about accuracy*)


I think that 0 is an incorrect result and worth a bug report.

NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z}, Ball[{0, 0, 0}, 1]],
Element[{x1, y1, z1}, Ball[{0, 0, 0}, 2]]] (*gives 71.6 with warnings*)


but I don't trust that result at all.

• Thanks! Later I'll look more into it. Unluckily I'd really need some analytic solution. But yeah, 0 is most defnitily wrong. I'll report the bug Commented Jul 7 at 15:46

The result is totally wrong. This integral isn't too hard to do in spherical coordinates.

Start with the interaction between two spherical shells of radii $$r$$ and $$r_1$$. Let $$\{x,y,z\}=\{0,0,r\}$$ (we pick the $$z$$-axis such that it contains the point on the first spherical shell) and $$\{x_1,y_1,z_1\}=r_1\{\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\}$$. The integrand becomes

1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2] /.
{x -> 0, y -> 0, z -> r,
x1 -> r1 Sin[θ] Cos[φ], y1 -> r1 Sin[θ] Sin[φ], z1 -> r1 Cos[θ]} // FullSimplify
(*    1/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]]    *)


We have a Jacobian of $$4\pi r^2$$ for the point on the first sphere (the point on the first sphere can be chosen from anywhere on the ball's surface) and $$2\pi r_1^2\sin\theta$$ for the point on the second sphere. The integral is therefore, depending on the ordering of the radii,

Assuming[0 < r < r1,
Integrate[((4 π r^2) (2 π r1^2 Sin[θ]))/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]],
{θ, 0, π}]]
(*    16 π^2 r^2 r1    *)

Assuming[0 < r1 < r,
Integrate[((4 π r^2) (2 π r1^2 Sin[θ]))/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]],
{θ, 0, π}]]
(*    16 π^2 r r1^2    *)


Together, this result is the shell–shell interaction:

shell[r_, r1_] = 16 π^2 r r1 Min[r, r1]


You see that this interaction is not singular when $$r=r_1$$, which helps with the integration over the balls:

Assuming[0 < R < R1,
Integrate[shell[r, r1], {r, 0, R}, {r1, 0, R1}]] // FullSimplify
(*    -(8/15) π^2 R^3 (R^2 - 5 R1^2)    *)

Assuming[0 < R1 < R,
Integrate[shell[r, r1], {r, 0, R}, {r1, 0, R1}]] // FullSimplify
(*    -(8/15) π^2 R1^3 (-5 R^2 + R1^2)    *)


Together, the ball–ball interaction is therefore

ball[R_, R1_] = 8/15 π^2 Min[R, R1]^3 (5 (R^2 + R1^2) - 6 Min[R, R1]^2);


So, definitely not zero.

Let's try it out numerically:

With[{R = 0.3, R1 = 0.7},
{NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
Element[{x, y, z}, Ball[{0, 0, 0}, R]],
Element[{x1, y1, z1}, Ball[{0, 0, 0}, R1]]],
ball[R, R1]}]

(* ...lots of complaining about accuracy... *)

(*    {0.336256, 0.335409}    *)


seems to work fine.

## The result is mathematically correct, therefore there is no bug in Mathematica.

What you are trying to compute is this integral $$\mathcal{I}=\int_{B_R^3}\int_{B_{R_1}^3}\frac{1}{|\mathbf{x}-\mathbf{y}|}d\mathbf{x}^3d\mathbf{y}^3$$ We can use Legendre expansion that reads $$\frac{1}{|\mathbf{x}-\mathbf{y}|}=\sum\limits_{n=0}^{\infty}\frac{r_{<}^{n}}{r_{>}^{n+1}}P_n(\cos(\gamma))$$ where $$r_{<}$$ is the smaller of $$|\mathbf{x}|$$ and $$|\mathbf{y}|$$ and $$r_{>}$$ is the larger.

Let us choose a fixed spherical coordinate system for $$\mathbf{x}$$ such that $$\theta,\phi$$ are the relevant spherical angles (and $$r_x$$ is the norm). In comparison, we will choose a relative spherical coordinate system for $$\mathbf{y}$$ such that its $$z-$$axis will always point towards $$\mathbf{x}$$: the relevant angles will then be $$\gamma$$ (angle between these vectors) and $$\rho$$ (and $$r_y$$ is the norm). In this coordinate system, we then have $$\mathcal{I}=\sum\limits_{n=0}^{\infty}\int\limits_{0}^{R}r_x^2dr_x\int\limits_{0}^{R_1}r_y^2dr_y\int\limits_{0}^\pi\sin(\theta)d\theta\int\limits_{0}^{2\pi}d\phi\int\limits_{0}^\pi\sin(\gamma)d\gamma\int\limits_{0}^{2\pi}d\rho\frac{r_{<}^{n}}{r_{>}^{n+1}}P_n(\cos(\gamma))$$ which is basically $$\mathcal{I}=\sum\limits_{n=0}^{\infty}\left(\text{gamma independent terms}\right)\int\limits_{0}^\pi\sin(\gamma)P_n(\cos(\gamma))d\gamma$$

This integral is zero. You can also check it explicitly for various $$n$$, or directly using Mathematica as well:

FullSimplify[Integrate[Sin[\[Gamma]] LegendreP[n, Cos[\[Gamma]]], {\[Gamma], 0, Pi}, Assumptions -> Element[n, NonNegativeIntegers]],  Assumptions -> Element[n, NonNegativeIntegers]]


Therefore, the given integral $$\mathcal{I}$$ is zero.

• If you integrate something greater than zero, the integral is greater than zero. Commented Jul 12 at 20:12
• Not really. The properties of integrand carries over to the integration only if the integrand is locally integrable. In fact, $\frac{1}{|x|}$ is not even a proper function but actually a generalized function (also called distribution) in the integration domain, therefore it is not necessarily locally integrable (and indeed it is not). In other words, it is perfectly fine (and mathematically correct) that the result is zero. Commented Jul 12 at 23:30
• I think you're totally off the mark (-1). Commented Jul 13 at 10:54
• $1/|x|$ is locally integrable in $\mathbb{R}^3$, you are forgetting the jacobian term. Moreover, when $1/|x|$ is not locally integrable (as in the one dimensional case), it is not a distribution. Indeed a distribution is a continuous functional on smooth compactly supported function, and if you make $1/|x|$ act on a function which is equal to, say, 1 in a neighborhood of 0, this is not defined. Commented Jul 13 at 10:54
• I have not access to that book, but I am reasonably sure that there is one of the following caveats for which "$1/|x|$ is a distribution" (there may be other arguments, but these two are the only ones coming to my mind right now). First possibility: define the distribution on a subspace of $C^\infty _c (\mathbb{R})$, I think that functions that are 0 in 0 are enough. Otherwise, you can define the Cauchy principal value of $1/|x|$, which is something entirely different (check wiki). As for the local integrability in $\mathbb{R}^3$: you pass to spherical coordinates and the singularity vanishes Commented Jul 13 at 13:28