Wilkinson’s perfidious polynomial

\(w(x) = (x - 1)(x - 2) \cdots (x - 20)\) has well separated integer roots, yet James Wilkinson found that subtracting \(2^{-23}\) from its \(x^{19}\) coefficient, \(-210\), sends ten of the roots far into the complex plane. The root condition numbers explain why: the roots near 16 amplify a relative change in that coefficient by a factor of about \(3 \times 10^{10}\).

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.numerical_analysis import root_condition_numbers, wilkinson_polynomial

Roots before and after a tiny perturbation

w = wilkinson_polynomial(20)
perturbed = w.copy()
perturbed[1] -= 2.0**-23
roots = np.roots(perturbed)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2))
ax1.scatter(np.arange(1, 21), np.zeros(20), color="black", label="roots of $w$")
ax1.scatter(roots.real, roots.imag, color="firebrick", marker="x", label=r"$a_{19} \to a_{19} - 2^{-23}$")
ax1.set_xlabel("real part")
ax1.set_ylabel("imaginary part")
ax1.set_title("Wilkinson's perturbation experiment")
ax1.legend(fontsize=8)

<matplotlib.legend.Legend object at 0x1196e3b60>

Condition number of each root with respect to \(a_{19}\)

kappa = [root_condition_numbers(w, float(r))[1] for r in range(1, 21)]
ax2.semilogy(np.arange(1, 21), kappa, "o-")
ax2.set_xlabel("root $r$")
ax2.set_ylabel(r"$\kappa_{19}(r) = |a_{19}|\, r^{18} / |w'(r)|$")
ax2.set_title("Relative sensitivity of each root")
fig.tight_layout()

print("roots with |imag| > 0.5:", np.sort_complex(roots[np.abs(roots.imag) > 0.5]).round(3))
print(f"largest condition number: {max(kappa):.2e} at r = {1 + int(np.argmax(kappa))}")

plt.show()

roots with |imag| > 0.5: [10.095-0.64j  10.095+0.64j  11.793-1.651j 11.793+1.651j 13.992-2.519j
 13.992+2.519j 16.731-2.813j 16.731+2.813j 19.502-1.94j  19.502+1.94j ]
largest condition number: 3.16e+10 at r = 16