Wilkinson's polynomial is defined as $$ W(x) = P_{N}(x) = \prod_{m=1}^{N}(x-m) $$ where $ N = 20 $

It is defined by the $(N+1) \;$ - dimensional array of coefficients $ p $ : $$ $$ $$ P_{N}(x) = p[0] x^{N} + p[1] x^{N-1} + \cdots + p[N] $$

In [1]:

import numpy as np

Initial value $P_1 (x) = (x - 1)$

In [3]:

p = np.array([1,-1])

Recursion: $p_n = [ p_{n-1},0 ] - [ 0, n \cdot p_{n-1} ]$

In [4]:

# arrays are zero-based
for n in range(2, N+1):
    p = np.append(p,[0]) - np.append([0], p*n)

In [5]:

r = np.roots(p)
print r

[ 10.   9.   8.   7.   6.   5.   4.   3.   2.   1.]

Wilkinson's polynomial $( N = 20 )$

The same with $ N = 20 $

In [6]:

N = 20
p = np.array([1,-1])
for n in range(2, N+1):
    p = np.append(p,[0]) - np.append([0], p*n)

In [7]:

x = np.arange(1,21)
print x

[ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20]

In [8]:

z = np.polyval(p,x)
print z

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

In [9]:

 r = np.roots(p)
print r

[ 87.00096132 +0.j          58.45406493+52.5151j      58.45406493-52.5151j
   2.58007091+54.91471744j   2.58007091-54.91471744j
   1.18198892 +0.28129905j   1.18198892 -0.28129905j
   0.90739576 +0.47880013j   0.90739576 -0.47880013j
   0.54219875 +0.8268119j    0.54219875 -0.8268119j
   0.11158039 +1.07480767j   0.11158039 -1.07480767j
  -0.76831677 +0.9416671j   -0.76831677 -0.9416671j
  -0.28453051 +0.86520925j  -0.28453051 -0.86520925j
  -0.86433849 +0.2809938j   -0.86433849 -0.2809938j   -0.72118908 +0.j        ]

Double precision (float - shorthand for float64)

In [10]:

p = np.array([1.0,-1.0],dtype=float)
for n in range(2, 21):
    p = np.append(p,[0.0]) - np.append([0.0], p*n)

In [12]:

x = np.arange(1,21,dtype=float)
print x

[  1.   2.   3.   4.   5.   6.   7.   8.   9.  10.  11.  12.  13.  14.  15.
  16.  17.  18.  19.  20.]

In [13]:

z = np.polyval(p,x)
print z

[  1.02400000e+03   8.19200000e+03   5.52960000e+04  -3.60448000e+05
  -1.53600000e+06  -9.51091200e+06  -2.10739200e+07  -5.47880960e+07
  -1.23918336e+08  -8.90880000e+07  -5.45177600e+08  -8.53770240e+08
  -1.35883571e+09  -1.93599078e+09  -3.76012800e+09  -5.96220314e+09
  -9.53860915e+09  -2.35116380e+10  -1.60700334e+10  -2.70295040e+10]

In [18]:

xx = np.arange(0.5, 20.5, 0.05, dtype=float)
yy = np.zeros_like(xx)
for i in range(xx.size):
    y = np.polyval(p,xx[i])
    yy[i] = sign(y)*log(1.0 + abs(y))

In [19]:

zz = np.zeros_like(xx)
plot(xx,yy, xx,zz)
title('Log of Wilkinson''s polynomial')
xlabel('$x$')

Out[19]:

[<matplotlib.lines.Line2D at 0x59f1c90>,
 <matplotlib.lines.Line2D at 0x59f1d90>]

In [20]:

r = np.roots(p)
print r

[ 19.99949308  19.00457196  17.98069544  17.04459974  15.92526766
  15.07956153  13.94558757  13.02424289  11.9984794   10.99484196
  10.00386198   8.99845889   8.00039949   6.99993101   6.00000798
   4.99999937   4.00000004   3.           2.           1.        ]

Now we modify the coefficient at $ x^{19}$

In [21]:

p[1] = p[1] - 2**(-23)

In [22]:

r = np.roots(p)
print r

[ 20.84689997+0.j          19.50242473+1.94032619j  19.50242473-1.94032619j
  16.73068765+2.8126196j   16.73068765-2.8126196j   13.99217151+2.51875869j
  13.99217151-2.51875869j   8.91849786+0.j          11.79322082+1.65176513j
  11.79322082-1.65176513j  10.09543537+0.64149236j  10.09543537-0.64149236j
   8.00699526+0.j           6.99972059+0.j           6.00000640+0.j
   4.99999986+0.j           4.00000000+0.j           3.00000000+0.j
   2.00000000+0.j           1.00000000+0.j        ]

In [25]:

xx = np.arange(0.5, 21.5, 0.05, dtype=float)
yy = np.zeros_like(xx)
for i in range(xx.size):
    y = np.polyval(p,xx[i])
    yy[i] = sign(y)*log(1.0 + abs(y))

In [26]:

zz = np.zeros_like(xx)
plot(xx,yy, xx,zz)
title('Modified Wilkinson''s polynomial')
xlabel('$x$')

Out[26]:

<matplotlib.text.Text at 0x5eb29f0>

Wilkinson's polynomial experiments

$N = 10$

Wilkinson's polynomial $( N = 20 )$