##---------------------------------------------------------------------- ## Title: PencilReduction ## Created: March 2010 ## Authors: Luu Ba Thang ## Advisor: Laurent Buse ## ---------------------------------------------------------------------- ## Description: Some functions dealing with linearization of a polynomial matrix, ## extracting the regular part of a non square pencil of matrices and computing eigenvalues. ##------------------------------------------------------------------------- ## Version last modified 04/2010 ##------------------------------------------------------------------------- __ := NULL: addtohelp := proc() global __; __ := __, args: end: with(LinearAlgebra): with(MatrixPolynomialAlgebra): with(PolynomialTools): with(Groebner): ######## Extracting a regular part of a pencil matrix A-tB ########## #---------------------------------------------------------------------- addtohelp(FirstCompanionMatrix): addtohelp(SecondCompanionMatrix): #---------------------------------------------------------------------- LinearizationMatrix := proc (M, var) local r, s, A, B, U, V, i; s := Degree(M, var); r := RowDimension(M); A := Transpose(Coeff(M, var, 0)); i := 1; while i < s do B := Transpose(Coeff(M, var, i)); A := `<|>`(A, B); i := i+1 end do; RETURN(A); end proc: ## HELP FirstCompanionMatrix; SecondCompanionMatrix ## CALLING SEQUENCE: ## FirstCompanionMatrix(M,var); SecondCompanionMatrix(M,var) ## DESCRIPTION: ## - Given a univariate polynomial matrix M. ## - var - : variable ## RETURN: A, B: Companion matrices of M. ## EXAMPLES: M := Matrix(2, 2, {(1, 1) = t^2, (1, 2) = t+1, (2, 1) = 1, (2, 2) = t}) ## FirstCompanionMatrix(M, var) ## Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, ## (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, ## (3, 1) = 0, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, ## (4, 1) = 1, (4, 2) = 0, (4, 3) = 1, (4, 4) = 1}) ## SecondCompanionMatrix(M, var) ## Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, ## (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, ## (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, ## (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}) FirstCompanionMatrix := proc (M, var) local m, d1, d, U, V, W, A; m := RowDimension(M); d := Degree(M, var); d1 := (d-1)*m; U := Matrix(d1, m); V := IdentityMatrix(d1); W := `<|>`(U, V); A := `<,>`(W, LinearizationMatrix(M, var)); RETURN(A); end proc: SecondCompanionMatrix := proc (M, var) local m, n, d, d1, U, V, W, T, B1, B2, B; m := RowDimension(M); n := ColumnDimension(M); d := Degree(M, var); d1 := (d-1)*m; U := Matrix(n, d1); V := IdentityMatrix(d1); W := Matrix(d1, m); T := -Transpose(Coeff(M, var, d)); B1 := `<|>`(V, W); B2 := `<|>`(U, T); B := `<,>`(B1, B2); RETURN(B); end proc: ReductionStep := proc (A, B) local m, n, p, r, s, B1, B2, B3, B4, B5, A1, A2, A3, A4, A5, C, C1, C2, P, L, U, P1, L1, U1, E1, E2, R, R1, E, F, M, N; m := RowDimension(B); n := ColumnDimension(B); r := Rank(B); M := A; N := B; while r < n do B1 := Transpose(N); P, L, U, R := LUDecomposition(B1, method = 'RREF'); C := Multiply(Multiply(P, L), U); C1 := Transpose(C); C2 := MatrixInverse(C1,method = LU); B2 := Multiply(N, C2); A1 := Multiply(M, C2); A2 := DeleteColumn(A1, 1 .. r); s := Rank(A2); P1, L1, U1, R1 := LUDecomposition(A2, method = 'RREF'); E1 := Multiply(Multiply(P1, L1), U1); E2 := MatrixInverse(E1,method = LU); A3 := Multiply(E2, A1); B3 := Multiply(Multiply(E2, N), C2); if 0 < s then A4 := DeleteRow(A3, 1 .. s); B4 := DeleteRow(B3, 1 .. s) else A4 := A3; B4 := B3 end if; A5 := DeleteColumn(A4, r+1 .. n); B5 := DeleteColumn(B4, r+1 .. n); M := A5; N := B5; r := Rank(B5); n := ColumnDimension(B5); m := RowDimension(B5) end do; RETURN(A5, B5); end proc: #---------------------------------------------------------------------- addtohelp(RegularPencil): #---------------------------------------------------------------------- ## HELP RegularPencil ## CALLING SEQUENCE: ## RegularPencil(A,B) ## DESCRIPTION: ## - Given a pencil matrix A-t*B ## - A, B: Matrices of size (m,n) with m>=n. ## RETURN: A pair matrix (A',B') such that A'-t*B' is regular pencil matrix. ## EXAMPLES: M := Matrix(2, 2, {(1, 1) = t^2, (1, 2) = t+1, (2, 1) = 1, (2, 2) = t}) ## A:= FirstCompanionMatrix(M, var) ## B:= SecondCompanionMatrix(M, var) ## RegularPencil(A,B) ## Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -1, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 1, (3, 2) = 0, (3, 3) = -1}), ## Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (3, 1) = 0, (3, 2) = -1, (3, 3) = 0}) RegularPencil := proc (A, B) local m, n, p, r, s, i, j, B1, B2, B3, B4, B5, B6, A1, A2, A3, A4, A5, A6, C, C1, C2, P, L, U, P1, L1, U1, E1, E2, R, R1, E, F, M, N; m := RowDimension(B); n := ColumnDimension(B); r := Rank(B); M := A; N := B; while r < n do B1 := Transpose(N); P, L, U, R := LUDecomposition(B1, method = 'RREF'); C := Multiply(Multiply(P, L), U); C1 := Transpose(C); C2 := MatrixInverse(C1,method = LU); B2 := Multiply(N, C2); A1 := Multiply(M, C2); A2 := DeleteColumn(A1, 1 .. r); s := Rank(A2); P1, L1, U1, R1 := LUDecomposition(A2, method = 'RREF'); E1 := Multiply(Multiply(P1, L1), U1); E2 := MatrixInverse(E1, method = LU); A3 := Multiply(E2, A1); B3 := Multiply(Multiply(E2, N), C2); if 0 < s then A4 := DeleteRow(A3, 1 .. s); B4 := DeleteRow(B3, 1 .. s) else A4 := A3; B4 := B3 end if; A5 := DeleteColumn(A4, r+1 .. n); B5 := DeleteColumn(B4, r+1 .. n); M := A5; N := B5; r := Rank(B5); n := ColumnDimension(B5); m := RowDimension(B5) end do; if r = m then RETURN(M, N) else A6 := Transpose(M); B6 := Transpose(N); M, N := ReductionStep(A6, B6) end if; RETURN(M, N); end proc: #---------------------------------------------------------------------- addtohelp(ComputationEigenvalues): #---------------------------------------------------------------------- ## HELP ComputationEigenvalues ## CALLING SEQUENCE: ## ComputationEigenvalues(M,var) ## DESCRIPTION: ## - Given a univariate polynomial matrix M(t) size (m,n) with m<=n. ## - var -: Variable t ## RETURN: List of value t such that rank M(t) drop. ## EXAMPLES:M := Matrix(4, 5, {(1, 1) = t^3, (1, 2) = 0, (1, 3) = t^2-1, (1, 4) = 2*t, (1, 5) = 0, ## (2, 1) = 0, (2, 2) = t^2+1, (2, 3) = 2*t^2, (2, 4) = 0, (2, 5) = 0, ## (3, 1) = 0, (3, 2) = t+1, (3, 3) = 0, (3, 4) = 1, (3, 5) = 1, ## (4, 1) = 0, (4, 2) = 0, (4, 3) = t^2-2, (4, 4) = 0, (4, 5) = 0}) ## ComputationEigenvalues(M,t) ## [0, sqrt(2), -sqrt(2), I, -I] ComputationEigenvalues := proc (M, var) local A, B, C, C1, C2; A := FirstCompanionMatrix(M, var); B := SecondCompanionMatrix(M, var); C1, C2 := RegularPencil(A, B); C := Eigenvalues(C1, C2, output = 'list'); RETURN(C); end proc: [__];