Wilkinson's polynomial


Solution P9.2.2

We can generate Wilkinson's polynomial from its roots using Polynomial.fromroots():

In [x]: import numpy as np
In [x]: Polynomial = np.polynomial.polynomial.Polynomial
In [x]: w = Polynomial.fromroots(range(1, 21))
In [x]: w.roots()
Out[x]: 
array([  1.00000000+0.j        ,   2.00000000+0.j        ,
         2.99999999+0.j        ,   4.00000021+0.j        ,
         4.99999682+0.j        ,   6.00003271+0.j        ,
         6.99974256+0.j        ,   8.00169121+0.j        ,
         8.99097171+0.j        ,  10.04201944+0.j        ,
        10.88097039+0.j        ,  12.41231512-0.30062093j,
        12.41231512+0.30062093j,  14.47999390-0.60517633j,
        14.47999390+0.60517633j,  16.59030707-0.43578888j,
        16.59030707+0.43578888j,  18.16420626+0.j        ,
        18.95026108+0.j        ,  20.00487543+0.j        ])

Note how the smaller roots are reproduced quite well (they are stable with respect to perturbation of the coefficients) but some of the larger ones are very poorly reproduced and become complex because of the finite precision of the representation of the polynomial's coefficients.