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

Question

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 missing motivation, 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
Okay, I retracted my close vote. Not sure if that was the right thing to do because this also borders on being a duplicate, in addition to being not particularly about Mathematica. But there is an answer by @Feyre that gives a good explanation of the behavior, so maybe best to keep it open. — Daniel Lichtblau, Aug 15 '16 at 18:21
@DanielLichtblau The fact that ans4 (from Simplify) is different from ans5 (from FullSimplify) makes this question very much about particulars of Mathematica. — Jerry Guern, Aug 15 '16 at 22:56
@DanielLichtblau To put it another way, someone who understands precision/roundoff issues perfectly would still need to find out somehow what MMA functions will cause/prevent those issues. — Jerry Guern, Aug 15 '16 at 23:00
What would be the best way to accomplish the simplification efficiently without causing this issue? — Jerry Guern, Aug 16 '16 at 02:12
I have to say that I( disagree with most of the recent comments. If you look at the results from Simplify et al, you see a structural breakdown into computations that involve basic arithmetic and root extraction. From there is is not so difficult to apply N[] to pieces and find out that some involve huge cancellation error. Also at that point, once the structure of the expression is made, it really becomes Mathematica-independent at least so far as double precision floating point arithmetic is concerned. — Daniel Lichtblau, Aug 16 '16 at 14:23
As for getting a form of result that is not prone to such problems, the thing that most comes to mind is RootReduce. It does not always give a "simplification", and it might be slow (but typically so is FullSimplify). Whatever else though, a Root object tends to be good for stability of numeric evaluation. — Daniel Lichtblau, Aug 16 '16 at 16:31
@DanielLichtblau A succinct statement of the issue would be: Simplify[] puts the expression into a form that causes N[] to return a wrong answer. Even after I saw the structure of what Simplify[] returns, it didn't occur to me that N[] would evaluate it naively. THAT is most definitely a Mathematica quirk that Mathematica users need to beware of and should be discussed on a Mathematica site. — Jerry Guern, Aug 16 '16 at 17:31
(1) That is absolute rot. SImplify puts the expression into a perfectly viable symbolic form. It is evaluated numerically pretty much the same way it would be done in any language. The fact that it happens to be susceptible to cancellation error is, once again, an issue of numerical math. A very well understood one, I might add. — Daniel Lichtblau, Aug 16 '16 at 18:20
(2) It's really not such a good idea to make the point that I never should have retracted the vote to close. — Daniel Lichtblau, Aug 16 '16 at 18:21

Feyre · Accepted Answer · 2016-08-15T13:04:38.650

6

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.

edited Aug 15 '16 at 13:04

answered Aug 15 '16 at 12:58

Feyre

8,597
2
27
46

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 02: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 https://mathematica.stackexchange.com/questions/123871/fullsimplify-providing-an-answer-not-valid-for-certain-parameters/123876#123876). 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 07:19

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

1 Answers1