# Numerical errors/inaccuracies in ProductLog

Context

In cosmology, a fairly accurate model to describe the gravitational potential, $\psi(r)$ of dark matter halos is given by $\psi( r)=\log(1+r)/r$.

Plot[Log[1 + r]/r, {r, 0.01, 4}] In this context it is of interest to find the radius at which the potential is equal to some energy.

Now Mathematica seems happy with solving this implicit equation:

r1 = Solve[y == Log[1 + r]/r, r][] The answer involves ProductLog. I can also be done numerically as:

r2[e_] := NSolve[e == Log[1 + r]/r, r][]

r /. r2[1/2]

(* ~ 2.5 *)


Question

But, then why does

r /. r1 /. y -> 1/2


return $0$? Why does the plot below returns complete non-sense?

Plot[r /. r1 // Evaluate, {y, 0, 1}] Finally, why do these two plots succeed and fail respectively?

Table[{e, r /. r2[e]}, {e, 1/10, 1 - 1/10, 1/30}] //
Quiet // ListLinePlot Table[{e, r /. r2[e]}, {e, 1/10, 1 - 1/10, 1/50}]//Quiet // ListLinePlot • For the first problematic plot, the issue seems to be WorkingPrecision. Set it to $MachinePrecision to resolve the problem. The other cases seem to be more subtle, perhaps to do with an inconsistent choice of branches in different evaluations. – Oleksandr R. Nov 4 '13 at 21:07 ## 1 Answer Result depends on branch cuts of Mathematica functions r1 = Solve[y == Log[1 + r]/r, r]  Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >> {{r -> (-y - ProductLog[-E^-y y])/y}}  Reduce[y == Log[1 + r]/r, r]  Reduce::useq: The answer found by Reduce contains unsolved equation(s) {0==y+Log[-(ProductLog[C,-Power[<<2>>] y]/y)]+ProductLog[C,-E^-y y]}. A likely reason for this is that the solution set depends on branch cuts of Mathematica functions. >> C \[Element] Integers && 1 + r != 0 && r != 0 && y != 0 && 0 == y + Log[-(ProductLog[C, -E^-y y]/y)] + ProductLog[C, -E^-y y] && r == (-y - ProductLog[C, -E^-y y])/y $0$-th branch returns 0 (-y - ProductLog[- E^-y y])/y /. y -> 1/2  0  However,$(-1)\$-th branch returns expected result

(-y - ProductLog[-1, -E^-y y])/y /. y -> 1/2 // N

2.51286