Message Boards Message Boards

GROUPS:

Visualizing simultaneous iteration methods for polynomial roots

Posted 1 year ago
2570 Views
|
3 Replies
|
9 Total Likes
|

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.

Szeg? curve

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]

Durand-Kerner

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]

Ehrlich-Aberth

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

POSTED BY: J. M.
Answer
3 Replies

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

Posted 1 year ago

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.

POSTED BY: J. M.
Answer

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract