\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra radeigen.spad} \author{Patrizia Gianni} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package REP RadicalEigenPackage} <>= )abbrev package REP RadicalEigenPackage ++ Author: P.Gianni ++ Date Created: Summer 1987 ++ Date Last Updated: October 1992 ++ Basic Functions: ++ Related Constructors: EigenPackage, RadicalSolve ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ Package for the computation of eigenvalues and eigenvectors. ++ This package works for matrices with coefficients which are ++ rational functions over the integers. ++ (see \spadtype{Fraction Polynomial Integer}). ++ The eigenvalues and eigenvectors are expressed in terms of radicals. RadicalEigenPackage() : C == T where R ==> Integer P ==> Polynomial R F ==> Fraction P RE ==> Expression R SE ==> Symbol() M ==> Matrix(F) MRE ==> Matrix(RE) ST ==> SuchThat(SE,P) NNI ==> NonNegativeInteger EigenForm ==> Record(eigval:Union(F,ST),eigmult:NNI,eigvec:List(M)) RadicalForm ==> Record(radval:RE,radmult:Integer,radvect:List(MRE)) C == with radicalEigenvectors : M -> List(RadicalForm) ++ radicalEigenvectors(m) computes ++ the eigenvalues and the corresponding eigenvectors of the ++ matrix m; ++ when possible, values are expressed in terms of radicals. radicalEigenvector : (RE,M) -> List(MRE) ++ radicalEigenvector(c,m) computes the eigenvector(s) of the ++ matrix m corresponding to the eigenvalue c; ++ when possible, values are ++ expressed in terms of radicals. radicalEigenvalues : M -> List RE ++ radicalEigenvalues(m) computes the eigenvalues of the matrix m; ++ when possible, the eigenvalues are expressed in terms of radicals. eigenMatrix : M -> Union(MRE,"failed") ++ eigenMatrix(m) returns the matrix b ++ such that \spad{b*m*(inverse b)} is diagonal, ++ or "failed" if no such b exists. normalise : MRE -> MRE ++ normalise(v) returns the column ++ vector v ++ divided by its euclidean norm; ++ when possible, the vector v is expressed in terms of radicals. gramschmidt : List(MRE) -> List(MRE) ++ gramschmidt(lv) converts the list of column vectors lv into ++ a set of orthogonal column vectors ++ of euclidean length 1 using the Gram-Schmidt algorithm. orthonormalBasis : M -> List(MRE) ++ orthonormalBasis(m) returns the orthogonal matrix b such that ++ \spad{b*m*(inverse b)} is diagonal. ++ Error: if m is not a symmetric matrix. T == add PI ==> PositiveInteger RSP := RadicalSolvePackage R import EigenPackage R ---- Local Functions ---- evalvect : (M,RE,SE) -> MRE innerprod : (MRE,MRE) -> RE ---- eval a vector of F in a radical expression ---- evalvect(vect:M,alg:RE,x:SE) : MRE == n:=nrows vect xx:=kernel(x)$Kernel(RE) w:MRE:=zero(n,1)$MRE for i in 1..n repeat v:=eval(vect(i,1) :: RE,xx,alg) setelt(w,i,1,v) w ---- inner product ---- innerprod(v1:MRE,v2:MRE): RE == (((transpose v1)* v2)::MRE)(1,1) ---- normalization of a vector ---- normalise(v:MRE) : MRE == normv:RE := sqrt(innerprod(v,v)) normv = 0$RE => v (1/normv)*v ---- Eigenvalues of the matrix A ---- radicalEigenvalues(A:M): List(RE) == x:SE :=new()$SE pol:= characteristicPolynomial(A,x) :: F radicalRoots(pol,x)$RSP ---- Eigenvectors belonging to a given eigenvalue ---- ---- expressed in terms of radicals ---- radicalEigenvector(alpha:RE,A:M) : List(MRE) == n:=nrows A B:MRE := zero(n,n)$MRE for i in 1..n repeat for j in 1..n repeat B(i,j):=(A(i,j))::RE B(i,i):= B(i,i) - alpha [v::MRE for v in nullSpace B] ---- eigenvectors and eigenvalues ---- radicalEigenvectors(A:M) : List(RadicalForm) == leig:List EigenForm := eigenvectors A n:=nrows A sln:List RadicalForm := empty() veclist: List MRE for eig in leig repeat eig.eigval case F => veclist := empty() for ll in eig.eigvec repeat m:MRE:=zero(n,1) for i in 1..n repeat m(i,1):=(ll(i,1))::RE veclist:=cons(m,veclist) sln:=cons([(eig.eigval)::F::RE,eig.eigmult,veclist]$RadicalForm,sln) sym := eig.eigval :: ST xx:= lhs sym lval : List RE := radicalRoots((rhs sym) :: F ,xx)$RSP for alg in lval repeat nsl:=[alg,eig.eigmult, [evalvect(ep,alg,xx) for ep in eig.eigvec]]$RadicalForm sln:=cons(nsl,sln) sln ---- orthonormalization of a list of vectors ---- ---- Grahm - Schmidt process ---- gramschmidt(lvect:List(MRE)) : List(MRE) == lvect=[] => [] v:=lvect.first n := nrows v RMR:=RectangularMatrix(n:PI,1,RE) orth:List(MRE):=[(normalise v)] for v: local in lvect.rest repeat pol:=((v:RMR)-(+/[(innerprod(w,v)*w):RMR for w in orth])):MRE orth:=cons(normalise pol,orth) orth ---- The matrix of eigenvectors ---- eigenMatrix(A:M) : Union(MRE,"failed") == lef:List(MRE):=[:eiv.radvect for eiv in radicalEigenvectors(A)] n:=nrows A #lef "failed" d:MRE:=copy(lef.first) for v in lef.rest repeat d:=(horizConcat(d,v))::MRE d ---- orthogonal basis for a symmetric matrix ---- orthonormalBasis(A:M):List(MRE) == not symmetric?(A) => error "the matrix is not symmetric" basis:List(MRE):=[] lvec:List(MRE) := [] alglist:List(RadicalForm):=radicalEigenvectors(A) n:=nrows A for alterm in alglist repeat if (lvec:=alterm.radvect)=[] then error "sorry " if #(lvec)>1 then lvec:= gramschmidt(lvec) basis:=[:lvec,:basis] else basis:=[normalise(lvec.first),:basis] basis @ \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}