Solutions of NDSolve and DSolve suddenly deviating

Question

I am trying to get solutions to the differential equation down in the code:

eom[m_] = f[xp]^2 + m f[xp]^(3/2) f'[xp] - f'[xp]^2 + f[xp] f''[xp];
m1 = 10;
v[xp_] = f[xp] /. DSolve[eom[m] == 0, f[xp], xp][[1]];
c2 = Solve[v[0] == 1, C[2]];
v2[xp_] = v[xp] /. c2[[2]];
c1 = Solve[v2'[0] == 0, C[1]][[1]];
vm[xp_, m_] = v2[xp] /. c1 // FullSimplify;

vnum[xp_] = 
  f[xp] /. NDSolve[{eom[m1] == 0, f[0] == 1, f'[0] == 0}, 
     f[xp], {xp, 0, 2}][[1]];

Plot[{vm[xp, m1], vnum[xp]}, {xp, 0, 2}]

While for $m1=1$ both solutions agree (I dragged one down a bit),

For $m1=10$, they deviate a lot.

I wonder what goes wrong here? For context: I tried implementing the initial data into DSolve, which did not work. Hence I had to do it manually.

P.S. DSolve above gives two distinct solutions for C2. I checked with

Plot[eom[m1] /. f -> Function[xp, vm[xp, m1]], {xp, 0, 2}]

that only the one chosen above solves the ODE. With the above solution and $m1=10$, the plot is

So obviously the solution is not correct anymore. But yeah I am confused. Shouldn't the analytic result be correct for all $m$?

Answer 1

The solution is still correct, it's just a matter of precision:

Plot[{vm[xp, m1], vnum[xp]}, {xp, 0, 2}, WorkingPrecision -> 16, 
 PlotStyle -> {Automatic, {Red, Dashed}}]

WorkingPrecision is an option that everybody doing numeric calculation in Mathematica should know. :)

Answer 2

Had the output of v2[xp] /. c1 // FullSimplify been included, some might have been able to suggest a solution by inspection. It contains Erf[m] - Erf[m - xp/2], a somewhat obvious problem to those who know Erf. Try this:

Erf[m] - Erf[m - xp/2] /. m -> m1 /. xp -> {0., 1., 2.}
(* {0., 0., 0. } *)

For large arguments, which is for m only around 6 or greater (a bit over $8\sigma$ in terms of the normal distribution), Erf[m] is so close to 1 that it rounds to 1. exactly at machine-precision. This is why good math libraries contain Erfc[z] == 1 - Erf[z]:

Plot[{vm[xp, m1] /. Erf -> (1 - Erfc[#] &) // Evaluate, 
  vnum[xp]}, {xp, 0, 2}, PlotStyle -> {Automatic, Dashed}]

You should probably add /. Erf -> (1 - Erfc[#] &) to the definition of vm, but beware: FullSimplify simplifies expression complexity, not numerical instability. Don't apply FullSimplify after the substitution, for it will undo it. For example, do something like this:

vm[xp_, m_] = FullSimplify[v2[xp] /. c1] /. Erf -> (1 - Erfc[#] &)
(*
  E^(2 m xp - xp^2/2) /
   (1 + E^m^2 m Sqrt[π] (-Erfc[m] + Erfc[m - xp/2]))^2
*)