Consider the problem of computing --- and, for example, plotting --- the eigenvalues of the matrix
If we wrote a Fortran program that accepted a numerical value for b and
then computed the eigenvalues of A by using, for example, some
standard library package for computing eigenvalues, we would consider
this a numerical approach. In Maple, we could use evalf(Eigenvals(A))
to numerically compute the eigenvalues of a numeric matrix A.
If instead we first found an `explicit formula' (whatever that means) for the eigenvalues, in terms of the parameter b, and then wrote a program to evaluate those eigenvalues using the formula, then we would consider this a symbolic approach.
The distinguishing feature between `numerics' and `symbolics' thus seems to be that numerical methods evaluate objects to numbers as soon as possible, while symbolic methods do so as late as possible (this was articulated by J. Rice at ISSAC '92). Deeper thought shows that the distinction is really artificial (a numeral is, after all, a symbol itself), and that there is a whole spectrum of methods available, depending on when quantities are evaluated.
Let us return to the matrix A. We can use Maple to compute explicit expressions for the eigenvalues by finding the characteristic polynomial and then finding its roots, as we were all taught to do in our first linear algebra class.
> with(linalg): Warning: new definition for norm Warning: new definition for trace > A := matrix([[ 0, 1], [1, 2*b]]);
> p := charpoly(A,lambda);
> answer := solve(p,lambda);
On the numerical end of the spectrum we can write (for example) a Matlab program to evaluate the eigenvalues given a numerical value of b:
function lambda = eiga(b)
A = [ 0 1
1 2*b];
lambda = eig(A);
We can use either the Maple output or the Matlab program for subsequent
analysis or plotting. For example, we can plot the eigenvalues for
e.g. 
using either representation.
Which approach is `better'? Surprisingly, neither is totally dominant. The symbolic approach really isn't explicit --- it involves taking a square root, an infinite algebraic operation, which we are really just `used to' --- we have come to accept finding square roots as `basic operations'. However, the computer implementation of `square root' is really quite complicated. ( Geometrically, of course, the idea is simple and even has a finite straightedge-and-compass construction --- hence, perhaps, our ingrained sense of `explicitness' can be rationalized).
But in Matlab, taking eigenvalues is a basic operation! Hence, in some sense, the Matlab program is also `explicit'! One just has to `get used to the idea'. So there is no clear philosophical superiority of one approach over the other.
Are there practical superiorities? Yes. The explicit formula is on the one hand more efficient than the numerical method, and on the other hand less stable than the numerical method (for large values of b).
The last point is not as well-known as it deserves to be. We illustrate
it by noting that the product of the two eigenvalues is identically
-1, for all values of b. When we plot
, (where each piece is calculated
as displayed), for large b we get the plot in Figure 2.
Figure 2: Numerical instability of the quadratic
formula.
We clearly see that something goes wrong for b between 
and
--- in fact, the explicit formula 
is
giving incorrect answers when evaluated (the rounding error made when
calculating 
is revealed by catastrophic cancellation of
the leading digits --- see [21] for details).
I emphasize that this shows symbolic expressions themselves can be unstable to evaluate --- and when this happens we don't care how fast they evaluate because they give the wrong answer.
In contrast, the Matlab program can be shown to be stable in the backward error sense: eig(A) will compute the exact eigenvalues of a matrix very near to A. Further, one can show that the eigenvalues of A are insensitive to small changes in A --- in numerical analysis parlance, the eigenvalues of A are well-conditioned [47].
Remark. Numerical experiments show that Matlab will in fact produce exact eigenvalues of a matrix of the form
where 
is less than the machine epsilon and 
is
usually small relative to b if b > 1 and small absolutely if
b < 1. This is an example of an experimental `structured backward
error result'. Structured backward error analysis is currently a
topic of much research interest. See [60] for an
introduction.
Can we, by judicious use of the art of scientific computing, create a program with the efficiency of the symbolic expression with the stability of the numerical eigenvalue calculation? The answer is yes, but the cost is the replacement of the explicit expressions for the eigenvalues with a `computation sequence' --- i.e. a small, `straight-line' program! Assuming that b is real, the following (well-known) computation sequence does the trick:
No catastrophic cancellation occurs here to reveal the rounding errors made. This computation sequence (in words, `compute the root with the largest magnitude first, since that can be done accurately in floating-point arithmetic, and then use that one to compute the small root') is efficient and stable. It is, if you like, a hybrid symbolic/numeric `recipe'.
But what about insight? Remember Hamming's famous motto, `The Purpose of
Computing is Insight, not Numbers'! Which formulas
give the most insight? Well, each contributes
something. The quadratic formula has a symmetry that the computation
sequence lacks, the computation sequence allows easier computation of
the asymptotics of 
as 
, and the eigenvalue
formulation retains important information on the sensitivity of the
problem to changes in the data. Depending on the problem context,
any one of the three formulations could give the required insight.
Of course, a good analyst would be aware of all three (and perhaps
more) and switch between representations as appropriate.
Final Remarks on this Example. Many ideas are discussed in this book using elementary examples. I will rarely even bother to say `the generalization of these ideas to larger examples is obvious'. In this case perhaps the generalizations do need emphasis.
The problem of finding eigenvalues of symmetric or unsymmetric n by
n matrices is well studied. The naive approach, the first learned
in a linear algebra course and hence the first thing one tries to
use in practice, is to first find the
characteristic polynomial and then to try
to find its roots. If n is less than 5, then there are explicit
formulas in terms of radicals, analogous to the quadratic formula
(see the excellent book [6]
for simple derivations, and see the exercises). These
formulas suffer from instability, just as the quadratic formula
does [67]. For n=5, there is an `explicit'
formula in terms of
elliptic functions [38] --- which can also be unstable.
There are no formulas in terms of radicals for 
, and hence
numerical methods must be used.
It has been known for more than thirty years that this naive approach of finding eigenvalues by first finding the characteristic polynomial is inherently unstable, and very bad in a practical sense (see the excellent expository essay by Wilkinson in [104]). Just because we can compute the characteristic polynomial doesn't mean that we should do so as a step to find the roots. Sadly, people still use this method (see e.g. [17] as an example of an otherwise adequate, even good, text that advocates this method). Perversely, part of the attraction of symbolic systems is the thought that this horrid approach can be used in practice, with very high precision used to find the roots of the characteristic polynomial to sufficient accuracy --- that is, one tries to `buy more accuracy by paying for more precision'. This fixes the problem, all right, but the cost can be horrendous. See the exercises.
Exercises
Find its eigenvalues by hand. Check your results with Maple.
A convenient way of entering this matrix in Maple is to use the following short procedure:
A := n -> linalg[toeplitz]([seq(a^k,k=0..n-1)]);See
?seq and ?toeplitz for a discussion of what these commands do
individually. Here, seq constructs the sequence of elements 
, and this is enclosed in list brackets
(square brackets) since toeplitz expects a list as input, from
which it constructs the Toeplitz matrix given by that list of
elements. Issuing the command
A3 := A(3);will reproduce the matrix of question 1, and assign it to the variable
A3 for later manipulation. Similarly A(5)
generates the 5-by-5 version of this matrix, etc.
charpoly and factor, find
the eigenvalues of 
. Is Maple's answer an `explicit' formula
or is it really a computation sequence? See ?label for a discussion
of how you can control the intermediate quantities %1, %2, etc.
.
Factor it. Given the `explicit' formula for fifth-degree
polynomials in terms of elliptic functions mentioned in the text, you
see that we could find explicit formulas for the roots of this
polynomial, but personally I doubt whether that would be useful.
is 
. Thus we expect some difficulty in numerically
computing
the roots of the characteristic polynomial for small a and large
n, because near-multiplicity of the roots makes them sensitive to
changes in a (or, indeed, roundoff error).
For 
, use realroot (see ?realroot for usage) to
isolate
all ten real roots of the characteristic polynomial of 
, to
an accuracy of 
.
(The roots are guaranteed to be all real because the matrix is
symmetric). Try to use fsolve or some other numerical rootfinder to
get the same results. Compare with direct computation of the
eigenvalues using evalf(Eigenvals(A10)). Repeat for
.
. Are the
eigenvalues sensitive to changes in a?