Computing Symbolic Eigenvalues and Eigenvectors in Maple

by Michael Monagan

    Abstract: How to compute exact symbolic values for eigenvalues and eigenvectors of both numeric and symbolic matrices, and computing with roots of polynomials.
    Application Areas/Subjects: Linear Algebra
    Keywords: Eigenvalues, eigenvectors, roots, RootOf, roots of polynomials
    See Also: MapleTech issue 9 1993

> restart:

Introduction

This is a partly a tutorial on how to compute symbolic eigenvalues and eigenvectors in Maple, and partly a collection of non-trivial examples to give an indication of the current limits of what can be done. The main difficulty in computing symbolically is in handling roots of polynomials. Maple has a facility called evala which computes with roots of polynomials represented using Maple's RootOf notation which often will yield better results, better in the sense that the results will be smaller, hence easier to understand, and also more efficiently computed. Some discussion of how to compute with these RootOf values is also given.

Using the Linalg Package

Let's begin by loading the linear algebra package, inputting a matrix, and doing some calculations

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> A := matrix([[1,2,3],[2,3,4],[4,5,6]]);

The eigenvalues of A are the roots of the characteristic polynomial which is the determinant of the matrix lambda I - A, i.e. a polynomial of degree n in the variable lambda. This polynomial is called the characteristic polynomial. One way of computing it in Maple is

> cp := det(lambda-A);

In this case, the characterisitc is a cubic polynomial with integer coefficients. We can use Maple's solve function to solve for the roots

> solve(cp,lambda);

Alternatively, we could have done this in one step using the builtin eigenvals command

> e := eigenvals(A);

Notice, that instead of getting decimal approximations for the eigenvalues, which is what a numerical piece of software would give, Maple has computed the exact algebraic values for the eigenvalues. You can get Maple to compute numerical eigenvalues and eigenvectors too if you wish. Maple will do this automatically for you if one or more of the inputs entries of the matrix is a decimal number. Let's just do this to see the difference

> B := map(evalf,A);

> eigenvals(B);

We see immediately one of the difficulties of working numerically. The eigenvalue 0 is not quite 0. In this tutorial we are mainly interested in computing exact formulae for the eigenvalues. Let's continue by computing the eigenvectors for A. They are the solutions to A x = lambda x, i.e. the nullspace of the matrix lambda I - A.

> nullspace(e[1]-A);

> nullspace(e[2]-A);

> nullspace(e[3]-A);

It is a little bit tricky to read because Maple does not put enough space between the entries in these vectors. Alternatively, one can use the builtin eigenvects command

> e := eigenvects(A,radical);

The output consists of a sequence of tripples of the form [e, m, b] where e is an eigenvalue, m it's multiplicity, and b a basis for the eigenspace for e, a set of vectors in Maple. We will explain the option `radical' in a moment. Let us now show that the results are correct, and in doing so, show how the user can access each of the eigenvectors. If x[i] is the eigenvector for the eigenvalue lambda[i], then x[i] satisfies A x = lambda x. Thus (lambda I - A) x = 0. The matrix lambda I - x is called the characteristic matrix. One way to construct it is to use the builtin function charmat(A,lambda). We'll use Maple directly. We pick off the first eigenvalue and it's corresponding eigenvector.

> e1 := e[1][1]; v1 := e[1][3][1];

> evalm( (e1 - A) &* v1 );

So far, so good. Infact, it all works for matrices containing formulae as well as numbers. For example, here is a matrix with 3 parameters, a,b, and c

> T := toeplitz( [a,b,c] );

> eigenvals(T,radical);

> e := eigenvects(T,radical);

A little bit more complicated, but actually, very nice. Let's just check the first eigenvector again.

Things however do get complicated when the characteristic polynomial does not factor into linear or quadratic factors. Because then the formulae for the roots of a cubic polynomial in general are messy, and a quartic polynomial are even worse. And of course, in general, there is no way to express the roots of a quintic or higher degree polynomial in terms of exact radicals. Let's consider a slight change, a purtubation to our A matrix

> T[2,1] := T[2,1]+epsilon;

> cp := det(lambda-T);

> factor(cp);

We see that the characteristic polynomial does not factor, and consequently, the roots for this cubic polynomial computed by Maple below are ``messy'', and of limited use.

> solve(cp,lambda);

Now the eigenvectors will be vectors whose entries are functions of these messy eigenvalues. Infact, you will probably not be surprised to hear that Maple will have difficulty in computing with these nested radicals. Infact, Maple will have difficulty in simply showing that these eigenvalues are roots of the characteristic polynomial. Computing with nested radicals is a known difficult computational problem. Now in order to compute the eigenvectors, we are are going to switch to a new notation if you like, for the roots of the characteristic polynomial. Instead of solving them in terms of exact radicals, and then proceeding to compute the eigenvectors in terms of these radicals, we will introduce an new variable alpha to denote the roots of the characteristic polynomial. In Maple

> alias( alpha=RootOf(cp,lambda) );

In one sense, we are cheating, we are just doing what every mathematician does whenever he gets an expression which is too big to work with, namely, he gives it a label. But infact, we are not cheating, because Maple knows how to compute with alpha, i.e. how to add, subtract, multiply and divide formulae involving alpha. It knows that alpha is a root of the characteristic polynomial.

> f := 1/(alpha);

> g := evala(f);

> evala(1/g);

Now when we compute the eigenvectors, Maple is going to compute them as a function of alpha. The output will be much simpler than it would be if we used the nested radical representation for the roots of the characteristic polynomial. The reader is encouraged to try instead eigenvects(A,radical); to force Maple to compute the eigenvectors in terms of exact radicals to see for himself how much simpler they are.

> eigenvects(T);

[This expression is several screens long - we recommend downloading the worksheet and viewing in Maple]

Now, it would appear that there is only one eigenvalue here, alpha. But it really stands for the three eigenvalues, the roots of the cubic characteristic polynomial. So if the user wants the three eigenvectors explicitly, what needs to be done is to solve the cubic polynomial, and substitute alpha for each of the three solutions. Let us now give an example which shows the limits of this symbolic approach, in the sense that this is a very nice result, but if the problem was slighly larger, a nice result could not be given. We are given the following matrix

> A := matrix([[-a1/a2, 1/a2, 0, 1/(R*a2)],[1, 0, 0, 0],
[-q1*a1/a2, q1/a2, 0, a2 + q1/(a2*R)],[q2, 0, 1, a1]] );

> cp := det(lambda-A);

> factor(cp);

So the characteristic polynomial does not factor. We are expecting very complicated roots for this polynomial if we try to solve for lambda. However, if you inspect the polynomial, you will notice that it is symmetric, i.e. coeff(cp,lambda,4) = coeff(cp,lambda,0) and coeff(cp,lambda,3) = coeff(cp,lambda,1). So solve is able to use a rational decomposition to find the roots in a relatively simple form.

> solve(cp,lambda);

The eigenvectors also turn out to be not too complicated when we allow Maple to use it's RootOf notation. The eigenvects command uses the RootOf notation by default.

> eigenvects(A);

Notice that the output is in terms of a RootOf. It is a polynomial of degree 4 in the variable _Z. We could substitute this for any of the 4 solutions that we computed above if we wanted to. Here is another worked example, of a matrix where the characteristic polynomial factors into three factors.

> B := matrix(6,6,
[0,0,0,0,0,a,
0,1,0,0,0,-b,
0,0,0,-a,b,0,
0,0,-c,0,0,0,
0,0,d,0,0,0,
c,-d,0,2,0,0] );

> cp := charpoly(B,lambda);

> f := factor(cp);

From the factorization of the characteristic polynomial, we see that one of the eigenvalues is 0, that two of the eigenvalues are + - sqrt(ac + bd), but the three others will be the roots of a cubic polynomial, so rather messy. Proceeding this time one step at a time, lets compute the eigenvector corresponding to the eigenvalue 0 by computing the nullspace of the characteristic matrix lambda I - A

> C := charmat(B,0);

> nullspace(C);

Now for the quadratic factor, here is the eigenvector for + sqrt(ac+bd)

> C := charmat(B, sqrt(a*c+b*d));

> nullspace(C);

Now for the cubic factor. Let's not use the radical notation here. Let's use the RootOf notation. First, let's give the roots of the cubic polynomial a label, beta

> alias( beta=RootOf(op(3,f),lambda) );

> C := charmat(B,beta);

> nullspace(C);

Finally, lets automate the entire procedure using the eigenvects command.

> eigenvects(B);

[This expression is several screens long - we recommend downloading the worksheet and viewing in Maple]

We can see that the eigenvector for beta is different from the one we computed previously. What has happened is that Maple has made a different entry equal to 1, the second entry, by dividing through. But the eigenvectors do agree up to multplication by a constant. The roots of the quadratic polynomial are displayed using a RootOf because we have not given an alias for this RootOf.

The limitations of Maple

We wish now to say something about the limitations of Maple. Will it work for complex matrices? Will it work for matrices which have algebraic numbers and functions in them, e.g. sqrt(2) or sqrt(1+x)? Yes it will. And the main facility that makes this possible is Maples factor command which can factor polynomials over algebraic number fields and, in Release 2, algebraic function fields. This example will show what we mean. Consider this matrix:

> T := toeplitz( [sqrt(2),x,omega] );

> cp := det(lambda-T);

The eigenvalues are of course the roots of this polynomial. We could use the general formula for a cubic polynomial to find them, but, it turns out that the polynomial factors. The factorization needs to be done over the algebraic number field with the sqrt(2) in it. We get

> factor(cp);

> eigenvals(T);

Some more examples with computing over algebraic number fields. Suppose we have the polynomial below and we try to factor it.

> p := x^4-4;

> factor(p);

Observe that the polynomial has not been factored into linear factors. The factor command factors a polynomial over the field implied by the coefficients in the input. In this example, the coefficients are integers. So factor finds all irreducible factors with integer coefficients. To factor the polynomial into linear factors, we have to introduce the sqrt(2) and the complex unit sqrt(-1) which is I in Maple.

> factor(p,sqrt(2));

> factor(p,I);

> factor(p,{sqrt(2),I});

Now this ability to factor polynomials over different number fields is exploited in the eigenvalues and eigenvectors commands to simplify the eigenvalues. The facilities in Maple also extend to algebraic functions, for example, sqrt(y). Consider this example

> p := x^4-4*y;

> factor(p);

> factor(p,sqrt(y));

Now here is an example where of a matrix where this factorization capability helps.

> T := toeplitz([sqrt(1+x), sqrt(2), 1+x]);

> cp := charpoly(T,lambda);

> factor(cp,sqrt(1+x));

> eigenvects(T);

Finally, an example from physics, also involving square roots.

> s := (2*t^2+2*t)^(1/2):
P := matrix(3,3,[t*(t+1), s, 0, s, t*(t+1)+2, s, 0, s, t*(t+1)]);

> cp := charpoly(P,lambda);

> factor(cp);

This example is rather interesting because the square roots cancelled out in the characteristic polynomial, which factors into 3 simple linear factors. So the eigenvalues are easy to compute. In computing the eigenvectors however, we will have to compute with these square roots

> eigenvects(P);

Now our examples have shown that Maple can work either using radicals or RootOf's, but actually, all the computations are done in the RootOf notation. I.e. the codes first convert from radicals too RootOfs. For example, the sqrt s in the example above converted into the RootOf notation is

> r := convert(s,RootOf);

You can see that indeed, s is a root of this polynomial in _Z. Converting back to the radical notation, we get

> convert(r,radical);

Now this square root can be factored into the form

> simplify(factor(s));

One could of course argue that this isn't really simpler because now we have two square roots. What we want to show now is that if we have an expression that involves both forms for s, Maple will complain and insist on using either one form, or the other, but not both at the same time. Here is what happens. Suppose we modify the 1,2 entry of the P matrix be

> P[1,2] := sqrt(2)*sqrt(t^2+t);

> print(P);

> eigenvects(P);

Error, (in evala) reducible RootOf
detected.  Substitutions are, {RootOf(_Z^2-t^2-t) =
-1/2*RootOf(_Z^2-2)*RootOf(_Z^2-2*t^2-2*t), RootOf(_Z^2-t^2-t) =
1/2*RootOf(_Z^2-2)*RootOf(_Z^2-2*t^2-2*t)}

We get an error. The first thing we notice is that Maple has converted the input radicals into the RootOf notation. What does the error mean. It means that Maple has discovered that the different RootOf's are not independent. Noting that sqrt(2) is RootOf(_Z^2-2) in the RootOf notation, Maple is telling us the relationship between them, namely that either.

> sqrt(t^2+t) = -sqrt(2)/2*s;

> sqrt(t^2+t) = +sqrt(2)/2*s;

Of course, it would be nice if Maple would automatically do the right simplification for us. But it doesn't. And the reason is that in general, because there is more than one possibility, it doesn't know which one it should choose. The user has to make sure himself that the input matrix does not have any such ``algebraic'' dependencies. We wish to conclude this tutorial by pointing out one further advantage of working algebraically, i.e. by doing all arithmetic exactly, and one disadvantage. The advantage is that by doing all the arithmetic exactly, there are no roundoff errors. So for matrices that have repeated eigenvalues, the eigenvectors will be computed correctly. Here is an example

> A := matrix(4,4,[0,1,0,0,-u^2,2*u,0,0,-u*v,v,u,0,-u*t,t,0,u]);

> cp := det(lambda-A);

> factor(cp);

> eigenvals(A);

> e := eigenvects(A);

This example shows a matrix in which the eigenvalue u is repeated 4 times. The eigenspace consists of a basis of more than one vector, 3 vectors. One normally expects to get 4 vectors in the basis. But it is possible for two vectors to conincide. We end by showing that the three eigenvectors computed are indeed correct. Instead of doing this one vector at a time, we build the matrix X of eigenvectors

> X := augment( e[3][1], e[3][2], e[3][3] );

> evalm( (u-A) &* X );

>