Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top


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}]

Mathematica graphics

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][[1]]

Mathematica graphics

The answer involves ProductLog. I can also be done numerically as:

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

r /. r2[1/2]

(* ~ 2.5 *)


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}]

Mathematica graphics

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

Mathematica graphics

 Table[{e, r /. r2[e]}, {e, 1/10, 1 - 1/10, 1/50}]//Quiet // ListLinePlot

Mathematica graphics

share|improve this question
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
Thanks for the tip! – chris Nov 4 '13 at 21:22
up vote 7 down vote accepted

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[1],-Power[<<2>>] y]/y)]+ProductLog[C[1],-E^-y y]}. A likely reason for this is that the solution set depends on branch cuts of Mathematica functions. >>

C[1] \[Element] Integers && 1 + r != 0 && r != 0 && y != 0 && 
  0 == y + Log[-(ProductLog[C[1], -E^-y y]/y)] + 
   ProductLog[C[1], -E^-y y] && r == (-y - ProductLog[C[1], -E^-y y])/y

$0$-th branch returns 0

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

However, $(-1)$-th branch returns expected result

(-y - ProductLog[-1, -E^-y y])/y /. y -> 1/2 // N
share|improve this answer
Thanks for the explanation! – chris Nov 4 '13 at 21:22

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.