Re: FindRoot

*To*: mathgroup at smc.vnet.net*Subject*: [mg111117] Re: FindRoot*From*: Bill Rowe <readnews at sbcglobal.net>*Date*: Wed, 21 Jul 2010 07:11:08 -0400 (EDT)

On 7/20/10 at 3:43 AM, tmatsoukas at 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.