I am trying to use Mathematica to check an integral from a book, and cannot get it done. The book talks about a spherically symmetric mass distribution in $R^3$ of the form

$\rho(r)=r^{-\gamma}, \gamma<3$

and it goes on by stating that the Newtonian potential to that distribution is given by

$\phi(r) = 4π\ln r + \mathrm{const.}, \gamma=2$


$\phi(r) = \frac{4π}{(2-\gamma)(3-\gamma)}r^{2-\gamma}+\mathrm{const.}, \gamma\neq2$ .

(There are a couple more constants in the book which I set to 1 here for simplicity.)

I wanted to verify this by directly evaluating the integral for the potential

$\phi(\mathbf x) = -\int_{R^3}\frac{\rho(\mathbf x')}{|\mathbf x-\mathbf x'|}\,d^3x' =-4π\int_0^\infty\frac{\rho(r')}{|r-r'|}r'^2\,dr'$ ,

and I am failing.

Here is my code to check the $\gamma=2$ case:

Assuming[{r > 0, g == 2}, -4 π Integrate[
   s^(-g) / Abs[r - s] s^2, {s, 0, ∞}]]

Mathematica gives me the message: "Integrate: Integral of 1/Abs[r-s] does not converge on {0,[Infinity]}"

My code for the other case is

Assuming[{r > 0, g != 2, g < 3}, -4 π Integrate[
   s^(-g) / Abs[r - s] s^2, {s, 0, ∞}]]

which also yields "Integrate: Integral of s^(2-g)/Abs[r-s] does not converge on {0, ∞}."

What am I missing?

  • 1
    $\begingroup$ One does problems of this sort with the help of the Gauss theorem. The way you try it, is much too difficult, and directly it will not work. $\endgroup$ – Alexei Boulbitch Mar 3 at 16:15
  • $\begingroup$ Hi @AlexeiBoulbitch , thanks for the answer! Sure, there might be simpler ways than directly evaluating the integral. I wanted to try this, since later I want to generalise this to more complex distributions and integrate numerically. But I found a mistake now: I cannot simply replace $|\mathbf x-\mathbf x'|$ in the 3D integrand by $|r-r'|$ in the 1D integrand! I will correct this and try again tomorrow and then report back! $\endgroup$ – Britzel Mar 3 at 22:42
  • $\begingroup$ FullSimplify[Integrate[s^(2-g)/Sqrt[(r-s)^2],s]] does this help? You can now try Limit and Series with Assuming. $\endgroup$ – Vitaliy Kaurov Mar 4 at 5:44
  • $\begingroup$ @Britzel If your (more complex) distributions are spherically- or cylindrically symmetric, the Gauss theorem approach is the way to implement. Do you know it? If you have an asymmetric problem you will have a bad time. I am just struggling against a comparable problem, where I need to integrate a Green function over a distribution. [Here][1] you can see the level of mathematical difficulties arising there. Have fun! [1]: mathematica.stackexchange.com/questions/215614/… $\endgroup$ – Alexei Boulbitch Mar 4 at 8:58
  • $\begingroup$ @AlexeiBoulbitch For the start they will be spherically symmetric, yes. Perhaps it is a good idea to go with a Gauß theorem approach for now, as you suggest. Eventually though they will become axisymmetric. But not for now. $\endgroup$ – Britzel Mar 4 at 13:16

The problem is with the expansion: to integrate $\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}$ over a volume with a spherically invariant factor (the charge density), the standard (and the most sensible trick) is to apply spherical harmonic decomposition:

$$\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)Y_{lm}^*(\theta',\phi')$$

With this expansion, the integral will be much simpler; in particular, the angular integral will immediately fix radial dependence. In fact, for spherically symmetric overall factor, the angular integration simply reads

$$\begin{align} \int d\Omega'\rho(r')\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=&\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)\rho(r')\int d\Omega'Y_{lm}^*(\theta',\phi')\\ =&\sum\limits_{l=0}^\infty\frac{r_<^l}{r_>^{l+1}}\left(\frac{4\pi}{2l+1}\right)\sum\limits_{m=-l}^{l}Y_{lm}(\theta,\phi)\rho(r')\sqrt{4\pi}\delta_l^0\delta_m^0\\ =&\frac{4\pi}{r_>}\rho(r') \end{align}$$

The radial integral is now straightforward:

$$\begin{align} \int d^3x'\frac{\rho(r')}{\lvert \mathbf{x'} - \mathbf{x}\rvert}=&\int r'^2 dr'\int d\Omega'\rho(r')\frac{1}{\lvert \mathbf{x'} - \mathbf{x}\rvert}\\ =&4\pi\int \rho(r')r' dr'\\ =&4\pi\int r'^{1-\gamma} dr' \end{align}$$

which yields your expected result.

| improve this answer | |
  • 1
    $\begingroup$ Fantastic, thanks! I will have to look back into spherical harmonics in order to understand the angular integration you performed. Will do so now! $\endgroup$ – Britzel Mar 4 at 13:18
  • 1
    $\begingroup$ Soner, can you show how you did this with the Wolfram Language? $\endgroup$ – CA Trevillian Mar 4 at 18:47
  • 1
    $\begingroup$ @CATrevillian I am unaware of how to do intermediate steps I wrote above with the Wolfram Language. $\endgroup$ – Soner Mar 5 at 4:47

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.