Critical polynomial roots bigger than 2

Question

In the question Why study critical polynomials?, it was said that each hyperbolic component of the Mandelbrot set contains a root for a critical polynomial. That is, the sequence $$\{0, 0^2 + c, c^2 + c, (c^2 + c)^2 + c, \dots, p_n, p_n^2 + c, \dots\}$$ the sequence of critical polynomials has roots corresponding to one of the blobs/mini-mandelbrots in the Mandelbrot set.

However, I also know that the mandelbrot set is contained inside $\{z \in \mathbb{C} : |z|<2\}$. So we should have that the roots of each polynomial is bounded inside this set.

In Sage I typed

sage: R.<c> = PolynomialRing(ZZ)
sage: R
Univariate Polynomial Ring in c over Integer Ring
sage: p = 0
sage: P = []
sage: for i in xrange(10):
....:         p = p^2 + c
....:         P.append( points( (p.complex_roots()), hue=i/10, size = 20 ))
sage: sum(P)

I get roots with size greater than 2. What do these roots mean? Have I messed up somewhere in the code? Are the degrees so big that sage can't handle it (unlikely).

This is for the first 7 polynomials

This is all 10

Answer 1

This is just numerical error - a standard thing that you've got to get used to dealing with in this context. On my computer, the errors become apparent around 7 or 8 iterations. If you expand the polynomial at 8 iterations and check the largest coefficient, you should get 2676118542978972739644, which is larger than $2^{63}$ or the maximum machine integer. Mathematica and Sage can both deal with this but it's not automatic. You'll have to trigger extra precision somehow.

In Mathematica, you could do something like so:

f[c_][z_] = z^2 + c;
F[c][z_] = Nest[f[c], z, 8];
pts = c /. NSolve[F[c][0] == 0, c,
    WorkingPrecision -> 20];
ListPlot[{Re[#], Im[#]} & /@ pts,
 PlotRange -> All]

Notice the WorkingPrecision option. If you remove it, you'll find you get roots outside the disk of radius two. It quickly gets worse as the number of iterates increases.

In sage there is a way to do this too:

sage: R.<c> = PolynomialRing(ZZ)
sage: p = 0
sage: P = []
sage: for i in xrange(10):
....:         p = p^2 + c
....:         P.append( points((p.roots(ring=ComplexField(400), multiplicities=False)), hue=i/10, size = 20 ))
sage: sum(P)

Where we change ComplexField(400) to be how precise we want it. In this case, it would be 400 bits of precision.

Doing this gives us the great pictures like

And here's a live version with the iteration count reduced by one to get it to run a little faster.