Visualizing simultaneous iteration methods for polynomial roots

Staff Picks
, , , ,
1

One of my longstanding research interests in numerical analysis has been the family of “simultaneous iteration” methods for finding polynomial roots. (See e.g. McNamee’s book for a comprehensive survey.) Briefly put, these are modified Newton-Raphson methods that allow one to find the roots of a polynomial all at once, as opposed to finding them one at a time and “deflating”.

I had the idea to visualize how these algorithms gradually proceed from initial approximations to the roots, up to convergence. After a number of experiments, I settled on using domain coloring for visualization. I have found that the logarithmic derivatives of polynomials gave particularly striking plots.

For this post, I have used the scaled exponential sum of degree 20:

\frac{20!}{20^{20}}\sum_{k=0}^{20}\frac{(20z)^k}{k!}

as the example polynomial whose roots we want to see. It is known that the zeroes of this polynomial asymptotically approach the so-called Szeg? curve as the polynomial degree goes to infinity.

vpliw.gif

expPoly[x_] = With[{n = 20}, Sum[(n! (n x)^k)/(k! n^n), {k, 0, n}]]

I will now look at two of the most popular simultaneous iteration methods. The first one is the (Weierstrass-)Durand-Kerner method,

x_i^{(k+1)}=x_i^{(k)}-\frac{p(x_i^{(k)})}{\prod\limits_{j\neq i} (x_i^{(k)}-x_j^{(k)})},\qquad i=1\dots n;\; k=0,1,\dots

which is (typically) quadratically convergent. (Note that in simultaneous iteration methods, the polynomials are always assumed to be monic (i.e., unit leading coefficient).)

Implementing the iteration is easy using FixedPointList[]. As is customary with these methods, we use as a starting approximation points equispaced around the unit circle, and slightly rotated:

ptsdk = FixedPointList[# - expPoly[#]/Table[Apply[Times, #[[k]] - Delete[#, k]],
                                            {k, Length[#]}] &,
                       N[Exp[2 ? I Range[0, 19]/20 - I ?/40]], 40,
                       SameTest -> (EuclideanDistance[##] < 1.*^-6 &)];

I use a loose convergence criterion that is good enough for visualization purposes.

For the domain coloring plot, I will use a slightly modified version of the DLMF color scheme, based on an idea of Quilez.

DLMFPhaseColor[u_, s_:1, b_:1] := Module[{rgb},
          rgb = Clip[{1, -1, -1} Abs[{8, 4, 8} Mod[u/(2 ?), 1] -
                     {9, 3, 11}/2] + {-3, 3, 5}/2, {0, 1}];
          rgb = (3 - 2 rgb) rgb^2;
          Apply[RGBColor, b (1 + s (rgb - 1))]]

I then use a simplified version of code originally written by user Heike on Mathematica Stack Exchange:

dcdk = RegionPlot[True, {x, -9/8, 9/8}, {y, -9/8, 9/8}, 
                  ColorFunction ->
                  Function[{x, y}, DLMFPhaseColor[Arg[Total[1/(x + I y - #)]]]], 
                  ColorFunctionScaling -> False, Frame -> False, 
                  PlotPoints -> 405] & /@ ptsdk;

(This takes some time, due to the high PlotPoints setting.)

We can now see an animation:

ListAnimate[dcdk]

dcdk.gif

The other method I will be looking at in this post is the (typically) cubically convergent Ehrlich-Aberth(-Maehly) method,

x_i^{(k+1)}=x_i^{(k)}-\frac{\tfrac{p(x_i^{(k)})}{p^\prime(x_i^{(k)})}}{1-\tfrac{p(x_i^{(k)})}{p^\prime(x_i^{(k)})}\sum\limits_{j\neq i} \tfrac1{x_i^{(k)}-x_j^{(k)}}},\qquad i=1\dots n;\; k=0,1,\dots

which is also one of the methods available in Mathematica’s NSolve[]/NRoots[].

Unfortunately, I have no way to get the iterates generated by NSolve[], so I had to reimplement the method myself. We can use essentially the same code as was used for Durand-Kerner, with a few changes:

ptsea = FixedPointList[With[{ld = expPoly[#]/expPoly'[#]},
                            # - ld/(1 - ld Table[Tr[1/(#[[k]] - Delete[#, k])],
                                                 {k, Length[#]}])] &,
                       N[Exp[2 ? I Range[0, 19]/20 - I ?/40]], 40,
                       SameTest -> (EuclideanDistance[##] < 1.*^-6 &)];

dcea = RegionPlot[True, {x, -9/8, 9/8}, {y, -9/8, 9/8}, 
                  ColorFunction ->
                  Function[{x, y}, DLMFPhaseColor[Arg[Total[1/(x + I y - #)]]]], 
                  ColorFunctionScaling -> False, Frame -> False, 
                  PlotPoints -> 405] & /@ ptsea;

ListAnimate[dcea]

dcea.gif

It would be interesting to use this visualization technique on other polynomials with interesting root structure, as well as other simultaneous iteration methods.

9 Likes
2

2015-11-02_10-07-11.png

3

Postscript: people using version 12 should now be able to use the new function ComplexPlot[] to visualize the iterations instead of the functions I wrote.

1 Like
4

Hi J.M.
Tanx for contributing your work.
I have learned a lot from it.
sasan