FindRoot
in [Mathematica]
|
Prev: What's wrong with this assuming?
Next: Scoping constructs Block, Module, ModuleBlock violate principle of
From: Themis Matsoukas on 20 Jul 2010 03:43 I have a cubic polynomial in z (ZEQ) that depends on a parameter p. The polynomial has three real roots when p is in the interval (p1, p2). I then define function F[X] to be the ratio of the Max-to-Min (real) roots of ZEQ when p has the value X. Of course, F[X] makes sense only when ZEQ has three real roots, which restricts p to be in the interval p1<p<p2. The problem: while F[X] is evaluated (and plotted) correctly, if I try to solve for, say F[X]=6, I get a wrong answer (the correct answer is X=0.4). I suspect the problem is that Mathematica does not restrict the evaluation of F[X] within the interval (p1,p2). If I don't want to write my own routine for the root, is there a way to restrict FindRoot to work only inside the interval (p1, p2)? ZEQ=-0.049973 p^2+(0.44137p-0.00873 p^2) z-z^2+z^3; {p1,p2}={0.133652,0.690323}; F[X_]:=( z123=Sort[z/.Solve[(ZEQ/.p->X)==0,z]]; z123[[3]]/z123[[2]] ); FindRoot[F[X]-6, {X, 0.40}] F[X/.%]-6 {X->0.586195-5.08604*10^-16 I} -3.40038 TableForm[Table[{X, F[X]-6}, {X, 0.38, 0.42, 0.01}]] 0.38 0.55303 0.39 0.266973 0.40 -0.00436745 0.41 -0.262155 0.42 -0.50744 Thanks Themis
From: Bill Rowe on 21 Jul 2010 07:10 On 7/20/10 at 3:43 AM, tmatsoukas(a)me.com (Themis Matsoukas) wrote: >I have a cubic polynomial in z (ZEQ) that depends on a >parameter p. The polynomial has three real roots when p is in >the interval (p1, >p2). I then define function F[X] to be the ratio of the Max-to-Min >(real) roots of ZEQ when p has the value X. Of course, F[X] makes >sense only when ZEQ has three real roots, which restricts p to be in >the interval p1<p<p2. >The problem: while F[X] is evaluated (and plotted) correctly, if I >try to solve for, say F[X]=6, I get a wrong answer (the correct >answer is X=0.4). I suspect the problem is that Mathematica does not >restrict the evaluation of F[X] within the interval (p1,p2). If I >don't want to write my own routine for the root, is there a way to >restrict FindRoot to work only inside the interval (p1, p2)? >ZEQ=-0.049973 p^2+(0.44137p-0.00873 p^2) z-z^2+z^3; >{p1,p2}={0.133652,0.690323}; >F[X_]:=( >z123=Sort[z/.Solve[(ZEQ/.p->X)==0,z]]; >z123[[3]]/z123[[2]] >); >FindRoot[F[X]-6, {X, 0.40}] >F[X/.%]-6 >{X->0.586195-5.08604*10^-16 I} >-3.40038 Yes, it is possible to restrict the range over which FindRoot works. By doing FindRoot[F[x]-6,{x,0.4,p1,p2}] But, restricting the range for x doesn't seem to be enough to allow FindRoot to get the desired solution. The problem seems to be your ZEQ is numerically unstable. Note the difference between In[18]:= F[.4] - 6 Out[18]= -0.00436745 and In[19]:= F[x] - 6 /. x -> .4 Out[19]= 4.34196-1.48672*10^-14 I With this last, Solve will generate an error message indicating it was unable to find a solution with inexact coefficients. Additionally, if you do F[x] and look at the coefficients you will see there are approximate numbers with both quite large and quite small absolute values. That is, you will get different results due to issues with machine precision arithmetic when making substitutions for the variables in different ways. Finally, the plot Plot[F[x] - 6, {x, .399, .401}] indicates the root is somewhat less than 0.4. But given the apparent numerical instability for this particular problem, I am not certain this result is more accurate than another result.
From: Themis Matsoukas on 22 Jul 2010 05:45
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 |