# Visualizing simultaneous iteration methods for polynomial roots

**URL:** <https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870>\
**Category:** Staff Picks\
**Tags:** equation-solving, graphics-and-visualization, numerical-computation, staff-pick, wolfram-repositories\
**Created:** [April 4, 2019, 4:38am UTC](https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870 "2019-04-04T04:38:26Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![pleasureoffiguring](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/pleasureoffiguring/32/18977_2.png) [@pleasureoffiguring](https://community.wolfram.com/u/pleasureoffiguring)\
**Post date:** [April 4, 2019, 4:38am UTC](https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870/1 "2019-04-04T04:38:26Z")

</div>

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](https://books.google.com/books?hl=en&id=4PMqxwG-eqQC&pg=PA67) 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](https://en.wikipedia.org/wiki/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](http://mathworld.wolfram.com/ExponentialSumFunction.html) 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](https://doi.org/10.1216/rmjm/1181072998) [known](https://doi.org/10.1007/BFb0087909) that the zeroes of this polynomial asymptotically approach the so-called Szeg? curve as the polynomial degree goes to infinity.

![vpliw.gif](https://global.discourse-cdn.com/wolfram/original/3X/2/d/2d3fa13b586c874d4f0073af263cbe5a9b91e352.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](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_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](https://dlmf.nist.gov/help/vrml/aboutcolor#S2), based on an idea of [Quilez](https://www.shadertoy.com/view/MsS3Wc).

```
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](https://mathematica.stackexchange.com/a/7293/) 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](https://global.discourse-cdn.com/wolfram/original/3X/2/0/20eb4efb70b6a1cade7857b0cedb19e4befa9d1d.gif)

The other method I will be looking at in this post is the (typically) cubically convergent [Ehrlich-Aberth(-Maehly) method](https://en.wikipedia.org/wiki/Aberth_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](https://global.discourse-cdn.com/wolfram/original/3X/1/6/161200d73313e3279a631482fe83ace051edbfc6.gif)

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

---

<div class="post-metadata">

**Author:** ![wolframcommunity](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/wolframcommunity/32/8335_2.png) [@wolframcommunity](https://community.wolfram.com/u/wolframcommunity)\
**Post date:** [April 4, 2019, 3:09pm UTC](https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870/2 "2019-04-04T15:09:03Z")

</div>

![2015-11-02_10-07-11.png](https://global.discourse-cdn.com/wolfram/original/3X/5/0/50e59c9ead62f729e7130d3ea1482555a75602b6.png)

- Congratulations! This post is now a [Staff Pick](http://community.wolfram.com/content?curTag=staff+picks) as distinguished by a badge on [your profile](http://community.wolfram.com/web/pleasureoffiguring)! Thank you, keep it coming, and consider contributing your work to the [The Notebook Archive](http://notebookarchive.org)!

---

<div class="post-metadata">

**Author:** ![pleasureoffiguring](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/pleasureoffiguring/32/18977_2.png) [@pleasureoffiguring](https://community.wolfram.com/u/pleasureoffiguring)\
**Post date:** [April 20, 2019, 12:16am UTC](https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870/3 "2019-04-20T00:16:37Z")

</div>

Postscript: people using version 12 should now be able to use the new function [ComplexPlot[]](http://reference.wolfram.com/language/ref/ComplexPlot.html) to visualize the iterations instead of the functions I wrote.

---

<div class="post-metadata">

**Author:** ![mostafasasan52](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/mostafasasan52/32/32212_2.png) [@mostafasasan52](https://community.wolfram.com/u/mostafasasan52)\
**Post date:** [May 14, 2019, 9:20pm UTC](https://community.wolfram.com/t/visualizing-simultaneous-iteration-methods-for-polynomial-roots/13870/4 "2019-05-14T21:20:00Z")

</div>

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