ExampleC.mws

>    with(linalg):

>    #I show a trade secret of how to produce complicated looking matrices having a nice Jordan form

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

A := matrix([[3, 1, 0, 0, 0, 0, 0], [0, 3, 1, 0, 0, 0, 0], [0, 0, 3, 0, 0, 0, 0], [0, 0, 0, 3, 1, 0, 0], [0, 0, 0, 0, 3, 0, 0], [0, 0, 0, 0, 0, 2, 1], [0, 0, 0, 0, 0, 0, 2]])

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

M := matrix([[1, 0, 3, 0, 3, 0, 3], [0, 1, 1, 1, 1, 1, 1], [0, 0, 1, 0, 3, 0, 3], [0, 0, 0, 1, 1, 1, 1], [0, 0, 0, 0, 1, 1, 1], [0, 0, 0, 0, 0, 1, 3], [0, 0, 0, 0, 0, 0, 1]])

>    #The fact that M has determinant 1 guarantees it's invertible over Z and so that the overall scrambling we are doing is not too bad numerically

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

S := matrix([[0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1], [1, 0, 0, 0, 0, 0, 0]])

>    #S is just a permuation matrix designed to kill the upper diagonal form of the rest of the matrices. And here comes the matrix we'll deal with:

>    B:=multiply(M, inverse(S), A, S, inverse(M));

B := matrix([[5, 0, -6, 3, 9, -9, 24], [1, 3, -2, 1, 2, -2, 5], [3, 0, -6, 1, 17, -15, 42], [1, 0, -3, 3, 6, -5, 14], [1, 0, -3, 0, 9, -5, 14], [3, 0, -9, 0, 18, -15, 51], [1, 0, -3, 0, 6, -6, 20]])

>    factor(charpoly(B, t));factor(minpoly(B, t));

(t-2)^2*(t-3)^5

(t-2)^2*(t-3)^3

>    #From this we know we deal with two spaces. Note that we know that there is one Jordan block corresponding to the eigenvalue 2 and it has size 2, however we don't know yet the basis in which we get that block. The situation for 3 is worse. We only know there's one Jordan block of size 3 and the total size is 5. We could have a priori either a (3, 2) situation or a (3, 1, 1) situation. We first deal with the space associated with the factor (t-3)^5.

>    L1:=kernel(B-3);nops(L1);

L1 := {vector([0, 1, 0, 0, 0, 0, 0]), vector([3, 0, 3, 1, 1, 0, 0])}

2

>    L2:=kernel((B-3)^2);nops(L2);

L2 := {vector([3, 0, 1, 0, 0, 0, 0]), vector([6, 0, 0, 0, 0, 1, 0]), vector([0, 1, 0, 0, 0, 0, 0]), vector([-6, 0, 0, 1, 1, 0, 0])}

4

>    L3:=kernel((B-3)^3);nops(L3);

L3 := {vector([0, 1, 0, 0, 0, 0, 0]), vector([0, 0, 0, 1, 0, 0, 0]), vector([-6, 0, 0, 0, 1, 0, 0]), vector([6, 0, 0, 0, 0, 1, 0]), vector([3, 0, 1, 0, 0, 0, 0])}

5

>    #We have calculated the dimensions of the kernels and bases. In our notation from class we should have C^3 one dimensional (=5-4), C^2 is two dimensionl(=4-2) and contains (B-3)C^3 and C^1 is 2 dimensional (=dim kernel B-3) and contains (B-3)C^2. We span C^3 using a vector in kernel (B-3)^3 and not in kernel (B_3)^2. We try:

>    multiply((B-3)^2, L3[1]);

vector([0, 0, 0, 0, 0, 0, 0])

>    multiply((B-3)^2, L3[2]);

vector([0, 0, 0, 0, 0, 0, 0])

>    multiply((B-3)^2, L3[3]);

vector([0, 1, 0, 0, 0, 0, 0])

>    multiply((B-3)^2, L3[4]);

vector([0, -1, 0, 0, 0, 0, 0])

>    multiply((B-3)^2, L3[5]);

vector([0, 0, 0, 0, 0, 0, 0])

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

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

>    #(Hopefully, when you run Maple again on this example it will indeed choose to include this vector in the basis L3....)

>    #Then (B-3)v3 is a vector in kernel (B-3)^2. It is

>    Uv3:=multiply(B-3, v3);

Uv3 := vector([-3, -4, -1, 0, 0, 0, 0])

>    #We are looking for one more vector v2 so that kernel (B-3)^2 is the direct sum of kernel(B-3) and Span((B-3)v3, v2). Namely, for a vector in kernel (B-3)^2 and not in the span of the union of L1 and (B-3)v3. We can find such a vector among L2. We are looking for a vector in L2 so that the matrix below has rank 4.

>    matrix([L1[1], L1[2], Uv3, L2[1]]);rank(%);

matrix([[0, 1, 0, 0, 0, 0, 0], [3, 0, 3, 1, 1, 0, 0], [-3, -4, -1, 0, 0, 0, 0], [3, 0, 1, 0, 0, 0, 0]])

3

>    matrix([L1[1], L1[2], Uv3, L2[2]]);rank(%);

matrix([[0, 1, 0, 0, 0, 0, 0], [3, 0, 3, 1, 1, 0, 0], [-3, -4, -1, 0, 0, 0, 0], [6, 0, 0, 0, 0, 1, 0]])

4

>    matrix([L1[1], L1[2], Uv3, L2[3]]);rank(%);

matrix([[0, 1, 0, 0, 0, 0, 0], [3, 0, 3, 1, 1, 0, 0], [-3, -4, -1, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0]])

3

>    matrix([L1[1], L1[2], Uv3, L2[4]]);rank(%);

matrix([[0, 1, 0, 0, 0, 0, 0], [3, 0, 3, 1, 1, 0, 0], [-3, -4, -1, 0, 0, 0, 0], [-6, 0, 0, 1, 1, 0, 0]])

3

>    #So we may take (again, provided Maple chooses the list L2 in the same way when you try it out...)

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

v2 := vector([6, 0, 0, 0, 0, 1, 0])

>    #At this point we are done with our basis for the t-3 factor. It is  

>    #  U2v3= (B-3)^2v3,   Uv3=(B-3)v3,    v3,    Uv2=(B-3)v2,   v2,   

>    #where

>    U2v3:=multiply((B-3)^2, v3);

U2v3 := vector([0, -1, 0, 0, 0, 0, 0])

>    Uv2:=multiply((B-3), v2);

Uv2 := vector([3, 4, 3, 1, 1, 0, 0])

>    #The t-2 factor is easier

>    H1:=kernel((B-2)^2);

H1 := {vector([0, 1, 3, 1, 1, 3, 1]), vector([1, 0, 0, 0, 0, 0, 0])}

>    H2:=kernel(B-2);

H2 := {vector([3, 1, 3, 1, 1, 3, 1])}

>    #We look for a vector in H1 and not in H2. For example, we can take

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

w2 := vector([1, 0, 0, 0, 0, 0, 0])

>    Uw2:=multiply(B-2, w2);

Uw2 := vector([3, 1, 3, 1, 1, 3, 1])

>    #Our basis is then the columns of the matrix

>    M:=transpose(matrix([U2v3, Uv3, v3, Uv2, v2, Uw2, w2]));

M := matrix([[0, -3, -6, 3, 6, 3, 1], [-1, -4, 0, 4, 0, 1, 0], [0, -1, 0, 3, 0, 3, 0], [0, 0, 0, 1, 0, 1, 0], [0, 0, 1, 1, 0, 1, 0], [0, 0, 0, 0, 1, 3, 0], [0, 0, 0, 0, 0, 1, 0]])

>    multiply(inverse(M), B, M);

matrix([[3, 1, 0, 0, 0, 0, 0], [0, 3, 1, 0, 0, 0, 0], [0, 0, 3, 0, 0, 0, 0], [0, 0, 0, 3, 1, 0, 0], [0, 0, 0, 0, 3, 0, 0], [0, 0, 0, 0, 0, 2, 1], [0, 0, 0, 0, 0, 0, 2]])

>    #DONE!!