Why does applying Simplify after Expand give an unexpected numerical result? [duplicate]

The last five lines of this code should give the same numerical result, but they don't.

z[n_, c_] := If[n > 1, z[n - 1, c]^2 + c, c];
expr = (z[6, c] - z[5, c])*c^-6;
b = Solve[expr == 0, c];
dz5 = D[z[5, c], c] /. b[[2]];
ans1 = dz5 // Expand // N
ans2 = dz5 // Simplify // N
ans3 = dz5 // N
ans4 = dz5 // Expand // Simplify // N
ans5 = dz5 // Expand // FullSimplify // N

(* -7.73335  -7.73335 -7.73335  0.  -7.73335*)


It's only the specific combination of Expand and Simplify that give a strage result. Is this a bug? Or is a numerics quirk I don't understand?

• Why not do RootReduce[dz5] first? In any event, ponder on the result of N[Apply[List, Numerator[dz5 // Expand // Simplify]/256]]. – J. M.'s discontentment Aug 15 '16 at 12:47
• I'm voting to close this question as off-topic because it is about (numerical) math and not Mathematica. – Daniel Lichtblau Aug 15 '16 at 13:58
• @DanielLichtblau In my opinion this seems like a good Mathematica question. I would vote to not close. – Jack LaVigne Aug 15 '16 at 16:46
• @JackLavigne Maybe I am missing something, but this to me appears to be a standard issue of numeric floating point cancellation. In which case it is in no way specific to Mathematica numerics. It has shown up here and here for example. But I am willing to be convinced this should stay open, if you care to elaborate. – Daniel Lichtblau Aug 15 '16 at 16:54
• @JackLavigne You mean the answer by Feyre? – Daniel Lichtblau Aug 15 '16 at 17:08

This is a matter of Precision. Your initial dz5 has infinite precision, the Expand[] and Simplify[] in that order only (which is why other orders of Expand[], Simplify[], FullSimplify[] work fine) generates three expressions of a - b Sqrt[c] in the numerator as part of the expression, which at the precision given by bare N[] are 0.Try calling with higher precision instead:

N[dz5 // Expand // Simplify, 50]


-7.733352194016861443064354917432174265254146586

Note that

N[dz5 // Expand // Simplify, 1]


even works (yields 8), but only because this forces N[] to work at higher precision in order to yield the single digit precision in the answer.

• Thanks for this answer. I appreciate the effort you put into editing and clarifying it. What do you think would be the best method for me to use, to accomplish what I was trying to accomplish, without running into this issue? – Jerry Guern Aug 16 '16 at 2:12
• @JerryGuern Because of cancellation and rounding errors, if you just want a numerical result, I think it's usually best not to use any functions like Expand or Simplify (partly because these can give different removable singularities too mathematica.stackexchange.com/questions/123871/…). If you want to have both a neat analytical expression and a numeric result, it's safest to take the numeric result from the original expression. I'd always use N[expr,n] rather than expr//N to be on the safe side. – Feyre Aug 16 '16 at 7:19