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


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


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


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$.

  • 4
    $\begingroup$ 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. $\endgroup$ – MarcoB Feb 26 at 18:52
  • $\begingroup$ 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... $\endgroup$ – Matt Stein Mar 1 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.