Re: FindRoot

*To*: mathgroup at smc.vnet.net*Subject*: [mg111171] Re: FindRoot*From*: Themis Matsoukas <tmatsoukas at me.com>*Date*: Thu, 22 Jul 2010 05:45:15 -0400 (EDT)

Bill, Your answer pointed me to the explanation for the wrong behavior, which I think is a bug in Mathematica. The clue is that F[X/.X->XX] and F[X]/.X->XX give different answers and the reason is that Mathematica's internal representation of F[X] as an analytic function of X is wrong. In the definition of F[X], the roots of the cubic are first sorted, then the function returns the ratio of root #3 to root #2 (in my original example I meant to use root #1 instead of #2 but this doesn't matter for the purposes of debugging). Apparently, the analytic form of F[X] (i.e. what you get when you execute X =.; F[X]]) does not sort the roots and instead takes the ratio of the wrong pair. At least, that's what the demo below shows. With F[X /. X -> 0.3], F is evaluated according to its definition at a specified numerical value of X and returns the correct answer (a number). With F[X]/.X->.3 a different thing happens: first, the definition of F is executed for a symbolic X and returns a symbolic expression for F in terms of X, which however is wrong; then the analytic form is evaluated at the specified value of X. The bottom line is this: F[X] is correct as long as X is a number; if it is an unevaluated symbol the answer is wrong. For this reason Plot, Table etc work but FindRoot, D, Integrate, etc do not. I would call this a bug in Mathematica since it produces an analytic expression for F that does not respect its definition and leads to apparent conflict between what should be equivalent representations, F[/.X->XX] and F[X]/.X->XX. ZEQ=-0.049973 p^2+(0.44137p-0.00873 p^2) z-z^2+z^3; {p1,p2}={0.133652,0.690320}; F[X_]:=( z123=Sort[z/.Solve[(ZEQ/.p->X)==0,z]]; z123[[3]]/z123[[2]]); XX=p1; dX=(p2-p1)/3.; imax=3; results={}; Do[ roots=Sort[z/.Solve[(ZEQ/.p->XX)==0,z]]; results=Append[results, Join[{XX}, roots, {F[XX], Re[F[X]/.X->XX]}]]; XX=XX+dX; ,{i,1,imax+1} ]; TableForm[results, TableHeadings->{None,{"XX", "z1", "z2", "z3","F[X]", "Re[F[X]/.X->XX]"}}] Themis