share/linearalgebra/linalg.usage

   1 Installation
   2
   3   To use linalg, you'll need to have nset version 1.203.  If you
   4   are using Maxima 5.9.1cvs or better, you already have the required
   5   version of nset. Otherwise you'll need to download and install nset.
   6   You may download nset from http://www.unk.edu/facstaff/profiles/willisb/.
   7
   8   You will need to append the path to linearalgebra to
   9   file_search_maxima and file_search_lisp. If the path to
  10   the directory linearalgebra is "c:/maxima/linearalgebra", you will
  11   use the Maxima commands
  12
  13     (%i1) file_search_maxima : cons("c:/maxima/linearalgebra/###.{mac}",file_search_maxima)$
  14     (%i2) file_search_lisp : cons("c:/maxima/linearalgebra/###.{lisp}",file_search_lisp)$
  15
  16   Now you should be able to load and use linalg.
  17
  18   (%i3) load("linalg");
  19         Warning - you are redefining the MACSYMA function EIGENVALUES
  20         Warning - you are redefining the MACSYMA function EIGENVECTORS
  21         Warning - you are redefining the MACSYMA function RANK
  22   (%o3) c:/maxima/linearalgebra/linalg.mac
  23   (%i4) m : matrix([1-z,2],[5,8-z]);
  24   (%o4) matrix([1-z,2],[5,8-z])
  25   (%i5) ptriangularize(m,z);
  26   (%o5) matrix([5,8-z],[0,-z^2/5+(9*z)/5+2/5])
  27
  28
  29 columnop(m,i,j,theta)
  30
  31   If m is a matrix, return the matrix that results from doing the
  32   column operation Ci <- Ci - theta * Cj. If m doesn't have a row
  33   i or j, signal an error.
  34
  35 columnswap(m,i,j)
  36
  37   If m is a matrix, swap columns i and j.  If m doesn't have a column
  38   i or j, signal an error.
  39
  40 columnspace(m)
  41
  42    If m is a matrix, return span(v1,v2,...,vn), where the set
  43    {v1,v2,...,vn} is a basis for the column space of m.  The span
  44    of the empty set is {0}. Thus, when the column space has only
  45    one member, return span().
  46
  47 ctranspose(m)
  48
  49    Return the conjugate of the transpose of the matrix m. If m is
  50    a block matrix, transpose and conjugate each block.
  51
  52 dotproduct(u,v)
  53
  54    Return the dotproduct of vectors u and v.  This is the same
  55    as conj(transpose(u)) . v.  The arguments u and v must be
  56    n x 1 matrices.
  57
  58 get_lu_factors(x)
  59
  60   When x = lu_factor(A), then get_lu_factors returns a list of the
  61   form [P, L, U], where P is a permutation matrix, L is lower triangular with
  62   ones on the digonal, and U is upper triangular. And A = P L U.
  63
  64 hilbert_matrix(n)
  65
  66   Return the n by n Hilbert matrix. When n isn't a nonnegative
  67   integer, signal an error.
  68
  69 kronecker_product(a,b)
  70
  71    Return the Kronecker product of the matrices a and b.
  72
  73 listp (e, {pred})
  74
  75   Given an optional argument pred, return true if e is
  76   a Maxima list and pred evaluates to true for every list element.
  77   When listp is not given the optional argument, return true if e is
  78   a Maxima list.  In all other cases, return false.
  79
  80 locate_matrix_entry(m, r1, c1, r2, c2, fn, rel)
  81
  82    The first argument must be a matrix; the arguments
  83    r1 through c2 determine a sub-matrix of m that consists of
  84    rows r1 through r2 and columns c1 through c2.
  85
  86    Find a entry in the sub-matrix m that satisfies some property.
  87    Three cases:
  88
  89      (1) rel = bool and fn a predicate:
  90
  91          Scan the sub-matrix from left to right then top to bottom,
  92          and return the index of the first entry that satisfies the
  93          predicate fn. If no matrix entry satisfies fn, return false.
  94
  95      (2) rel = 'max and fn real-valued:
  96
  97          Scan the sub-matrix looking for an entry that maximizes fn.
  98          Return the index of a maximizing entry.
  99
 100      (3) rel = 'min and fn real-valued:
 101
 102          Scan the sub-matrix looking for an entry that minimizes fn.
 103          Return the index of a minimizing entry.
 104
 105 lu_backsub(m, b fld)
 106
 107   When m = lu_factor(A), then lu_backsub(m,b, fld) solves the linear
 108   systme A x = b.
 109
 110 lu_factor(m, generalfield)
 111
 112   Return a list of the form [LU, perm, {cnd}], where
 113
 114   (1) The matrix LU contains the factorization of m in a packed form. Packed
 115       form means three things: First, the rows of LU are permuted according to the
 116       list perm.  If, for example, perm is the list [3,2,1], the actual first row
 117       of the LU factorization is the third row of the matrix LU. Second,
 118       the lower triangular factor of m is the lower triangular part of LU with the
 119       diagonal entries replaced by all ones. Third, the upper triangular factor of
 120       m is the upper triangular part of LU.
 121
 122   (2) When the field is either 'floatfield' or 'complexfloatfield,'
 123       the number 'cnd' is an upper bound for the infinity norm condition number of m.
 124       For all other fields, the condition number isn't estimated. For these fields,
 125       lu_factor returns a two item list.
 126
 127   The argument m must be a square matrix.
 128
 129   The argument 'fld' must be a symbol that determines a field. The pre-defined
 130   fields are:
 131
 132     (a) generalfield -- the field of Maxima expressions,
 133     (b) floatfield --  the field of floating point numbers of the type double,
 134     (c) complexfloatfield --  the field of complex floating point numbers of the
 135         type double,
 136     (d) crefield  -- the field of Maxima CRE expressions,
 137     (e) rationalfiled -- the field of rational numbers.
 138
 139   When the field is floatfield or complexfloatfield, the algorithm uses partial
 140   pivoting; when the field is generalfield, rows are switched only when needed to
 141   avoid a zero pivot.
 142
 143   Floating point addition arithmetic isn't associative, so the meaning of 'field'
 144   differs from the mathematical definition.
 145
 146   There is no user-interface for defining a new field. A user that is familiar with
 147   Common Lisp should be able to define a new field.
 148
 149   To compute the factorization, the first task is to convert each matrix entry
 150   to a member of the indicated field. When conversion isn't possible, the factorization
 151   halts with an error message. Members of the field needn't be Maxima expressions.
 152   Members of the complexfloatfield, for example, are Common Lisp complex numbers. Thus
 153   after computing the factorization, the matrix entries must be converted to Maxima
 154   expressions.
 155
 156   See also get_lu_factors.
 157
 158   Examples
 159
 160    (%i1) w[i,j] := ?random(1.0) + %i * ?random(1.0)$
 161            Evaluation took 0.00 seconds (0.00 elapsed)
 162    (%i2) m : genmatrix(w,100,100)$
 163            Evaluation took 00.30 seconds (00.30 elapsed)
 164    (%i3) lu_factor(m,complexfloatfield)$
 165            Evaluation took 1.45 seconds (1.45 elapsed)
 166    (%i4) lu_factor(m,generalfield)$
 167            Evaluation took 12.10 seconds (12.10 elapsed)
 168
 169    (%i1) m : matrix([1-z,3],[3,8-z]);
 170    (%o1) matrix([1-z,3],[3,8-z])
 171    (%i2) lu_factor(m,generalfield);
 172    (%o2) [matrix([1-z,3],[3/(1-z),-z-9/(1-z)+8]),[1,2]]
 173    (%i3) get_lu_factors(%);
 174    (%o3) [matrix([1,0],[0,1]),matrix([1,0],[3/(1-z),1]),matrix([1-z,3],[0,-z-9/(1-z)+8])]
 175    (%i4) %[1] . %[2] . %[3];
 176    (%o4) matrix([1-z,3],[3,8-z])
 177
 178 mat_norm(m, p)
 179
 180   Return the matrix p-norm of the matrix m.  The allowed values for p are
 181   1, inf, and frobenius (the Frobenius matrix norm). The matrix m should be
 182   an unblocked matrix.
 183
 184 matrixp (e, {pred})
 185
 186   Given an optional argument pred, return true if e is
 187   a matrix and pred evaluates to true for every matrix element.
 188   When matrixp is not given an optional argument, return true
 189   if e is a matrix.  In all other cases, return false.
 190
 191 nonnegintegerp (n)
 192
 193   Return true if and only if n >= 0 and n is an integer.
 194
 195 nullspace(m)
 196
 197    If m is a matrix, return span(v1,v2,...,vn), where the set {v1,v2,...,vn}
 198    is a basis for the nullspace of m.  The span of the empty set is {0}.
 199    Thus, when the nullspace has only one member, return span().
 200
 201 nullity(m)
 202
 203    If m is a matrix, return the dimension of the nullspace of m.
 204
 205 orthogonal_complement(v1,v2,...,vn)
 206
 207    Return span(u1,u2,...,um), where the set {u1,u2,...,um} is a
 208    basis for the orthogonal complement of the set(v1,v2,...,vn).
 209
 210    Each vector v1 through vn must be a n x 1 matrix.
 211
 212 polynomialp(p, l,  {(coeffp 'constantp), (exponp 'nonnegintegerp)})
 213
 214   Return true if p is a polynomial in the variables in the Maxmia list l,
 215   The predicate coeffp (default constantp) must evaluate to true for each
 216   coefficient, and the predicate exponp must evaluate to true for all
 217   exponents of the variables in the list l. If you want to use a non-default
 218   value for exponp, you must supply coeffp with a value even if you want
 219   to use the default for coeffp.
 220
 221   The polynomial needn't be expanded:
 222
 223     (%i1) polynomialp((x+1)*(x+2),[x]);
 224     (%o1) true
 225
 226     (%i2) polynomialp((x+1)*(x+2)^a,[x]);
 227     (%o2) false
 228
 229   An example using non-default values for coeffp and exponp:
 230
 231     (%i3) polynomialp((x+1)*(x+2)^(3/2),[x],numberp,numberp);
 232     (%o3) true
 233
 234     (%i4) polynomialp((x^(1/2)+1)*(x+2)^(3/2),[x],numberp,numberp);
 235     (%o4) true
 236
 237   Polynomials with two variables:
 238
 239     (%i5) polynomialp(x^2 + 5*x*y + y^2,[x]);
 240     (%o5) false
 241
 242     (%i6) polynomialp(x^2 + 5*x*y + y^2,[x,y]);
 243     (%o6) true
 244
 245 polytocompanion(p,x)
 246
 247   If p is a polynomial in x, return the companion matrix of p. For
 248   a monic polynomial p of degree n, we have p = (-1)^n charpoly(polytocompanion(p, x)).
 249
 250   When p isn't a polynomial in x, signal an error.
 251
 252 ptriangularize(m,v)
 253
 254    If m is a matrix with each entry a polynomial in v, return
 255    a matrix m' such that
 256
 257      (1) m' is upper triangular,
 258      (2) m' = E_n ... E_1 m, where E1 through En are elementary matrices
 259          whose entries are polynomials in v,
 260      (3) |det(m)| = |det(m')|,
 261
 262    Note: This function doesn't check that every entry is a polynomial in v.
 263
 264 rowop(m,i,j,theta)
 265
 266   If m is a matrix, return the matrix that results from doing the
 267   row operation Ri <- Ri - theta * Rj. If m doesn't have a row
 268   i or j, signal an error.
 269
 270 rowswap(m,i,j)
 271
 272   If m is a matrix, swap rows i and j. If m doesn't have a row
 273   i or j, signal an error.
 274
 275 mat_unblocker(m)
 276
 277   If m is a block matrix, return its unblocked form. If m is a matrix,
 278   return m; otherwise, signal an error.  If you use block matrices,
 279   most likely you'll want to set matrix_element_mult : "." and
 280   matrix_element_transpose : 'transpose.
 281
 282   Example
 283
 284     (%i1) a : matrix([1,2],[3,4])$
 285     (%i2) b : matrix([7,8],[9,10])$
 286     (%i3) matrix([a,b]);
 287     (%o3) matrix([matrix([1,2],[3,4]),matrix([7,8],[9,10])])
 288     (%i4) mat_unblocker(%);
 289     (%o4) matrix([1,2,7,8],[3,4,9,10])
 290
 291 vandermonde_matrix([x1,x2,...,xn])
 292
 293   Return a n by n matrix whose i-th row is
 294   [1, xi, xi^2, ... xi^(n-1)].
 295
 296 zeromatrixp(m)
 297
 298   If is(equal(e,0)) is true for each entry of the matrix m, return true; otherwise,
 299   return false.  When m is not a matrix, signal an error.
 300
 301
 302 Short demo:
 303
 304 (%i1) batch("linalg.demo");
 305 batching #pc:/maxima/linearalgebra/linalg.demo
 306 (%i2) load(linalg)
 307 (%o2) c:/maxima/linearalgebra/linalg.mac
 308 (%i3) m:matrix([1,2],[1,2])
 309 (%o3) matrix([1,2],[1,2])
 310 (%i4) nullspace(m)
 311 (%o4) span(matrix([1],[-1/2]))
 312 (%i5) columnspace(m)
 313 (%o5) span(matrix([1],[1]))
 314 (%i6) ptriangularize(m-z*ident(2),z)
 315 (%o6) matrix([1,2-z],[0,3*z-z^2])
 316 (%i7) m:matrix([1,2,3],[4,5,6],[7,8,9])-z*ident(3)
 317 (%o7) matrix([1-z,2,3],[4,5-z,6],[7,8,9-z])
 318 (%i8) mm:ptriangularize(m,z)
 319 (%o8) matrix([4,5-z,6],[0,66/49,-z^2/7+(102*z)/49+132/49],[0,0,(49*z^3)/264-(245*z^2)/88-(147*z)/44])
 320 (%i9) algebraic:true
 321 (%o9) true
 322 (%i10) tellrat(mm[3,3])
 323 (%o10) [z^3-15*z^2-18*z]
 324 (%i11) mm:ratsimp(mm)
 325 (%o11) matrix([4,5-z,6],[0,66/49,-(7*z^2-102*z-132)/49],[0,0,0])
 326 (%i12) nullspace(mm)
 327 (%o12) span(matrix([1],[(z^2-14*z-16)/8],[-(z^2-18*z-12)/12]))
 328 (%i13) m:matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16])
 329 (%o13) matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16])
 330 (%i14) column_space(m)
 331 (%o14) column_space(matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]))
 332 (%i15) apply('orthogonal_complment,args(nullspace(transpose(m))))
 333 (%o15) orthogonal_complment(matrix([0],[1],[-2],[1]),matrix([1],[0],[-3],[2]))
 334
 335