Mathematics with Maple

S.Duzhin

Lesson 7: Matrices and Determinants

Mathematics: matrices, determinants, vectors, eigenvalues.

Maple: "linalg" package.

Today, we will review some basic notions of linear algebra with the help of the Maple package "linalg".

We begin by loading this package:

> with(linalg);

Warning, new definition for cond

Warning, new definition for norm

Warning, new definition for trace

In this word list, you can notice many words that you studied in the course of "Linear Algebra" last year. Here is a list of several most important words (form in parentheses are Maple abbreviations of the mathematical terms). Try to recall their meaning:

vector

basis

matrix

determinant (det)

minor

rank

eigenvalues (eigenvals)

eigenvectors (eigenvects)

In this lesson, we will basically study operations with matrices.

First of all, let us define several matrices, using the function "matrix":

> A1:=matrix([[1,0],[0,1],[a,b]]);

> A2:=matrix([[-1,x],[0,2]]);

You see that elements of matrices can be numbers and/or any expressions that contain letters.

We have two matrices. Let us try to do some operations with them.

> add(A1,A2);

Error, wrong number (or type) of parameters in function add

No go! You cannot add a 3x2 matrix to a 2x2 matrix. This operation is undefined.

Let us try multiplication:

> multiply(A1,A2);

Good, this works!

Let us try to multiply the same matrices in another order:

> multiply(A2,A1);

Error, (in multiply) non matching dimensions for vector/matrix product

This multiplication is impossible.

The general rule is that if A is a k x l matrix, and B is an m x n matrix, then multiply(A,B) makes sense if l=m, while add(A,B) works only if both k=m, l=n.

Exercise 1. Try the operations "add" and "multiply" on matrices A1,A2,A3,A4,A5, where:

> A3:=matrix([[1,1/2],[-1/2,1]]);

> A4:=matrix([[1,0,0],[0,2,1],[0,0,3]]);

> A5:=matrix([[0,1,2],[1,2,3],[2,3,0],[3,0,1]]);

Write a complete list of all such operations which can be performed.

Let us study more closely the case of SQUARE matrices, i.e. matrices that have the same number n of rows and columns: #(rows)=#(columns).

Any two such matrices can be added together and multiplied, so that the set of all matrices of size n x n makes a RING. Another nice thing about square matrices is that they have determinants.

Let us make some computations in the case n=2.

> S:=matrix([[0,-1],[1,0]]);

> T:=matrix([[1,1],[0,1]]);

Let us compute the two products S*T and T*S.

> ST:=multiply(S,T);

> TS:=multiply(T,S);

You may notice that the two results are different: matrix multiplication, in general, is not commutative.

Let us do some more multiplications:

> SST:=multiply(S,ST);

> SSTT:=multiply(SST,T);

> TSSTT:=multiply(T,SSTT);

> STSSTT:=multiply(S,TSSTT);

> TSTSSTT:=multiply(T,STSSTT);

> TTSTSSTT:=multiply(T,TSTSSTT);

You can proceed in this way and compute any monomial in variables S and T, or, in other words, any word written in the alphabet {S,T}. You will obtain many different matrices, but all of them have a common property: their determinant is always equal to1.

> det(TSSTT);

> det(TSTSSTT);

> det(TTSTSSTT);

This can be explained as follows:

You know that if A and B are two square matrices of the same size, then

det(A*B) = det(A)*det(B).

Therefore, if we only check that det(S)=det(T)=1:

> det(S);

> det(T);

-- then all the products ST, STT, STTS, etc. will automatically have determinant 1.

Exercise 2.

Express the following matrices as products of S and T (do not forget to check that the determinant is 1):

> matrix([[0,1],[-1,0]]);

> matrix([[1,3],[0,1]]);

> matrix([[1,-1],[0,1]]);

> matrix([[-1,4],[2,-7]]);

> matrix([[2,3],[3,5]]);

Using Maple, you can also find the general formula for the determinant of an arbitrary square matrix. For 2x2 matrices:

> det(matrix([[a11,a12],[a21,a22]]));

For 3x3 matrices:

> det(matrix([[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]]));

Exercise 3:

Find the formula of the determinant of a general matrix 4x4.

The next thing we will discuss is characteristic polynomial.

First, let us choose a nice matrix:

> M:=matrix([[-1/3,2/3,2/3],[2/3,-1/3,2/3],[2/3,2/3,-1/3]]);

The characteristic polynomial is the determinant of the characteristic matrix x*E-M, where E is the unit matrix (1 on the diagonal, 0 elsewhere):

> Mchar:=charmat(M,x);

> p:= det(Mchar);

It is also possible to compute the characteristic polynomial directly, without the characteristic matrix:

> p:=charpoly(M,x);

One nice thing about the characteristic polynomial is that if you substitute the matrix itself instead of the variable x in this polynomial, you will always obtain 0. Thus, the matrix is a root of its characteristic equation!

> subs(x=M,p);

> evalm(%);

The last command "evalm" is "evaluate matrix", it forces all operations to be performed as matrix operations, e.g. M^3 = multiply(M,multiply(M,M)).

By the way, there is another polynomial, called "minimal polynomial of a matrix", which is the polynomial q of a _minimal_ degree such that the matrix is its root: q(M)=0.

> q:=minpoly(M,x);

> evalm(subs(x=M,q));

Exercise 4.

For the following matrix, find the characteristic and the minimal polynomials and check that the matrix itself is a root of either of them:

> matrix([[9, 7, -1, -5], [-22, -16, 2, 10], [18, 16, 0, -10], [-19, -13, 1, 7]]);

The main application of characteristic polynomials is their use in the problem of finding eigenvalues and eigenvectors of a matrix.

Actually, eigenvalues are nothing but the roots of the characteristic equation p(x)=0.

If x=a is one of such roots, then eigenvectors corresponding to this value a are such vectors v that M*v=a*v, where M*v is matrix product of a matrix and a vector, while a*v is multiplication of a vector by a number.

Let us see this on our example:

> p;

> solve(p=0,x);

You see that there is a simple root 1 and a multiple root -1 of multiplicity 2.

The same result could be achived by one command:

> eigenvals(M);

Now, the eigenvectors:

> eigenvects(M);

Let me explain the meaning of this line.

1, 1, {[1 1 1]} means: 1 is an eigenvalue of multiplicity 1 with eigenvector [1 1 1],

-1, 2, {[-1 0 1], [-1 1 0]} means: -1 is an eigenvalue of multiplicity 2 with 2 linearly independent eigenvectors [-1 0 1] and [-1 1 0].

We can check this result by multiplication:

> v1:=vector([1,1,1]);

> v2:=vector([-1, 0, 1]);

> v3:=vector([-1,1,0]);

> multiply(M,v1);

> multiply(M,v2);

> multiply(M,v3);

Exercise 5.

Find the eigenvalues and eigenvectors of the matrix

> matrix([[9, 7, -1, -5], [-22, -16, 2, 10], [18, 16, 0, -10], [-19, -13, 1, 7]]);

and check the result by multiplication, following the previous example.

In the set of all square matrices, there are several remarkable subsets.

Let us recall just three of them: symmetric matrices, skew-symmetric matrices and orthogonal matrices.

To check that a matrix is symmetric, you can do as follows:

> evalm(M);

> equal(M,transpose(M));

This means that our matrix is symmetric. (You can also see this directly with your eyes!)

Similarly, to chack whether a matrix A is skew-symmetric, the following test can be used:

> iszero(M+transpose(M));

A nice property of skew-symmetric matrices is that their determinant is always a complete square of some polynomial (called "Pfaffian"):

> a:='a';f:='f';

> M:=matrix([[0,a,b,c],[-a,0,d,e],[-b,-d,0,f],[-c,-e,-f,0]]);

> det(M);

> factor(%);

Finally, orthogonal matrices are very important because they describe the motion of rigid bodies: if you multiply all vectors in the space by one and the same orthogonal matrix, the lengths of vectors and the angles between them do not change.

An orthogonal matrix can be defined as a matrix whose columns (or rows) are vectors of unit length orthogonal to each other.

This property holds if and only if the product M*transpose(M) is the unit matrix E.

In Maple, there is a special Boolean function orthog to check this property:

> orthog(M);

Exercise 6.

Complete the following matrix definition with numbers 1/2 and -1/2 in such a way that you get an orthogonal matrix. Check the answer by applying the function "orthog":

matrix([[1/2,1/2,1/2,1/2],[1/2,...,...,...],[1/2,...,...,...],[1/2,...,...,...]]);

Exercise 7.

Let

> A:=matrix([[sqrt(3)/2,1/2],[-1/2,sqrt(3)/2]]);

Using the Maple construction evalm(A^n) with different natural numbers n, find the minimal number n such that A^n is the unit matrix E = matrix([[1,0],[0,1]]). Explain the result.

Exercise 8.

Find the determinant of the Vandermonde matrix of size 3 and 4, then factorize it using the Maple function "factor".

> V3 := vandermonde([x1,x2,x3]);

> V4 := vandermonde([x1,x2,x3,x4]);