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