Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

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 *)

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

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

1 Answer 1

up vote 6 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
0

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

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

Your Answer

 
discard

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.