The algorithms employed for `poly`

and `roots`

illustrate
an interesting aspect of the modern approach to eigenvalue computation. `poly(A)`

generates
the characteristic polynomial of `A`

, and `roots(poly(A))`

finds
the roots of that polynomial, which are the eigenvalues of `A`

.
But both `poly`

and `roots`

use `eig`

,
which is based on similarity transformations. The classical approach,
which characterizes eigenvalues as roots of the characteristic polynomial,
is actually reversed.

If `A`

is an `n`

-by-`n`

matrix, `poly(A)`

produces
the coefficients `p(1)`

through `p(n+1)`

,
with `p(1)`

`=`

`1`

,
in

$$\mathrm{det}\left(\lambda I-A\right)={p}_{1}{\lambda}^{n}+\dots +{p}_{n}\lambda +{p}_{n+1}\text{\hspace{0.17em}}.$$

The algorithm is

z = eig(A);
p = zeros(n+1,1);
p(1) = 1;
for j = 1:n
p(2:j+1) = p(2:j+1)-z(j)*p(1:j);
end

This recursion is derived by expanding the product,

$$(\lambda -{\lambda}_{1})(\lambda -{\lambda}_{2})\dots (\lambda -{\lambda}_{n})\text{\hspace{0.17em}}.$$

It is possible to prove that `poly(A)`

produces
the coefficients in the characteristic polynomial of a matrix within
roundoff error of `A`

. This is true even if the eigenvalues
of `A`

are badly conditioned. The traditional algorithms
for obtaining the characteristic polynomial do not use the eigenvalues,
and do not have such satisfactory numerical properties.