Triple Integral $ \iiint\limits_S\frac1{\sqrt{\left(x-a\right)^{2}+\left(y-b\right)^{2}+\left(z-c\right)^{2}}}\;dx\;dy\;dz $

I used spherical coordinates and couldn't calculate this triple integral. $$ \iiint_S\frac1{\sqrt{\left(x-a\right)^{2}+\left(y-b\right)^{2}+\left(z-c\right)^{2}}}\;dx\;dy\;dz $$

where $S$ is a solid sphere of radius $R$ and center at the origin, and $\left ( a,b,c \right )$ is a fixed point outside this sphere.

Method 1

This integral is actually very easy to evaluate because for any $\vec{p} \in \mathbb{R}^3$, the function $\phi(\vec{x}) = \frac{1}{|\vec{x}-\vec{p}|}$ satisfies the Laplace's equation away from $\vec{p}$.

$$\nabla^2 \phi(\vec{x}) = 0 \quad\text{ for } x \in \mathbb{R}^3\setminus \{\vec{p}\}$$

This sort of function is called Harmonic function and one of its most useful property is the so called mean value property:

If $\phi$ is a function harmonic over a ball $B(\vec{q},R)$ centered at $\vec{q}$ with radius $R$, then

$$\phi(\vec{q}) = \frac{1}{4\pi R^2} \int_{\partial B(\vec{q},R)} \phi(\vec{x}) dS(\vec{x}) = \frac{3}{4\pi R^3} \int_{B(\vec{q},R)} \phi(\vec{x}) dV(\vec{x})$$

i.e. the average of $\phi(\vec{x})$ over the surface and the body of the ball $B(\vec{q},R)$ are both equal to the value of $\phi(\vec{x})$ at the center. Apply this to the integral at hand, one can read off the value of the integral as:

$$\frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$

Method 2

If one insists to evaluate the integral using polar-coordinates, one way to do it is to use the so called multipole expansion.

Let $\vec{p} = (a,b,c)$ and $\theta$ be the angle between the two vectors $\vec{x}$ and $\vec{p}$.
Let $r_{<} = \min( |\vec{x}|, |\vec{p}| )$ and $r_{>} = \max(|\vec{x}|,|\vec{p}|)$, then $$\frac{1}{|\vec{x} - \vec{p}|} = \sum_{n=0}^{\infty}\frac{r_{<}^n}{r_{>}^{n+1}} P_n(\cos\theta)$$ where $P_n(\cos\theta)$ are the Legendre's polynomials.

In our cases, $\vec{p}$ is outside the sphere, the triple integral can be rewritten as

$$\sum_{n=0}^{\infty} \frac{2\pi}{|\vec{p}|^{n+1}} \left( \int_0^R r^{n+2} dr \right) \left( \int_{0}^{\pi} P_n(\cos\theta) \sin\theta d\theta \right)$$

Using $$\int_{0}^{\pi} P_n(\cos\theta)\sin\theta d\theta = \begin{cases}2,&n = 0\\0,&n > 0\end{cases}$$

one can evaluate the triple integral as:

$$\frac{4\pi}{|\vec{p}|} \left( \int_0^R r^2 dr \right) = \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$

Method 3

If one really want to evaluate the integral in the most elementary manner, then one should use rotational symmetry to simplify the problem. We will assume $(a, b, c) = (0,0,-p)$ and in terms of cylindrical polar coordinate, the triple integral becomes:

$$\begin{align} & 2 \pi \int_{-R}^{R} \left( \int_0^{\sqrt{R^2-z^2}} \frac{rdr}{\sqrt{r^2 + (z+p)^2}} \right) dz\\ = & 2 \pi \int_{-R}^{R} \left[\sqrt{r^2+(z+p)^2}\right]_0^{\sqrt{R^2-z^2}} dz\\ = & 2 \pi \int_{-R}^{R} \left(\sqrt{R^2+2zp+p^2} - (z+p)\right) dz\\ = & 2 \pi \left[ \frac{1}{3p} \sqrt{R^2+2zp+p^2}^3 - (\frac{z^2}{2} + p z)\right]_{-R}^R\\ = & 2 \pi \left( \frac{1}{3p}\left( (p+R)^3 - (p-R)^3 \right) - 2pR \right)\\ = & \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}} \end{align}$$