\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra eigen.spad} \author{Patrizia Gianni, Barry Trager} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package EP EigenPackage} <>= )abbrev package EP EigenPackage ++ Author: P. Gianni ++ Date Created: summer 1986 ++ Date Last Updated: October 1992 ++ Basic Functions: ++ Related Constructors: NumericRealEigenPackage, NumericComplexEigenPackage, ++ RadicalEigenPackage ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This is a package for the exact computation of eigenvalues and eigenvectors. ++ This package can be made to work for matrices with coefficients which are ++ rational functions over a ring where we can factor polynomials. ++ Rational eigenvalues are always explicitly computed while the ++ non-rational ones are expressed in terms of their minimal ++ polynomial. -- Functions for the numeric computation of eigenvalues and eigenvectors -- are in numeigen spad. EigenPackage(R) : C == T where R : GcdDomain P ==> Polynomial R F ==> Fraction P SE ==> Symbol() SUP ==> SparseUnivariatePolynomial(P) SUF ==> SparseUnivariatePolynomial(F) M ==> Matrix(F) NNI ==> NonNegativeInteger ST ==> SuchThat(SE,P) Eigenvalue ==> Union(F,ST) EigenForm ==> Record(eigval:Eigenvalue,eigmult:NNI,eigvec : List M) GenEigen ==> Record(eigval:Eigenvalue,geneigvec:List M) C == with characteristicPolynomial : (M,Symbol) -> P ++ characteristicPolynomial(m,var) returns the ++ characteristicPolynomial of the matrix m using ++ the symbol var as the main variable. characteristicPolynomial : M -> P ++ characteristicPolynomial(m) returns the ++ characteristicPolynomial of the matrix m using ++ a new generated symbol symbol as the main variable. eigenvalues : M -> List Eigenvalue ++ eigenvalues(m) returns the ++ eigenvalues of the matrix m which are expressible ++ as rational functions over the rational numbers. eigenvector : (Eigenvalue,M) -> List M ++ eigenvector(eigval,m) returns the ++ eigenvectors belonging to the eigenvalue eigval ++ for the matrix m. generalizedEigenvector : (Eigenvalue,M,NNI,NNI) -> List M ++ generalizedEigenvector(alpha,m,k,g) ++ returns the generalized eigenvectors ++ of the matrix relative to the eigenvalue alpha. ++ The integers k and g are respectively the algebraic and the ++ geometric multiplicity of tye eigenvalue alpha. ++ alpha can be either rational or not. ++ In the seconda case apha is the minimal polynomial of the ++ eigenvalue. generalizedEigenvector : (EigenForm,M) -> List M ++ generalizedEigenvector(eigen,m) ++ returns the generalized eigenvectors ++ of the matrix relative to the eigenvalue eigen, as ++ returned by the function eigenvectors. generalizedEigenvectors : M -> List GenEigen ++ generalizedEigenvectors(m) ++ returns the generalized eigenvectors ++ of the matrix m. eigenvectors : M -> List(EigenForm) ++ eigenvectors(m) returns the eigenvalues and eigenvectors ++ for the matrix m. ++ The rational eigenvalues and the correspondent eigenvectors ++ are explicitely computed, while the non rational ones ++ are given via their minimal polynomial and the corresponding ++ eigenvectors are expressed in terms of a "generic" root of ++ such a polynomial. T == add PI ==> PositiveInteger MF := GeneralizedMultivariateFactorize(SE,IndexedExponents SE,R,R,P) UPCF2:= UnivariatePolynomialCategoryFunctions2(P,SUP,F,SUF) ---- Local Functions ---- tff : (SUF,SE) -> F fft : (SUF,SE) -> F charpol : (M,SE) -> F intRatEig : (F,M,NNI) -> List M intAlgEig : (ST,M,NNI) -> List M genEigForm : (EigenForm,M) -> GenEigen ---- next functions needed for defining ModularField ---- reduction(u:SUF,p:SUF):SUF == u rem p merge(p:SUF,q:SUF):Union(SUF,"failed") == p = q => p p = 0 => q q = 0 => p "failed" exactquo(u:SUF,v:SUF,p:SUF):Union(SUF,"failed") == val:=extendedEuclidean(v,p,u) val case "failed" => "failed" val.coef1 ---- functions for conversions ---- fft(t:SUF,x:SE):F == n:=degree(t) cf:=monomial(1,x,n)$P :: F cf * leadingCoefficient t tff(p:SUF,x:SE) : F == degree p=0 => leadingCoefficient p r:F:=0$F while not zero? p repeat r:=r+fft(p,x) p := reductum p r ---- generalized eigenvectors associated to a given eigenvalue --- genEigForm(eigen : EigenForm,A:M) : GenEigen == alpha:=eigen.eigval k:=eigen.eigmult g:=#(eigen.eigvec) k = g => [alpha,eigen.eigvec] [alpha,generalizedEigenvector(alpha,A,k,g)] ---- characteristic polynomial ---- charpol(A:M,x:SE) : F == dimA :PI := (nrows A):PI dimA ~= ncols A => error " The matrix is not square" B:M:=zero(dimA,dimA) for i in 1..dimA repeat for j in 1..dimA repeat B(i,j):=A(i,j) B(i,i) := B(i,i) - monomial(1$P,x,1)::F determinant B -------- EXPORTED FUNCTIONS -------- ---- characteristic polynomial of a matrix A ---- characteristicPolynomial(A:M):P == x:SE:=new()$SE numer charpol(A,x) ---- characteristic polynomial of a matrix A ---- characteristicPolynomial(A:M,x:SE) : P == numer charpol(A,x) ---- Eigenvalues of the matrix A ---- eigenvalues(A:M): List Eigenvalue == x:=new()$SE pol:= charpol(A,x) lrat:List F :=empty() lsym:List ST :=empty() for eq in solve(pol,x)$SystemSolvePackage(R) repeat alg:=numer lhs eq degree(alg, x)=1 => lrat:=cons(rhs eq,lrat) lsym:=cons([x,alg],lsym) append([lr::Eigenvalue for lr in lrat], [ls::Eigenvalue for ls in lsym]) ---- Eigenvectors belonging to a given eigenvalue ---- ---- the eigenvalue must be exact ---- eigenvector(alpha:Eigenvalue,A:M) : List M == alpha case F => intRatEig(alpha::F,A,1$NNI) intAlgEig(alpha::ST,A,1$NNI) ---- Eigenvectors belonging to a given rational eigenvalue ---- ---- Internal function ----- intRatEig(alpha:F,A:M,m:NNI) : List M == n:=nrows A B:M := zero(n,n)$M for i in 1..n repeat for j in 1..n repeat B(i,j):=A(i,j) B(i,i):= B(i,i) - alpha [v::M for v in nullSpace(B**m)] ---- Eigenvectors belonging to a given algebraic eigenvalue ---- ------ Internal Function ----- intAlgEig(alpha:ST,A:M,m:NNI) : List M == n:=nrows A MM := ModularField(SUF,SUF,reduction,merge,exactquo) AM:=Matrix MM x:SE:=lhs alpha pol:SUF:=unitCanonical map(coerce,univariate(rhs alpha,x))$UPCF2 alg:MM:=reduce(monomial(1,1),pol) B:AM := zero(n,n) for i in 1..n repeat for j in 1..n repeat B(i,j):=reduce(A(i,j)::SUF,pol) B(i,i):= B(i,i) - alg sol: List M :=empty() for vec in nullSpace(B**m) repeat w:M:=zero(n,1) for i in 1..n repeat w(i,1):=tff((vec.i)::SUF,x) sol:=cons(w,sol) sol ---- Generalized Eigenvectors belonging to a given eigenvalue ---- generalizedEigenvector(alpha:Eigenvalue,A:M,k:NNI,g:NNI) : List M == alpha case F => intRatEig(alpha::F,A,(1+k-g)::NNI) intAlgEig(alpha::ST,A,(1+k-g)::NNI) ---- Generalized Eigenvectors belonging to a given eigenvalue ---- generalizedEigenvector(eigen :EigenForm,A:M) : List M == generalizedEigenvector(eigen.eigval,A,eigen.eigmult,# eigen.eigvec) ---- Generalized Eigenvectors ----- generalizedEigenvectors(A:M) : List GenEigen == n:= nrows A leig:=eigenvectors A [genEigForm(leg,A) for leg in leig] ---- eigenvectors and eigenvalues ---- eigenvectors(A:M):List(EigenForm) == n:=nrows A x:=new()$SE p:=numer charpol(A,x) MM := ModularField(SUF,SUF,reduction,merge,exactquo) AM:=Matrix(MM) ratSol : List EigenForm := empty() algSol : List EigenForm := empty() lff:=factors factor p for fact in lff repeat pol:=fact.factor degree(pol,x)=1 => vec:F :=-coefficient(pol,x,0)/coefficient(pol,x,degree(pol,x)) ratSol:=cons([vec,fact.exponent :: NNI, intRatEig(vec,A,1$NNI)]$EigenForm,ratSol) alpha:ST:=[x,pol] algSol:=cons([alpha,fact.exponent :: NNI, intAlgEig(alpha,A,1$NNI)]$EigenForm,algSol) append(ratSol,algSol) @ \section{package CHARPOL CharacteristicPolynomialPackage} <>= )abbrev package CHARPOL CharacteristicPolynomialPackage ++ Author: Barry Trager ++ Date Created: ++ Date Last Updated: ++ Basic Functions: ++ Related Constructors: ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package provides a characteristicPolynomial function ++ for any matrix over a commutative ring. CharacteristicPolynomialPackage(R:CommutativeRing):C == T where PI ==> PositiveInteger M ==> Matrix R C == with characteristicPolynomial: (M, R) -> R ++ characteristicPolynomial(m,r) computes the characteristic ++ polynomial of the matrix m evaluated at the point r. ++ In particular, if r is the polynomial 'x, then it returns ++ the characteristic polynomial expressed as a polynomial in 'x. T == add ---- characteristic polynomial ---- characteristicPolynomial(A:M,v:R) : R == dimA :PI := (nrows A):PI dimA ~= ncols A => error " The matrix is not square" B:M:=zero(dimA,dimA) for i in 1..dimA repeat for j in 1..dimA repeat B(i,j):=A(i,j) B(i,i) := B(i,i) - v determinant B @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}