Método 1
Esta integral es realmente muy fácil de evaluar porque para cualquier $\vec{p} \in \mathbb{R}^3$,
la función de $\phi(\vec{x}) = \frac{1}{|\vec{x}-\vec{p}|}$ satisface la
De Laplace de la ecuación de distancia de $\vec{p}$.
$$\nabla^2 \phi(\vec{x}) = 0 \quad\text{ for } x \in \mathbb{R}^3\setminus \{\vec{p}\}$$
Este tipo de función se llama función Armónica y uno de sus más útiles de la propiedad es el llamado valor medio de la propiedad:
Si $\phi$ es una función armónica a través de una bola de $B(\vec{q},R)$ centrado
en $\vec{q}$ radio $R$, luego
$$\phi(\vec{p}) = \frac{1}{4\pi R^2} \int_{\parcial B(\vec{p},R)}
\phi(\vec{x}) dS(\vec{x}) = \frac{3}{4\pi R^3} \int_{B(\vec{p},R)}
\phi(\vec{x}) dV(\vec{x})$$
es decir, la media de $\phi(\vec{x})$ más de la superficie y el cuerpo de la bola de $B(\vec{q},R)$
son iguales al valor de $\phi(\vec{x})$ en el centro. Aplicar esto a la integral
en la mano, uno puede leer el valor de la integral como:
$$\frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$
Método 2
Si uno insiste en que para evaluar la integral utilizando polar de coordenadas, una manera de hacerlo es utilizar
el llamado multipolo de expansión.
Deje $\vec{p} = (a,b,c)$ $\theta$ ser el ángulo entre los dos
vectores $\vec{x}$$\vec{p}$.
Deje $r_{<} = \min( |\vec{x}|, |\vec{p}| )$$r_{>} = \max(|\vec{x}|,|\vec{p}|)$, luego
$$\frac{1}{|\vec{x} - \vec{p}|} = \sum_{n=0}^{\infty}\frac{r_{<}^n}{r_{>}^{n+1}} P_n(\cos\theta)$$
donde $P_n(\cos\theta)$ son los de Legendre de la
polinomios.
En nuestros casos, $\vec{p}$ está fuera de la esfera, el triple integral puede escribirse como
$$\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)$$
El uso de
$$\int_{0}^{\pi} P_n(\cos\theta)\sin\theta d\theta = \begin{cases}2,&n = 0\\0,&n > 0\end{cases}$$
se puede evaluar la integral triple como:
$$\frac{4\pi}{|\vec{p}|} \left( \int_0^R^2 dr \right)
= \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$
Método 3
Si uno realmente desea evaluar la integral en la más elemental forma, entonces uno
debe usar la simetría rotacional para simplificar el problema. Supondremos $(a, b, c) = (0,0,-p)$ y en términos de cilíndrico de coordenadas polares, el triple de la integral se convierte en:
$$\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}$$