4
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Nov 4 '13 at 21:07
7
$\begingroup$

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
$\endgroup$

Your Answer

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

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