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.
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]
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]
It would be interesting to use this visualization technique on other polynomials with interesting root structure, as well as other simultaneous iteration methods.
