Basic Linear Algebra with Mathematica
Part 2
Version 0.89
©2002 by Roland Steinbauer <roland.steinbauer@univie.ac.at>, <http://www.mat.univie.ac.at/~stein/>

This short notebook provides an introduction into using Mathematica in the context of some more advanced concepts of a first year's university course on linear algebra. It is designed as a second part to the notebook: Basic Linear Algebra with Mathematica.

Basic Operations on Matrices

Recall that a matrix in Mathematica is a list of lists, where each individual list represents a row and the matrix is represented as a list of these rows. Alternatively you may use the Input...Create Table/Matrix/Palette menu item to enter a matrix into the program. To start with let us define

A := {{1, 0}, {1, 1}}
MatrixForm[A]

B:=(({{1, 2, 3}, {1, 3, 2}, {2, 2, 4}}))

The basic matrix operations can be performed using the obvious commands given in the box below.

({{"\!\(\*
StyleBox[\"\\\"Transpose[M]\\\"\",\n\"MR\"]\)\!\(\* StyleBox[\"\\\"\\\"\",\n\"TI\"]\)\!\(\* StyleBox[\"\\\"\\\"\",\n\"MR\"]\)", "transpose "}, {"\!\(\*
StyleBox[\"\\\"Inverse[\\\"\",\n\"MR\"]\)\!\(\* StyleBox[\"M\",\n\"MR\"]\)\!\(\* StyleBox[\"]\",\n\"TI\"]\) ", matrix inverse}, {"\!\(\*
StyleBox[\"\\\"Det[\\\"\",\n\"MR\"]\)\!\(\* StyleBox[\"M\",\n\"MR\"]\)\!\(\* StyleBox[\"]\",\n\"MR\"]\)", "determinant "}, {Tr[M], "trace "}})

Observe, that the commands for calculating the determinat and the trace are not consistent with the "overall convention" to completely spell out mathematical expressions!
Now let's play a little with these commands

Transpose[A] // MatrixForm

Inverse[B] // MatrixForm

To check this we may use

B . % // MatrixForm

Now let's see how Mathematica reacts to obviously stupid things.

M:=(({{1, 0, 1, 2}, {3, 4, 5, 7}, {1, 2, 0, 0}}))
Inverse[M]

Det[A]
Det[B]

So both of our matrices are nonsingular---so if we rowreduce them we should end up with the unity matrix, right?

RowReduce[A] // MatrixForm
RowReduce[B] // MatrixForm

Now let's try this one

Det[(({{1, 1, 2}, {2, 1, 2}, {3, 2, 4}}))]

Clearly the  the third row is the sum of the first and second ones. So the nullspace should be one-dimensional, right?

NullSpace[(({{1, 1, 2}, {2, 1, 2}, {3, 2, 4}}))]

What about its inverse?

Inverse[(({{1, 1, 2}, {2, 1, 2}, {3, 2, 4}}))]

Finally we try the trace command

Tr[(({{1, 1, 2}, {2, 1, 2}, {3, 2, 4}}))]

Can you explain what is happening here?

A:=(1/(Sqrt[2]))(({{1, 0, 1}, {0, (Sqrt[2]), 0}, {-1, 0, 1}}))
A . Transpose[A] // 
      MatrixForm

Does this aid to your understanding?

Inverse[A] // MatrixForm

To finish this section we return to the Inverse command

A := {{a, b}, {c, d}}
MatrixForm[A]

B := Inverse[A]
MatrixForm[B]

A . B

This doesn't look too nice. To simplify the above matrix we use the Together command which puts terms in a sum over a common denominator, and cancels factors in the result.

Together[%] // MatrixForm

Clear[A, B, M]

Eigenvalues, Eigenvectors

The basic Mathematica-commands to calculate eigenvalues and eigenvectors of a matrix M are (what a surprise) Eigenvalues[M] and Eigenvectors[M] . More precisely
({{"\!\(\*StyleBox[\"\\\"Eigenvalues[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"m\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "a list of the eigenvalues of \!\(\*StyleBox[\"\\\"m\\\"\", \"TI\"]\) "}, {"\!\(\* StyleBox[\"\\\"Eigenvectors[\\\"\",\n\"MR\"]\)\!\(\* StyleBox[\"\\\"m\\\"\",\n\"TI\"]\)\!\(\* StyleBox[\"\\\"]\\\"\",\n\"MR\"]\) ", "a list of the eigenvectors of \!\(\*StyleBox[\"\\\"m\\\"\", \"TI\"]\) "}})
We start out with the first one of these commands.

A := {{1, 0}, {1, 1}}
MatrixForm[A]
B:=(({{1, 2, 3}, {0, 3, 2}, {2, 0, 4}}))

Eigenvalues[A]
Eigenvalues[B]

As you can see Eigenvalues[M] returns the eigenvalues in form of a list. If some eigenvalue is a multiple zero of the characteristic polynomial (i.e., if its algebraic multiplicity is bigger than one) it appears multiple in the list. Of course the numer of its appearences equals its algebraic multiplicity. 
There is of course Mathematica-command which gives you the charcteristic polynomial of a matrix, i.e.,
({{CharacteristicPolynomial[M, t]", \(computes\ the\ characteristic\ polynomial\ of\ the\ matrix\ M\ in\ the\ variable\ t\)}})
But you can also easily put things together using the commands of the previous chapter and the IdentityMatrix command, which does the obvious thing.

IdentityMatrix[4] // MatrixForm

pA[λ_] := 
      Det[A - λ IdentityMatrix[2]]
pA[λ]
pB[λ_] := 
      CharacteristicPolynomial[B, λ]
pB[t]

You can of course solve for the zeroes of these polynomials using the Solve command or factor the ploynomial using the Factor command.

Solve[pA[λ] == 0, λ]
Factor[pA[λ]]
Solve[pB[λ] == 0, λ]
Factor[pB[λ]]

By default the Factor command assumes that all coefficients in the polynomials they produce must involve only rational numbers. But by setting the option Extension you can extend the domain of coefficients that will be allowed; in this case you have to use (for more information see the Helpbrowser...Factor).

Factor[pB[λ], Extension -> Sqrt[2]]

What about complex eigenvalues?

A1 := {{1, 1}, {-3/2, 0}}
Eigenvalues[A1]

If you want to know how I did choose this matrix A1 then consult Exercise 141!  To factor A1's characteristic polynomial you can use the following rather clumsy command.

Factor[CharacteristicPolynomial[A1, t], GaussianIntegers -> True, 
      Extension -> Sqrt[11]]

To calculate the eigenvectors of a matrix you can use the Eigenvectors command mentiones above. However, in many cases it will be more convenient to use the following command which computes eigenvalues and eigenvectors in one step.
({{"\!\(\* StyleBox[\"\\\"Eigensystem[\\\"\",\n\"MR\"]\)\!\(\* StyleBox[\"\\\"M\\\"\",\n\"TI\"]\)\!\(\* StyleBox[\"\\\"]\\\"\",\n\"MR\"]\) ", ""gives a list of the form "\!\(\*
StyleBox[\"{\",\n\"InlineFormula\"]\)""eigenvalues","eigenvectors""\!\(\*
StyleBox[\"}\",\n\"InlineFormula\"]\)"   for the matrix  M}})

Eigenvalues[A1]
Eigenvectors[A1]
Eigensystem[A1]

In the last line, i.e. the output of the Eigensystem command, the first list gives the two (distinct) eigenvalues, the second list gives the eigenvectors. Since the eigenvalues are distinct there are two one-dimensional eigenspaces and hence the geometrical multiplicity of both eigenvalues matches their algebaric multiplicity (which is one). (Is A1 diagonizable?). 
Let's check that these results are indeed correct.

v = Eigenvectors[A1][[1]]
λ  = Eigenvalues[A1][[1]]

A1 . v == λ v

Hmm...That doesn't look too good. So we try it this way

Simplify[A1 . v] == Simplify[λ v]

Let's see what happens if things get more involved.

A := {{1, 1, 0}, {0, 1, 0}, {0, 0, 2}}
MatrixForm[A]
Eigensystem[A]

Here we have a double eigenvalue λ\(TraditionalForm\`\(\(\_1\)\(=\)\)\)1 (i.e., algebraic multiplicity 2) and an eigenvalue (Subscript[λ, 2])=2 with algebraic multiplicity 1. In the list of eigenvectors, i.e. the second list in the output, the first entry is the eigenvector (1,0,0) to the eigenvalue (Subscript[λ, 1]). The next entry is the zero-vector! This is how Mathematica tells you that there is no eigenvector to (Subscript[λ, 1]) which is linearly independent of (1,0,0), i.e., that the eigenspace of (Subscript[λ, 1])is one-dimensional, i.e. that the geometric multiplicity of (Subscript[λ, 1])is one. The third entry (0,0,1) then is a basis of the eigensapce of (Subscript[λ, 1.]) 
Lets check this by direct calculation. We first calculate the eigenspace of (Subscript[λ, 1])by solving the respective system of linear equations.

RowReduce[A - 1*IdentityMatrix[3]] // MatrixForm

This now tells us that the space of solutions of the linear system, i.e. the eigenspace to (Subscript[λ, 1])is one-dimensional (just one free parameter for (Subscript[x, 1])for which we miss a leading one) generated by the vector (1,0,0). We of course can solved this problem directly by using the NullSpace command.

Now for (Subscript[λ, 2]).

NullSpace[A - 2*IdentityMatrix[3]]

To sum up the Eigensystem command returns (as the first list) the eigenvalues of the respective matrix. Then a basis of the eigensapce of the frist eigenvalue. If the geometric multplicity of this eigenvector is less than its algebaric one there is a zero- vector added for each missing dimesion of the eigenspace.   
Let's do some more examples.

A:=(({{0, 3, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 5}}))
Eigensystem[A]

Here the eigenspace of the eigenvalue (Subscript[λ, 1]) = 0 is generated by the vector (1,0,0,0,0) hence has algebraic multiplicity 2 and geometric multiplicity 1. Next a basis of the eigenspace of (Subscript[λ, 2]) = 1 is given by (0,0,0,1,0), (0,0,1,0,0). Hence the algebaric multiplicity of (Subscript[λ, 2])coincides with its geometric one which is 2. Finally the eigenspace for the single eigenvalue (Subscript[λ, 3]) = 5 is spanned by (0,0,0,0,1).
In general the calculation of the einegvalues and eigenvectors can be pretty messy as the next example shows.

B:=(({{1, 0, 1}, {2, 4, -9}, {0, -6, 0}}))
Eigensystem[B]

In such cases might be advisable to calculate the eigensystem numerically as follows...

Eigensystem[N[B]]

Clear[A, B, A1]

Diagonalization, Jordan-Decomposition

To decide wheter a given square matrix is diagonalizable or not and to calculate its diagonalization we can use the tools introduced above. Let's do an example.

A1:=(({{1, 2, 3}, {0, 2, 0}, {4, 1, 2}}))
Eigensystem[A1]

So we have three distinct eigenvalues hence the matrix is diagonalizable, i.e., we may find an invertible matrix C such that A=TΔT\(TraditionalForm\`\^\(-1\)\),  where the columns of T are given by the eigenvectors of A and Δ is a diagonal matrix built form the eigenvalues of A.
Let's now calculate C and check that everything has been done correctly. We hereby use the command DiagonalMatrix which does the obvious thing.

T1 := Transpose[Eigenvectors[A1]]
MatrixForm[T1]
Δ1 := 
      DiagonalMatrix[Eigenvalues[A1]]
MatrixForm[Δ1]

MatrixForm[T1 . Δ1 . Inverse[T1]]

yep! 
Let's do another example

A2:=(({{-1, -1, 4}, {-3, 0, 4}, {-3, 3, 1}}))
Eigensystem[A2]

Again we have three distinct complex eigenvalues and the matrix can be diagonalized over the complex numbers.

T2 := Transpose[Eigenvectors[A2]]
MatrixForm[Simplify[Inverse[T2] . A2 . T2]]

Yet another try

A:=(({{1, 0, 1}, {-1, 2, 1}, {-1, 0, 3}}))
Eigensystem[A]

Ups! This matrix cannot be diagonalized! Why?
However, it still can be triangulized over the reals. Why?
To do so we have to complete the system consisting of the two (linearly independent!) eigenvectors of A to a Basis of R\(TraditionalForm\`\^3\). Watchout this gets a little clumsy...

<< LinearAlgebra`MatrixManipulation`
T := Transpose[
        AppendColumns[
          Table[Eigenvectors[A][[j]], {j, 1, 2}], {{0, 0, 
              1}}]]
MatrixForm[T]

MatrixForm[Simplify[Inverse[T] . A . T]]

This is, however, not as good as it can bee. We can transform A to Jordan-form (over the reals, why?). We can either read the Jordan-form off from the result of the Eigensystem command or calculate it using the JordanDecomposition command.

JordanDecomposition[A]

The command returns two matrices. The second one is the actual Jordan-form J of the matrix A while the first one is the according matrix of the basis change R, i.e., we have A=RJR\(TraditionalForm\`\^\(-1\)\). Let's check this.

J := JordanDecomposition[A][[2]]
R := JordanDecomposition[A][[1]]
MatrixForm[J]
MatrixForm[R]

MatrixForm[R . J . Inverse[R]]

Of course we can use the JordanDecomposition command also to calculate the diagonal form of a matrix---if the latter is achievable (!) since in this case the Jordan-form of the matrix is diagonal.

JordanDecomposition[A1]
JordanDecomposition[A2]

Now go ahead and try some examples by yourself!