# Solutions of NDSolve and DSolve suddenly deviating

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][];
c2 = Solve[v == 1, C];
v2[xp_] = v[xp] /. c2[];
c1 = Solve[v2' == 0, C][];
vm[xp_, m_] = v2[xp] /. c1 // FullSimplify;

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

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$$?

• I deleted my previous reply as I was able to verify your claim. I am truly sorry. I don't know what was messed up in the code (very red face)
– user49048
Mar 16, 2022 at 19:21

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. :)

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