Creating worksheets with nice numbers and nice diagonalization Mike May, S.J. Saint Louis University For teaching purposes we want to create matrices that are big enough to be make exercises nontrivial but nice enough to do and with nice enough numbers so that the students are not distracted by the arithmetic. For example, we may want 3 by 3 matrices that diagonalize nicely. (Randomly chosen 3 by 3 matrices often have problems with the step that solves the cubic characteristic polynomial. The answers may quite ugly.) The trick is to work backwards. Start with the matrix in nice form and then find a conjugation with small enough numbers. with(LinearAlgebra):

We start with the diagonal matrix that gives the eigenvalues. MatSize := 3: MaxDiag := 6: f1 := (i,j) -> rand(-MaxDiag..MaxDiag)(): D1 := Matrix(MatSize,MatSize,f1,shape=diagonal); We next create a nice matrix to conjugate D with. For the time being, I think a nice matrix has small integer valued entries and has determinant plus or minus 1 so that its inverse is also nice. We use the method of looking at random matrices until we find a nice one. test := 0: count := 0: MaxValP := 3: while test = 0 do: P1 := RandomMatrix(MatSize,MatSize,generator=rand(-MaxValP..MaxValP)): d1 := abs(Determinant(P1)): if d1 = 1 then test := 1 end if: count := count + 1: end do: count; P1; Q1 := P1^(-1); Now we conjugate to get our matrix. Just to verify, we compute eigenvalues and eigenvectors. The eigenvalues should be the entries of the diagonal matrix. The eigenvectors are scalar multiples of the columns of P1. M1 := P1.D1.Q1; Eigenvectors(M1); It is useful to put all the code in a single block that can be repeated until the answer is sufficiently pretty. MatSize := 3: MaxDiag := 6: MaxValP := 3: f1 := (i,j) -> rand(-MaxDiag..MaxDiag)(): D1 := Matrix(MatSize,MatSize,f1,shape=diagonal): test := 0: while test = 0 do: P1 := RandomMatrix(MatSize,MatSize,generator=rand(-MaxValP..MaxValP)): d1 := abs(Determinant(P1)): if d1 = 1 then test := 1 end if: end do: Q1 := P1^(-1): M1 := P1.D1.Q1: print("D1"=D1, "P1" = P1, "Q1"=Q1, "M1"=M1); Eigenvectors(M1); The diagonal matrix is set as M1. Repeat until you get a matrix that looks nice.

We can vary the technique above if we are interested in matrices that are not diagonalizable, either because they have complex eigenvalues or because of non-trivial Jordan blocks. In those cases we can use the CompanionMatrix command to get the nice form that we start with. Poly1 := (x-1)^2*(x-2); Poly2 := expand(Poly1); deg := degree(Poly2,x); N1 := CompanionMatrix(Poly2,x); Eigenvectors(N1); Once again we put this in a block of code to give us a nice but disguised example. Poly1 := (x-1)^2*(x-2): MaxValP := 6: Poly2 := expand(Poly1): MatSize := degree(Poly2,x): N1 := CompanionMatrix(Poly2,x): test := 0: while test = 0 do: P1 := RandomMatrix(MatSize,MatSize,generator=rand(-MaxValP..MaxValP)): d1 := abs(Determinant(P1)): if d1 = 1 then test := 1 end if: end do: Q1 := P1^(-1): M1 := P1.N1.Q1: print("N1"=N1, "P1" = P1, "Q1"=Q1, "M1"=M1); Eigenvectors(M1); CharacteristicPolynomial(M1,x); Executing the block above with Poly1 = (x-1)^2*(x-2) gives a matrix M1 where Rank(M1-I)=2 and Rank((M1-I)^2)=1. Executing the block above with Poly1 = (x^2-x+1)*(x-2) gives a matrix M1 with complex eigenvectors

There are a number of commands for Matrix construction that are quite useful for special cases. The commands are gathered here. They are simple enough to understand that an example will suffice. I3 := IdentityMatrix(3); Z3 := ZeroMatrix(3); D1 := DiagonalMatrix([2,3,4]); ScalarMatrix(x,2,3); ScalarMatrix(x,2); M2 := JordanBlockMatrix([[1,2],[3,1]]); RandomMatrix(3,4); RandomMatrix(3,4,generator=rand(10)); RandomMatrix(3,4,generator=rand(-10..10)); The general matrix constructor is quite powerful in that you can specify a function of i and j to fill in the entries. f1 := (i,j)-> rand(-5..5)(): Matrix(3,3,f1); With a modification of f1 we can produce triangular and unitriangular matrices f1 := (i,j)-> if i>j then 0 else rand(-5..5)()end if: Matrix(3,3,f1); f1 := (i,j)-> if i>j then 0 elif i=j then 1 else rand(-5..5)()end if: Matrix(3,3,f1);