(*:Version: Mathematica 2.0 *) (*:Name: KnoxPackages`LinearAlgebra` *) (*:Author: Dennis M. Schneider with assistance from Robby S. Villegas *) (*:Summary: Defines some functions useful for linear algebra. *) (* Copyright 1990-92 by Dennis M. Schneider *) (* Work on this package partially supported by grants from Pew Charitable Trusts National Science Foundation NSF-ILI Grant # USE-9050757 NSF Grant # USE-9153249 Lilly Foundation Knox College *) BeginPackage["KnoxPackages`LinearAlgebra`",{"KnoxPackages`CommonFunctions`","Utilities`FilterOptions`"}]; PivotColumns::usage = "PivotColumns[M] returns the column numbers of the pivot columns of the matrix M." ; ColumnSpace::usage = "ColumnSpace[M] returns a list of vectors that make up a basis for the column space of M." ; Rank::usage = "Rank[M] returns the number of nonzero rows in the reduced echelon form of M." ; RowSpace::usage = "RowSpace[M] returns the list of the nonzero rows of the reduced echelon form of M." ; AppendColumn::usage = "AppendColumn[matrix,vector] appends the vector as the last column of the matrix. AppendColumn[matrix,{vector1,vector2,...}] appends vector1, vector2, ... as the last columns of the matrix."; AppendRow::usage = "AppendRow[matrix,vector] appends the vector as the last row of the matrix. AppendRow[matrix,{vector1,vector2,...}] appends vector1, vector2, ... as rows of the matrix."; GramSchmidt::usage = "GramSchmidt[{v1, v2, ...}] returns the the orthogonal set of vectors that results from applying the standard Gram Schmidt algorithm to the list of vectors. GramSchmidt[{v1, v2, ...},Normalized -> True] returns the orthonormal set of vectors that results from applying GramSchmidt. The inner product can be set using the option InnerProduct. Possible options are: NDot, FunctionDot, NFunctionDot, WeightedDot, NWeightedDot, ComplexDot, NComplexDot."; NDot::usage = "InnerProduct -> NDot[v1,v2] uses the inner product N[v1].N[v2]."; FunctionDot::usage = "InnerProduct -> FunctionDot[{x,a,b},weight] uses the inner product Integrate[fun1 fun2 weight,{x,a,b}]."; NFunctionDot::usage = "InnerProduct -> NFunctionDot[{x,a,b},weight] uses the inner product NIntegrate[fun1 fun2 weight,{x,a,b}]."; WeightedDot::usage = "InnerProduct -> WeightedDot[weight] uses the inner product v1.v2.weight."; NWeightedDot::usage = "InnerProduct -> NWeightedDot[weight] uses the inner product N[v1].N[v2].N[weight]."; ComplexDot::usage = "InnerProduct -> ComplexDot uses the inner product v1.Conjugate[v2]."; NComplexDot::usage = "InnerProduct -> NComplexDot uses the inner product N[v1].N[Conjugate[v2]]."; Project::usage = "Project[b,{v1, v2, ...}] returns the projection of b onto the subspace spanned by v1, v2, ... ."; ProjectionMatrix::usage = "ProjectionMatrix[{v1,v2,...}] returns the projection matrix for the projection of R^m onto the subspace spanned by v1, v2, ... ."; LeastSquaresSoln::usage = "LeastSquaresSoln[matrix,vector] returns the solution of the normal equations."; DiagQuad::usage = "DiagQuad[quadratic,{varlist},{newvarlist}] returns a list consisting of two elements: the diagonalized quadratic in the variables specified in newvarlist and the eigenvectors of the quadratic. If newvarlist is not specified, variables not previously used in the current Mathematica session will be chosen."; NDiagQuad::usage = "NDiagQuad[quadratic,{varlist},{newvarlist}] is similar to DiagQuad except that numerical approximations are used in all calculations."; CompleteSquare::usage = "CompleteSquare[expr, var] completes the square in the variable specified by var. CompleteSquare[expr,varlist] completes the square in all the variables specified in varlist."; HilbertMatrix::usage = "HilbertMatrix[m,n] gives the m x n Hilbert matrix. HilbertMatrix[n] gives an nx n Hilbert matrix."; RandomVector::usage = "RandomVector[n] returns an n-vector whose components are Random[]. As in the built-in Random, the type and range of the random numbers can be specified, using the form RandomVector[n, type, range]." ; RandomMatrix::usage = "RandomMatrix[{n1,n2}] generates a matrix of dimensions n1 by n2 with entries Random[]. RandomMatrix[{n1,n2}, type, range, (opts)] fills the entries of the matrix with Random[type, range]. The option MatRank->k returns a random matrix with rank k. RandomMatrix[n] returns a square matrix whose entries are random. RandomMatrix[n,type,range,(opts)] returns a square matrix of size n whose entries are Random[type,range]." ; RandomSymmetricMatrix::usage = "RandomSymmetricMatrix[n] returns an nxn symmetric matrix whose entries are Random[]. RandomSymmetricMatrix[n, type, range] returns an nxn symmetric matrix whose entries are Random[type, range]."; RandomTriangularMatrix::usage = "RandomTriangularMatrix[n] returns an nxn upper triangular matrix whose entries are Random[]. RandomTriangularMatrix[n, type, range] returns an nxn upper triangular matrix whose entries are Random[type, range]."; RandomOrthogonalMatrix::usage = "RandomOrthogonalMatrix[n] returns an nxn orthogonal matrix whose entries are Random[]. RandomOrthogonalMatrix[n, type, range] returns an nxn orthogonal matrix whose entries are Random[type, range]."; RandomNonsingularMatrix::usage = "RandomNonsingularMatrix[n] returns an nxn nonsingular matrix whose entries are Random[]. RandomNonsingularMatrix[n, type, range] returns an nxn nonsingular matrix whose entries are Random[type, range]."; RandomSingularMatrix::usage = "RandomSingularMatrix[n] returns an nxn singular matrix whose entries are Random[]. RandomSingularMatrix[n, type, range] returns an nxn singular matrix whose entries are Random[type, range]."; RandomPermutationMatrix::usage = "RandomPermutationMatrix[n] returns an nxn permutation matrix."; Options[RandomMatrix] = {Rank->{}}; Options[GramSchmidt] = {Normalized -> False, InnerProduct -> Dot, Tolerance -> 10^-10}; Begin["`Private`"]; Poly::notpoly = "The input expression `1` is not a polynomial in the specified variable(s) `2`."; Poly::degree2 = "The input polynomial `1` is not a polynomial of degree 2 in the specified variable(s) `2`."; GramSchmidt::Normalized = "The setting `` for the Normalized option is invalid; this option must be True or False."; GramSchmidt::Tolerance = "The specified tolerance `` is invalid; Tolerance must be or symbolically represent a positive real number."; NDot := Dot[N[#1],N[#2]]& FunctionDot[{x_,a_,b_},wt___] := Integrate[#1 #2 wt,{x,a,b}]& NFunctionDot[{x_,a_,b_},wt___] := NIntegrate[#1 #2 wt,{x,a,b}]& WeightedDot[wt_] := Apply[Plus,Times[#1,#2,wt]]& NWeightedDot[wt_] := Apply[Plus,Times[N[#1],N[#2],N[wt]]]& ComplexDot := Dot[#1,Conjugate[#2]]& NComplexDot := Dot[N[#1],Conjugate[N[#2]]]& nonzerotest[tolerance___] := ( Chop[N[#],tolerance] =!= (0 #) )& firstnonzero[matrix_?MatrixQ] := Map[Position[#,_?(# != 0 &),{1},1,Heads->False]&,matrix]//Flatten PivotColumns[matrix_?MatrixQ ] := firstnonzero[RowReduce[matrix]] ColumnSpace[matrix_?MatrixQ] := Transpose[matrix] [[ PivotColumns[matrix] ]] Rank[matrix_?MatrixQ] := Length[firstnonzero[RowReduce[matrix]] ] RowSpace[matrix_?MatrixQ] := Block[{reduced = RowReduce[matrix]}, reduced[[ Range[Length[firstnonzero[reduced]]] ]] ] AppendColumn[mat:{__List},column:{__List}] := Transpose[Join[Transpose[mat],column]] AppendColumn[mat:{__List},column_List] := AppendColumn[mat,{column}] AppendRow[mat:{__List},row:{__List}] := Join[mat,row] AppendRow[mat:{__List},row_List] := AppendRow[mat,{row}] CompleteSquare[poly_,vars_] := Module[{newpoly,varslist}, varslist = If[Head[vars] === List,vars,{vars}]; If[!PolynomialQ[poly,vars],Return @ Message[Poly::notpoly,poly,vars]]; If[PolynomialDegree[poly,varslist] > 2, Return @ Message[Poly::degree2,poly,vars]]; newpoly = Collect[poly,varslist]; Fold[#1/. a_. #2^2 + b_. #2 -> a (#2 + Simplify[b/2/a])^2 - Simplify[b^2/4/a] &, newpoly,varslist] ] GramSchmidt[spanset:{__},opts___Rule] := Module[{normalized,ip,tolerance,newspanset,orthostep,orthogbasis = {},newvector}, {normalized,ip,tolerance} = {Normalized,InnerProduct,Tolerance} /. {opts} /.Options[GramSchmidt]; If[!MatchQ[normalized, True | False], Return @ Message[GramSchmidt::Normalized, normalized] ]; If[Head[N[tolerance]] =!= Real || N[tolerance] <= 0, Return @ Message[GramSchmidt::Tolerance, tolerance] ]; newspanset = Select[spanset,nonzerotest[tolerance]]; If[newspanset === {},Return[{}]]; If[First @ Characters[ToString[ip]] === "N", newspanset = Chop[N[newspanset],tolerance 10^-2]]; orthostep = Function[{list,y}, newvector = Simplify[y - Apply[Plus,Map[ip[y,#]/ip[#,#] #&,list]]]; If[nonzerotest[tolerance][newvector], orthogbasis = Flatten[{orthogbasis, {newvector}},1],orthogbasis]]; orthogbasis = Chop[Fold[orthostep,{},newspanset],tolerance]; If[normalized, Map[(#/Sqrt[Simplify[ip[#,#]]])&,orthogbasis],orthogbasis] ] Project[x_,spanset:{__},opts___Rule] := Module[{innerprod,orthobasis}, innerprod = InnerProduct /. {opts} /. Options[GramSchmidt]; orthobasis = GramSchmidt[spanset, Normalized -> True, InnerProduct->innerprod,opts]; Apply[Plus,Map[Simplify[(innerprod[x,#] #)]&,orthobasis]] ] ProjectionMatrix[spanset_?MatrixQ] := Module[{matrix = GramSchmidt[spanset,Normalized->True]}, Transpose[matrix].matrix ] LeastSquaresSoln[matrix_?MatrixQ,vector_?VectorQ] := Module[{transposemat = Transpose[matrix]}, LinearSolve[transposemat.matrix, transposemat.vector] ] DiagQuad::varlists = "The original variables `1` and new variables `2` are lists of unequal lengths, `3` and `4` respectively."; DiagQuad::commonvars = "The specified original and new variables have `` in common, when they should be disjoint."; (* 12-11-92 Changed output to {diagquad,evalues,evectors} *) DiagQuad[quadratic_, vars_List, newvars_List] := Block[{quad = quadratic//Expand, length = Length[vars], symmat,evalues,evectors,linearterms,transformlinear,diagquad,newquad}, (* Error-checking of the input: *) If[!PolynomialQ[quadratic, vars], Return @ Message[Poly::notpoly, quadratic, vars] ]; If[PolynomialDegree[quadratic, vars] > 2, Return @ Message[Poly::degree2, quadratic,vars] ]; If[Length[vars] != Length[newvars], Return @ Message[DiagQuad::varlists, vars, newvars, length,Length[newvars]] ]; If[Intersection[vars, newvars] =!= {}, Return @ Message[DiagQuad::commonvars, Intersection[vars, newvars]] ]; (* Begin diagonalization routine: *) symmat = Table[ If[i != j, Coefficient[quad,vars[[i]] vars[[j]]]/2, Coefficient[quad,vars[[i]] vars[[j]]]], {i,length},{j,length}]; {evalues,evectors} = Eigensystem[symmat]//Chop; evectors = GramSchmidt[evectors,Normalized->True]; linearterms = quad - Expand[vars.symmat.vars]; transformlinear = Map[Coefficient[linearterms,#]&,vars]. Transpose[evectors].newvars; diagquad = evalues.(newvars^2)//Cancel; newquad = {diagquad + transformlinear + linearterms/.Map[#->0&,vars],evalues,evectors}//Chop ] DiagQuad[quadratic_, vars_List] := Block[{newvars = Map[Unique[ToString[#]]&,vars]}, DiagQuad[quadratic, vars, newvars]] NDiagQuad[quadratic_, vars_List, newvars_List] := DiagQuad[1. quadratic//Expand//N, vars//N, newvars] NDiagQuad[quadratic_, vars_List] := Block[{newvars = Map[Unique[ToString[#]]&,vars]}, DiagQuad[1. quadratic//Expand//N, vars//N, newvars]] HilbertMatrix[m_Integer?Positive,n_Integer?Positive] := Table[1/(i + j),{i,m},{j,n}] HilbertMatrix[n_Integer?Positive] := HilbertMatrix[n,n] RandomVector[n_Integer?Positive, type_:Real, range_:{0,1}] := Table[Random[type, range], {n}] RandomVector[0, type_:Real, range_:{0,1}] := {} RandomMatrix[{dim1_,dim2_}, type_:Real, range_:{0,1}, opts___] := Block[{matrank,mat}, matrank = Rank /. {opts} /. Options[RandomMatrix]; mat = Table[Random[type, range],{dim1},{dim2} ]; If[matrank === {},mat, While[Rank[mat] != matrank, mat = Table[Random[type, range],{dim1},{dim2} ]]; mat]] RandomMatrix[dim_Integer?Positive, type_:Real, range_:{0,1}, opts___] := RandomMatrix[{dim,dim}, type, range, opts] RandomSymmetricMatrix[n_Integer?Positive, type_:Real, range_:{0,1}] := Module[ {lower, diagonal}, lower = Table[ RandomVector[j - 1, type, range] ~Join~ Table[0, {n - (j - 1)}], {j, 1, n}]; diagonal = DiagonalMatrix[RandomVector[n, type, range]]; lower + Transpose[lower] + diagonal ] RandomTriangularMatrix[n_Integer?Positive, type_:Real, range_:{0,1}] := Table[ Table[0,{j - 1}] ~ Join ~ RandomVector[n - (j - 1), type, range], {j,1,n}] RandomNonsingularMatrix[n_Integer?Positive,type_:Real,range_:{0,1}] := Block[{mat = RandomMatrix[n,type,range]}, While[Rank[mat] < n,mat = RandomMatrix[n,type,range]]; mat ] RandomSingularMatrix[n_Integer?Positive, type_:Real, range_:{0,1}] := Block[{mat = RandomMatrix[n,type,range]}, While[Rank[mat] == n,mat = RandomMatrix[n,type,range]]; mat ] RandomOrthogonalMatrix[n_, type_:Real, range_:{0,1}] := GramSchmidt[RandomNonsingularMatrix[n, type, range], Normalized->True] RandomPermutationMatrix[size_?IntegerQ] := Part[IdentityMatrix[size], Map[Last,Sort[Table[{Random[],i},{i,size}]]]] End[]; Protect[ PivotColumns, ColumnSpace, Rank, RowSpace, AppendColumn, AppendRow, GramSchmidt, Project, ProjectionMatrix, LeastSquaresSoln, DiagQuad, NDiagQuad, CompleteSquare, HilbertMatrix, RandomVector, RandomMatrix, RandomSymmetricMatrix, RandomTriangularMatrix, RandomNonsingularMatrix, RandomOrthogonalMatrix, RandomSingularMatrix, RandomPermutationMatrix, MatRank, Normalized, Dot, NDot, FunctionDot, NFunctionDot, WeightedDot, NWeightedDot, ComplexDot, NComplexDot ]; EndPackage[];