How to find and express all polynomial roots with greatest precision?

Question

I've already plowed through the documentation and a bunch of other Questions on this, but I've only found partial or tangentially related answers.

I need numerical results for all of the complex roots of high-order polynomials with real coefficients. But Solve[] give really bad results, for example:

z[n_, c_] := If[n > 1, z[n - 1, c]^2 + c, c];
expr = z[8, c] // Expand; (*Example polynomial*)
numericRoots = Solve[expr == 0, c] // N;
expr /. numericRoots // Abs // Max

The result should be zero, but the numerical solutions numericRoots are so far off that the result is:

(* 1.5218*10^34 *)

Because of the high orders, I need to MMa to return the most precise results possible and express them as precisely as possible. I thought this was the way:

r = SetPrecision[NSolve[1 + 2*c^2 + 3*c^3 + 3*c^4 + 3*c^5 + c^6 == 0, c, WorkingPrecision -> Infinity], Infinity]

But that doesn't actually solve the polynomial for some reason:

(* {{c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 1]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 2]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 3]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 4]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 5]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 6]}} *)

I can get numerical results by doing this:

c /. r // N

But I'm not sure that all those high-precision instructions I made in the NSolve[] command above actually made it into these numbers. This isn't a big issue for 6th order polynomials, but it becomes one for 1000th order or so.

So, why didn't my NSolve[] command solve what I wanted it to solve? And how can I get the highest-precision numerical roots I want?

The code in my first example (which is supposed to return 0) is a good test of whether or not a proposed solutions is working right.

Thanks!

Answer 1

Your numericRoots are OK. The core of your problem is that you use Expand to put your polynomial into power series form. High order power series are numerically unstable. The form your z function yields by itself is much better:

z[8, c] /. numericRoots // Abs // Max

(* 9.87841*10^-11 *)

Expanding on my answer here:

"I can see that what you're saying is correct, but it leaves me even more confused than before. z[8,c] and z[8,c]//Expand are exactly the same because the coefficients are integers."

But once you apply one of your numericRoots rules, the expressions are not exactly the same. Your rules substitute "machine numbers" for c. When you do mixed arithmetic with integers and machine numbers, you wind up with machine numbers. Machine arithmetic is sensitive to the order of operations. Mathematically equivalent expressions that are different in form will generally yield different results in machine arithmetic.

terms = expr /. Plus -> List

(* large expression that looks like an alluvial fan *)

You can see that while the coefficients are indeed integers, some are very large, too large for accurate representation as machine numbers. Then, even though the roots are modest, you're raising them to high powers, so the numbers you're adding are very large indeed:

terms /. numericRoots[[9]] // Abs // Max

(* 3.02339*10^49 *)

When you're adding a bunch of numbers of this magnitude in machine precision (~16 decimal digits), a residual of 10^34 is about as close to zero as you can reasonably expect.

"if simply changing the polynomial to a different but equivalent form creates a residual that humongous, how can I have any confidence in the root Mathematica gave me?"

Computation using machine numbers is fast, but (like other programming systems) Mathematica does not track the resulting precision. Sometimes, as in your formulation, the result contains no significant digits. There are, however, ways you can check.

One way is to use exact roots:

roots = Solve[expr == 0, c];

Don't look at them: they are mostly "algebraic numbers" expressed as complicated Root objects. But they are true, precise numbers. They are expressed in notation that is unfamiliar to most, but there is no general way to express algebraic numbers in familiar notation. Be happy that Mathematica can do this: in some situations it's very powerful. Mathematica is a bit lazier about computing with general algebraic numbers than with other kinds of numbers, but it can here if you coax it, and wait a few minutes:

Simplify[expr /. roots] // Abs // Max

(* 0 *)

So Mathematica asserts that its results are perfect ;-)

Another approach is to use arbitrary precision instead of machine precision:

numericRoots = N[Solve[expr == 0, c], 50];
expr /. numericRoots // Abs // Max // FullForm

(* 0``-2.568868432106492 *)

In arbitrary precision arithmetic, Mathematica tracks the precision of the result. Here, the result is 0, but with an accuracy of -2.56... decimal digits, so it might be 300 or so. That's a lot better than 10^34, and Mathematica informs you of the accuracy. You can do much better than a precision of 50 for your roots, of course:

numericRoots = N[Solve[expr == 0, c], 1000];
expr /. numericRoots // Abs // Max // FullForm

0``947.4311315678933

The residuals are zero to better than 947 decimal places.