MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Kolmogorov-Smirnov 2-sample test
  • Next by Date: Cases vs. StringCases vs. Select and StringMatchQ vs. StringFreeQ
  • Previous by thread: Re: FindRoot
  • Next by thread: Scoping constructs Block, Module, ModuleBlock violate principle of