Re: Getting the (correct) roots of a polynomial

• To: mathgroup at smc.vnet.net
• Subject: [mg67274] Re: [mg67247] Getting the (correct) roots of a polynomial
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Thu, 15 Jun 2006 03:26:01 -0400 (EDT)
• References: <200606141029.GAA28180@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```nickc8514 at yahoo.com wrote:
> Hi, I have what seems like a simple question, but I haven't been able
> to find a satisfactory answer.  I wish to find the complex roots of the
> quartic polynomial
>
> (z - I u) z (z + I v) (z - I v) + (g^2) z (z + I v) + (g^2) z (z - I v)
> + (g^2)(z + I v)(z - I v)
>
> where z is the variable of the polynomial, u, v, and g are parameters
> whose values will be specified later, and I is the imaginary unit.  In
> the end, I wish to be able to fix the values of u and g and be able to
> examine (say, by plotting) how the values of the roots depend on v.
> I've tried several (possibly equivalent) approaches, but I seem to keep
> getting answers that are not actually roots.  I am using Mathematica
> 5.2 on SunOS.
>
> I assigned the above expression to the name Poly.  I then tried
> Solve[Poly==0,z].  If I then substitute in the particular values
> TestVals = {u -> 1/10,  v -> 10^4, g -> 1, w -> 10^7}, I find that the
> roots produced by Solve yield significantly non-zero values when
> plugged back into poly, meaning that the following does not produce a
> list of essentially zero values:
>
> Soln = Solve[Poly==0,z];
> Table[N[(Poly /. Soln[[j]]) /. TestVals],{j,1,4}]
>
> The values returned were of order 10^9.  I also tried plugging in the
> parameter values first (i.e. applying the rules in TestVals) before
> using Solve, and this also produced bad solutions.
>
> I also tried using the Root object, thinking that this might produce
> different results.  In this case, I converted Poly to a function of a
> single variable and applied the Root function to get each root.
>
> Table[Root[Evaluate[Poly /. z->#], j],{j,1,4}]
>
> Plugging these values back into Poly and plugging in the parameter
> values (again, by applying the TestVals rules), I get two values that
> are essentially zero (order 10^(-8) ) and two values that are
> significantly non-zero (order 1).  Again, I tried pluggin in the
> parameter values first and then using root, but this doesn't seem to
> change the result (which is what I expected).  After looking at the
> help for the Root function, I also tried setting the option
> ExactRootIsolation->True inside the Root function, but this did not
> seem to improve my results.
>
> I'd appreciate any guidance on how to get the correct roots to this
> polynomial.  I think I either simply don't understand enough about the
> way in which these functions (particularly root) work to understand the
> origin of the problem, or else it's some issue of numerics having to do
> with rounding errors or the like.  Again, I'd like to be able to
> produce the correct roots for any given values of the parameters, and
> particularly for a large range of values of v.  It's probably also
> worth noting that given the other parameter values used above, if v is
> set to be a small value (say, v = 5), then the solutions given by any
> of the methods above seem to be quite good.
>
> Many thanks,
>
> Nick

What you see is an artifact of using machine arithmetic. You need to
evaluate numerically at higher precision, using software significance
arithmetic. The following will show that the solutions are fine. First
form solutions and exact residuals.

poly = g^2*z*((-I)*v + z) + g^2*z*(I*v + z) +
g^2*((-I)*v + z)*(I*v + z) +
z*((-I)*u + z)*((-I)*v + z)*(I*v + z);
soln = Solve[poly==0, z];
testVals = {u->1/10, v->10^4, g->1, w->10^7};
residuals = poly /. soln /. testVals;

This shows the residuals are zero to around 49 places of accuracy.

In[16]:= N[residuals, 25] // InputForm

N::meprec: Internal precision limit \$MaxExtraPrecision = 50.
I    Sqrt[<<1>>]   Sqrt[<<1>>]
reached while evaluating {(-- - ----------- - -----------) <<1>> +
40        2             2
<<3>>, <<3>>}.

Out[16]//InputForm=
{0``48.87974661497277 + 0``48.87974661510448*I,
0``48.87975530086154 + 0``48.87975530099325*I,
0``48.87962522398726 + 0``48.87962522398726*I,
0``48.87962522529048 + 0``48.879625225290475*I}

Often it is preferable to avoid Cardano formulas and work with Root
objects. This might be done as below.

SetOptions[Roots, Cubics->False, Quartics->False];

InputForm[soln2 = Solve[poly==0, z]]

Out[23]//InputForm=
{{z -> Root[g^2*v^2 - I*u*v^2*#1 + 3*g^2*#1^2 + v^2*#1^2 - I*u*#1^3 +
#1^4 & , 1]},
{z -> Root[g^2*v^2 - I*u*v^2*#1 + 3*g^2*#1^2 + v^2*#1^2 - I*u*#1^3 +
#1^4 & , 2]},
{z -> Root[g^2*v^2 - I*u*v^2*#1 + 3*g^2*#1^2 + v^2*#1^2 - I*u*#1^3 +
#1^4 & , 3]},
{z -> Root[g^2*v^2 - I*u*v^2*#1 + 3*g^2*#1^2 + v^2*#1^2 - I*u*#1^3 +
#1^4 & , 4]}}

In[24]:= residuals2 = poly /. soln2 /. testVals;

Here we see that numerical evaluation of the residuals with machine
arithmetic still gives (somewhat smaller) nonzero values.

In[26]:= InputForm[N[residuals2]]

Out[26]//InputForm=
{0.31519216299057007 + 0.*I, -1.546086996793747 + 0.*I,
-3.0824594432488084*^-8 + 0.*I, 3.14321368932724*^-9 + 0.*I}

Repeating at higher precision using significance arithmetic shows that
they are in fact zero.

In[27]:= InputForm[N[residuals2, 25]]

N::meprec: Internal precision limit \$MaxExtraPrecision = 50.
reached while evaluating {Root[10<<16>>0 + <<4>> & , 1] <<1>> + <<3>>,
<<3>>}.

Out[27]//InputForm=
{0``19.864646654512317 + 0``1.980775265487835*I,
0``19.86465186604562 + 0``1.98078395137652*I,
0``66.25225198045372 + 0``66.5371330235564*I,
0``66.2641880947502 + 0``66.54840337518279*I}

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: Getting the (correct) roots of a polynomial
• Next by Date: Re: Re: Resolve/Reduce is taking forever
• Previous by thread: Getting the (correct) roots of a polynomial
• Next by thread: Re: Getting the (correct) roots of a polynomial