Numerical integration: gravitational force

What I would like to accomplish is to numerically integrate the following equation:

$$\nabla^{2}\Phi(r)=\rho(r)$$,

to get $$\nabla\Phi(r)$$. $$r$$ is the radial distance from the origin. Suppressing the $$r$$-direction unit vector for a moment, we have

$$r^{2}\nabla\Phi(t)-r_{i}^{2}\nabla\Phi(r_{i})=\int_{r_{i}}^{r}\rho(r^{\prime})r^{\prime2}dr^{\prime}$$,

where $$r_i$$ is some arbitrary point to give a boundary condition. For concreteness, let us consider $$\rho$$ given as

$$\rho(r)=\frac{1}{r(1+r)^2}$$.

This can be integrated analytically, so I can check that I am losing any precision or not.

Because of the physical problem that I am interested in, I am using a logarithmically spaced grid in the $$r$$-direction. Also, I want to integrate from the outer radius to the inner radius by giving a boundary condition at the outermost grid point.

Using the trapezoidal rule, integration may be done as

ri = -3.;
rf = 1.;
NR = 1000;
δR = (rf - ri)/NR;
ρNFW[r_] := 1/((r)*(1 + r)^2);
ϕi[r_] := (-1 + 1/(1 + r) + Log[1 + r])/r^2;
Do[ρsol[j] = ρNFW[10^(ri + δR*j)], {j, 0, NR}];
Do[ϕgrid[j] = 10^(3*(ri + δR*j))*ρsol[j]*Log[10],{j,0, NR}];
Do[ϕintgrnd[j] = δR*1/2*(ϕgrid[j] + ϕgrid[j + 1]), {j, 0, NR - 1}];
ϕtest[NR] = 10^(2*(ri + δR*NR))*ϕi[10^(ri + δR*NR)];
Do[ϕtest[NR - j] = (ϕtest[NR - (j - 1)] - ϕintgrnd[NR - j]), {j, 1, NR}];
Do[ϕsol[j] = 1/10^(2*(ri + δR*j))*ϕtest[j], {j, 0, NR}];
dataϕ = Table[{ri + δR*j, ϕsol[j]}, {j, 0, NR}];
ϕint = Interpolation[dataϕ, InterpolationOrder -> 2];
LogPlot[{ϕi[10^a], ϕint[a]}, {a, ri, rf}, PlotRange -> All, AxesLabel -> {"r", "ϕ"}, LabelStyle -> {Black, Bold, Medium}]


We can see that this method starts to deviate from the analytic result around $$r=10^{-3}$$. This does not happen if I integrate from inside to outside. Origin of this problem is that when I integrate from outside to inside, a delicate cancelation of two large numbers is involved. This gets more severe as we go to a smaller radius.

For the problem I am dealing with, it is better to integrate from outside to inside. Would there be a good numerical method to integrate to a radius as small as $$r=10^{-5}$$? Due to the limited RAM, it is better for me to keep the number of grid points smaller than $$NR=1000$$.

• Would it not be possible to set up your problem within the NDSolve framework? You may have better numerical stability with the built in method. – MarcoB Feb 26 at 18:52
• Looks like you have an issue with parts prior to about the 45th element in your lists (which is obvious from your plot). I've worked it back to [Phi]test[j] which you can see is negative from 0 to about 45, then crosses over. Is this what you intended? I don't have time to back it out all the way, since this function relies on all the others... – Matt Stein Mar 1 at 16:04