Wilkinson.nb

Wilkinson

Definition

P (x) = Underoverscript[&prod;, m = 1, arg3] (x - m)

Factored form

p20[x_] = Product[(x - i), {i, 20}]

(-20 + x) (-19 + x) (-18 + x) (-17 + x) (-16 + x) (-15 + x) (-14 + x) (-13 + x) (-12 + x) (-11 + x) (-10 + x) (-9 + x) (-8 + x) (-7 + x) (-6 + x) (-5 + x) (-4 + x) (-3 + x) (-2 + x) (-1 + x)

Expanded form

w20[x_] = Expand[Product[(x - i), {i, 20}]]

Obviously p20[x] = w20[x] .

Horner form

Later we will compare the above forms with the Horner form .

&lt;&lt;Algebra`Horner`

h20[x_] = Horner[w20[x], x]

Checking the roots

Integer roots

The roots are also integer :

r20i = Table[{x&rarr;i}, {i, 20}]

Applying the rules to the factored form

p20[x]/.r20i

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

and to the expanded form

w20[x]/.r20i

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

we got the same results .

Floating&nbsp;&nbsp;point roots

Define similar rules in the real (floating point) representation .

r20 = Table[{x&rarr;i}, {i, 1., 20., 1.}]

Applying these rules to p20[x] and w20[x] gives quite different results :

p20[x]/.r20

{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

w20[x]/.r20

This is the effect of round - off errors .

Horner form :

h20[x]/.r20

Numerical root finding

The rules r20 or r20i represent the exact roots of the Wilkinson

r = NSolve[w20[x], x]

Or in the list form :

x/.r

{1., 2., 3., 4., 5., 6., 7.00002, 7.99987, 9.00072, 9.99709, 11.0089, 11.9801, 13.0357, 13.9537, 15.0462, 15.965, 17.018, 17.9933, 19.0014, 19.9999}

With higher accuracy

NumberForm[%, {20, 17}]

The rules r applied to w20[x] give

w20[x]/.r

- still not zeroes .

Try the Horner form

rh = NSolve[h20[x], x]

Sensitivity experiments

pm[x_] = w20[x] - 2^(-23) x^19

The solution is

rm = NSolve[pm[x], x]

x/.rm

Substitut the roots back to the modified polinomial

pm[x]/.rm

Plot[p[x], {x, 0.9, 20.9}]

[Graphics:HTMLFiles/Wilkinson_71.gif]

-Graphics -

q[x_] = Sign[x] Log[1 + Abs[x]]

Log[1 + Abs[x]] Sign[x]

Plot[q[w20[x]], {x, 0.9, 20.1}]

[Graphics:HTMLFiles/Wilkinson_76.gif]

-Graphics -

Plot[q[p[x]], {x, 0.9, 21}]

[Graphics:HTMLFiles/Wilkinson_79.gif]

-Graphics -