diff options
| author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
|---|---|---|
| committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
| commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
| tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/radeigen.spad.pamphlet | |
| download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz | |
Initial population.
Diffstat (limited to 'src/algebra/radeigen.spad.pamphlet')
| -rw-r--r-- | src/algebra/radeigen.spad.pamphlet | 235 |
1 files changed, 235 insertions, 0 deletions
diff --git a/src/algebra/radeigen.spad.pamphlet b/src/algebra/radeigen.spad.pamphlet new file mode 100644 index 00000000..04e006fb --- /dev/null +++ b/src/algebra/radeigen.spad.pamphlet @@ -0,0 +1,235 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra radeigen.spad} +\author{Patrizia Gianni} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{package REP RadicalEigenPackage} +<<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 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 <n => "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) == + ^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} +<<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. +@ +<<*>>= +<<license>> + +<<package REP RadicalEigenPackage>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |