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.