From ab8cc85adde879fb963c94d15675783f2cf4b183 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Tue, 14 Aug 2007 05:14:52 +0000 Subject: Initial population. --- src/algebra/eigen.spad.pamphlet | 340 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 340 insertions(+) create mode 100644 src/algebra/eigen.spad.pamphlet (limited to 'src/algebra/eigen.spad.pamphlet') diff --git a/src/algebra/eigen.spad.pamphlet b/src/algebra/eigen.spad.pamphlet new file mode 100644 index 00000000..193b58f9 --- /dev/null +++ b/src/algebra/eigen.spad.pamphlet @@ -0,0 +1,340 @@ +\documentclass{article} +\usepackage{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 p^=0 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} -- cgit v1.2.3