6
$\begingroup$

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), m1=1

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

enter image description here

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

enter image description here

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

$\endgroup$
1
  • $\begingroup$ 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) $\endgroup$
    – user49048
    Mar 16, 2022 at 19:21

2 Answers 2

6
+50
$\begingroup$

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

Mathematica graphics

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

$\endgroup$
5
$\begingroup$

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
*)
$\endgroup$
2
  • $\begingroup$ Thanks for your helpful comment! $\endgroup$
    – korni1990
    Mar 20, 2022 at 21:59
  • $\begingroup$ @korni1990 You're welcome! It's good to have numerically sound code. $\endgroup$
    – Michael E2
    Mar 21, 2022 at 0:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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