STA: "example.Diag.2x2.txt"

(use-ring Rational-ring)

;; Let's create a 2x2 diagonal matrix with eVals 3 and -1.
(setq D (mat-make-diag '(3 -1)))
    [ 3   0 ]
    [ 0  -1 ]

;; We now disguise the matrix.
;; I used   (mat-make-random-detone 2 2)
;; to make t.fol matrix:
(setq InvC (mat-make-square '(-1 -2 1 1)))
    [ -1  -2 ]
    [  1   1 ]

(setq C (mat-inverse? InvC))
    [  1   2 ]
    [ -1  -1 ]


;; Matrix G is similar [conjugate-to] matrix D.
(setq G (mat-mul C D InvC))
    [ -5  -8 ]
    [  4   7 ]

;; We're told that 3 and -1 are the eVals of G, and that our goal
;; is to diagonalize G.
;;   We'll subtract the corresponding scalar matrices from G,
;; then compute bases for the corresponding  nullspaces, which are
;; bases for the corresponding eigenspaces of G.
(setq I (mat-eye 2))
    [ 1  0 ]
    [ 0  1 ]


;;;;;;;;;;;;;;;;
;; We examine matrix   G - 3I.
;;;;;;;;;;;;;;;;
(setq Gm3I (mat-sub G (mat-scal-mult 3 I)))
    [  4   4 ]
    [ -8  -8 ]

;; The following "scrutinize-matrix" RREFs the given matrix.
;; If the matrix is square, then prints
;; info about the matrix, e.g, determinant
;; and bases of various subspaces.
(scrutinize-matrix Gm3I)

  JK: Found 1 pivots before the second column.

Matrix Gm3I,
    [ -8  -8 ]
    [  4   4 ]

has determinant 0 and consequently has no inverse-matrix.

    Dultiplying Gm3I by its RowOp matrix

    [ -1/8  0 ]
    [  1/2  1 ]

produces   R := RREF(Gm3I)   matrix

    [ 1  1 ]
    [ 0  0 ]

with pivot columns  0, and free columns  1.

   A basis for LeftNull(Gm3I) is
    [ -1 ]
    [  1 ]


   A basis for ColSpn(Gm3I) is
    [ -8 ]
    [  4 ]

   A basis for RowSpn(Gm3I) is
    [ 1  1 ]


;;;;;;;;;;;;;;;;
;; We now examine matrix   G - (-1)*I.
;;;;;;;;;;;;;;;;
(setq G[-1]I (mat-sub G (mat-scal-mult -1 I)))
    [ -4  -8 ]
    [  4   8 ]

(scrutinize-matrix G[-1]I)

  JK: Found 1 pivots before the second column.

Matrix G[-1]I,
    [ -4  -8 ]
    [  4   8 ]

has determinant 0 and consequently has no inverse-matrix.

    Multiplying G[-1]I by its RowOp matrix

    [ -1/4  0 ]
    [    1  1 ]

produces   R := RREF(G[-1]I)   matrix

    [ 1  2 ]
    [ 0  0 ]

with pivot columns  0, and free columns  1.

   A basis for LeftNull(G[-1]I) is
    [ -2 ]
    [  1 ]


   A basis for ColSpn(G[-1]I) is
    [ -4 ]
    [  4 ]

   A basis for RowSpn(G[-1]I) is
    [ 1  2 ]



;; Gluing together the bases of
;;   LNul( Gm3I )    and   LNul( G[-1]I )
;; gives us an eigenbasis for our original matrix G.
;; For some reason, I chose to call this CoB matrix "InvC";
;; I don't remember why.

(setq InvC (mat-make-square '(-1 -2 1 1)))
    [ -1  -2 ]
    [  1   1 ]

(setq C (mat-inverse? InvC))
    [  1   2 ]
    [ -1  -1 ]

(setq DiagonalizedG (mat-mul C G InvC))
    [ 3   0 ]
    [ 0  -1 ]

;; Yay!  We succeeded in diagonalizing G, recovering
the diagonal matrix D that we started with.


;;;; ================
;; Attempting to diagonalize a 4x4 matrix.
(use-ring Rational-ring)

(setq G (mat-make-square '(-1 -12 8 4 1 8 -4 -3 -2 -2 3 -2 1 3 -2 2)))
    [ -1  -12   8   4 ]
    [  1    8  -4  -3 ]
    [ -2   -2   3  -2 ]
    [  1    3  -2   2 ]

(setq *POLY-var* "t")
(setq f (mat-CharPoly G))
  t^4 - 12t^3 + 54t^2 - 108t + 81

;; Maybe some small integer is an f-root.
(iter (for v from -4 to 4) (format t "~% f(~2D) = ~3D" v (poly-eval f v)) )

 f(-4) = 2401
 f(-3) = 1296
 f(-2) = 625
 f(-1) = 256
 f( 0) =  81
 f( 1) =  16
 f( 2) =   1
 f( 3) =   0  ;; Yes! Value 3 is an f-root.
 f( 4) =   1

(setq I (mat-eye 4))
    [ 1  0  0  0 ]
    [ 0  1  0  0 ]
    [ 0  0  1  0 ]
    [ 0  0  0  1 ]

(setq G3 (mat-sub G (mat-scal-mult 3 I)))
    [ -4  -12   8   4 ]
    [  1    5  -4  -3 ]
    [ -2   -2   0  -2 ]
    [  1    3  -2  -1 ]

(scrutinize-matrix G3)

  JK: Found 2 pivots before the fourth column.

Matrix G3,

    [ -4  -12   8   4 ]
    [  1    5  -4  -3 ]
    [ -2   -2   0  -2 ]
    [  1    3  -2  -1 ]

has determinant 0 and consequently has no inverse-matrix.

    Multiplying G3 by its RowOp matrix

    [ -5/8  -3/2  0  0 ]
    [  1/8   1/2  0  0 ]
    [   -1    -2  1  0 ]
    [  1/4     0  0  1 ]

produces   R := RREF(G3)   matrix

    [ 1  0   1   2 ]
    [ 0  1  -1  -1 ]
    [ 0  0   0   0 ]
    [ 0  0   0   0 ]

with pivot columns  0, 1, and free columns  2, 3.

   A basis for LeftNull(G3) is
    [ -1  -2 ]
    [  1   1 ]
    [  1   0 ]
    [  0   1 ]

;; I erased the printout of bases for ColSpn and RowSpn

;; Let's store that above nullspace-basis, as a 4x2 matrix.
(setq G3LNull (MTB-LNullBsis *Tab*))
    [ -1  -2 ]
    [  1   1 ]
    [  1   0 ]
    [  0   1 ]

;; We prepare to glue on two more column,
;; to extend the above to a basis for the whole space.
(setq 4by2 (mat-make-optval 4 2 0)) => [ 0  0 ]
                                       [ 0  0 ]
                                       [ 0  0 ]
                                       [ 0  0 ]

(setf (mat-E 4by2 2 0) 1) (setf (mat-E 4by2 3 1) 1)
4by2  =>  [ 0  0 ]
          [ 0  0 ]
          [ 1  0 ]
          [ 0  1 ]

(setq InvD (mat-Horiz-concat G3LNull 4by2))
    [ -1  -2  0  0 ]
    [  1   1  0  0 ]
    [  1   0  1  0 ]
    [  0   1  0  1 ]

(setq D (mat-inverse? InvD))
    [  1   2  0  0 ]
    [ -1  -1  0  0 ]
    [ -1  -2  1  0 ]
    [  1   1  0  1 ]

;; We have /partly/  diagonalized G.
(setq GDiagABit (mat-mul D G InvD))
    [ 3  0   0  -2 ]
    [ 0  3  -4  -1 ]
    [ 0  0   3   0 ]
    [ 0  0   2   3 ]


;;  We don't know if  G is diagonalizable.  Let's
;; explore further.
;; Recall the CharPoly of G, made by (setq f (mat-CharPoly G))
f  =>   t^4 - 12t^3 + 54t^2 - 108t + 81


;; The GeometricMultiplicity is GeoMult_G(3) = 2.
;; Hence AlgMult_G(3) >= 2.  

;;   Let's now divide   [t - 3]^2   into  f(t), 
;; then try to factor the quotient polynomial.
(setq h (poly-build 1 -3))  =>  t-3

(setq DoubleRoot3 (poly-mult-two h h)) =>  t^2 - 6t + 9

;; Before we compute the quotient polynomial q(), let's
;; figure out the possibilities.
;;
;;   If q() has roots b and c such that {3,b,c} are three distinct
;; values, then G is similar to Diag(3,3,b,c)
;;
;;   Alternatively, if q() has a double-root b, with b =/= 3,
;; then G may or may-not be diagonalizable, depending on
;; whether GeoMult_G(b) is 2 or is only 1.
;; 
;;  Lastly, if 3 is either a 
;;            single or double root of q(),
;; then G is NOT diagonalizable, as AlgMult_G(3) is 
;;         either 3  or   4.
;; In either case, AlgMult_G(3) exceeds GeoMult_G(3).


;; Let's now compute the quotient polynomial:
(setq QuoPoly (poly-divide f DoubleRoot3))
  t^2 - 6t + 9         <ZIP poly>
;; Quotient.             Remainder.

;; Whoa!  So QuoPoly equals f, hence 
;;         AlgMult_G(3) = 4.
;; The Jordan Canonical Form theory (which we don't know yet)
;; says that G is conjugate-to either

    [ 3  1  0  0 ]      [ 3  1  0  0 ]
    [ 0  3  0  0 ]      [ 0  3  1  0 ]
    [ 0  0  3  1 ]  or  [ 0  0  3  0 ]
    [ 0  0  0  3 ]      [ 0  0  0  3 ]



;; ================
;; Example 7 of Section 5.1 (P.252).
;; Let \V be the VS of 3-topped polya
;;   Ax^2 + Bx + C.
;;
;; LinTrn T carries poly f to
;;    x |->  f(x) + [x+1]*f'(x) .
;;
;; Our goal is to determine if T is diagonalizable.
;;
;;; Soln:  W.r.t basis  {1, x, x^2}, our T has matrix

M  =>  [ 1  1  0 ]
       [ 0  2  2 ]
       [ 0  0  3 ]

;; This upper-triangular matrix has distinct eVals, so
;; we know that M is conjugate-to
       [ 1       ]
       [    2    ]
       [       3 ]

;; ... but let's pretend that we don't know this.

(setq MCharPly (mat-CharPoly M)   f MCharPly)
  -t^3 + 6t^2 -11t + 6

;; We search for roots of MCharPly:
(iter (for j below 5) (format t "~% f(~D) = ~D"  j (poly-eval f j)))

 f(0) = 6
 f(1) = 0   ;; A root.
 f(2) = 0   ;; Another, and
 f(3) = 0   ;; another.
 f(4) = -6

;; So  f(t) = [1 - t][2 - t][3 - t] .


(setq I (mat-eye 3))
    [ 1  0  0 ]
    [ 0  1  0 ]
    [ 0  0  1 ]

(setq M1 (mat-sub M I))
    [ 0  1  0 ]
    [ 0  1  2 ]
    [ 0  0  2 ]

(setq M2 (mat-sub M (mat-scal-mult 2 I)))
    [ -1  1  0 ]
    [  0  0  2 ]
    [  0  0  1 ]

(setq M3 (mat-sub M (mat-scal-mult 3 I)))
    [ -2   1  0 ]
    [  0  -1  2 ]
    [  0   0  0 ]


(scrutinize-matrix M1)        [ 1 ]
 A basis for LeftNull(M1) is  [ 0 ]
                              [ 0 ]
    
(scrutinize-matrix M2)        [ 1 ]
 A basis for LeftNull(M2) is  [ 1 ]
                              [ 0 ]

(scrutinize-matrix M3)        [ 1 ]
 A basis for LeftNull(M3) is  [ 2 ]
                              [ 1 ]
    
 ;; We define C-inverse to be
InvC  =>  [ 1  1  1 ]
          [ 0  1  2 ]
          [ 0  0  1 ]
                                [ 1  -1   1 ]
(setq C (mat-inverse? InvC)) => [ 0   1  -2 ]
                                [ 0   0   1 ]
    
(mat-mul C M InvC)  =>  [ 1  0  0 ]
                        [ 0  2  0 ]
                        [ 0  0  3 ]

;; ...so we've diagonalized T, the map
;;
;;  f   |->   [ x |->  f(x) + [x+1]*f'(x) ]


END: "example.Diag.2x2.txt"