Where is this difference between the solutions of DSolve and NDSolve coming from?

Question

I am trying to solve this DE:

DSolve[{r'[t] == r[t]^2/3 (1 - r[t]), r[0] == 2}, r[t], 
  t] // Simplify

For which I get the solution

InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][
 1/2 + I \[Pi] - t/3 - Log[2]]

After plotting the solution

Plot[InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][
  1/2 + I \[Pi] - t/3 - Log[2]], {t, 0, 20}, PlotRange -> All, 
 PlotPoints -> 100]

I get

However, the numerical solution using NDSolve results in:

How can I fix my DSolve result? I suspect the root finding process in the inverse function is a candidate? Because when I run

Table[{t, 
  InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][
   1/2 + I \[Pi] - t/3 - Log[2]]}, {t, 0, 20, 0.5}]

some of the elements are not evaluated:

{{0., InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][-0.193147 + 3.14159 I]}, {0.5, InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][-0.359814 + 3.14159 I]},...

I am running MMA 10.3 on Windows 7.

Answer 1

Indeed, the numerical inversion of the DSolve solution,

s = DSolveValue[{r'[t] == r[t]^2/3 (1 - r[t]), r[0] == 2}, r[t], t] // Simplify
(* InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][1/2 + I π - t/3 - Log[2]] *)

used to plot r as a function of t in the question is having difficulties. An alternative approach is to plot t as a function of r. First, extract the two arguments of the InverseFunction.

fun = Head[s][[1]][r]
(* 1/r + Log[1 - r] - Log[r] *)

arg = s[[1]]
(* 1/2 + I π - t/3 - Log[2] *)

equate the two, and solve for t.

a = (t /. Flatten@Solve[fun == arg, t]) // Expand
(* 3/2 + 3 I π - 3/r - 3 Log[2] - 3 Log[1 - r] + 3 Log[r] *)

Finally,

ParametricPlot[{a, r}, {r, 1, 2}, AspectRatio -> 1/GoldenRatio, 
    PlotRange -> {{0, 20}, All}, AxesOrigin -> {0, .95}]

which reproduces the numerical solution in the question, as desired.

Answer 2

It could also be view as a precision problem:

Plot[
 InverseFunction[Log[1 - #1] - Log[#1] + 1/#1 &][1/2 + I π - t/3 - Log[2]],
 {t, 0, 20},
 PlotRange -> All, PlotPoints -> 100, WorkingPrecision -> 50, Frame -> True]

Of course InverseFunction is much slower than @bbgodfrey's parametric approach.

Answer 3

You can get a direct anylytical solution for t[r] if you make use of the identiy

tsol = t /. First@DSolve[{1/t'[r] == r^2/3 (1 - r), t[2] == 0}, t, r]

(*   Function[{r}, (1/(2 r))
3 (-2 + r + 2 I \[Pi] r - 2 r Log[2] - 2 r Log[1 - r] + 2 r Log[r])]   *)

tsol[r] // FullSimplify[#, 1 < r < 2] &

(*   3/2 - 3/r - 3 Log[-1 + r] + Log[r^3/8]   *)

ParametricPlot[{tsol[r], r}, {r, 1, 2}, AspectRatio -> .8, 
  PlotRange -> {{0, 20}, {.95, 2.05}}, AxesOrigin -> {0, .95}
]