share/contrib/diag.mac

   1 /*
   2    diag.mac
   3
   4    Contains the functions:
   5      (1) diag
   6      (2) JF
   7      (3) jordan
   8      (4) minimalPoly
   9      (5) dispJordan
  10      (6) ModeMatrix
  11      (7) mat_function
  12
  13    plus some internal utility funcions. If you make use of one of
  14    them, consider asking on the mailing list: we can try to make a
  15    more useful API and give it a sensible name (not prefixed by
  16    "diag_"), that we won't change.
  17
  18    See the manual or the comments above the relevant function for more
  19    information on how each works.
  20
  21    As with the rest of Maxima, this code is distributed under the GPL,
  22    version 2+
  23 */
  24 load("eigen")$
  25
  26 /*
  27   If EXPR is a matrix then return it unchanged. Otherwise, return a
  28   1x1 matrix whose sole element is EXPR.
  29 */
  30 diag_matrixify (expr) :=
  31   if matrixp (expr) then expr else matrix ([expr])$
  32
  33 /*
  34   Pad R with zeros before and after to make it WIDTH wide using INDENT
  35   zeros on the left hand side.
  36 */
  37 diag_zeropad_row (r, indent, width) :=
  38   block([right: width - (left + length(r))],
  39     if is(right < 0) then
  40       error(r, "does not fit in a row of length", width,
  41             "with indent", indent),
  42     append (makelist (0, indent), r, makelist (0, right)))$
  43
  44 /*
  45   Return the direct sum of the elements of LST, made into 1x1 matrices
  46   if they weren't already matrices.
  47 */
  48 diag (lst) :=
  49   block ([lst: map (diag_matrixify, lst),
  50           width, height, left: 0, rows: []],
  51     width: lsum (length (first (A)), A, lst),
  52     height: lsum (length (A), A, lst),
  53     for A in lst do (
  54       for r in args(A) do
  55         rows: cons (diag_zeropad_row (r, left, width), rows),
  56       left: left + length (first (A))),
  57     apply (matrix, reverse (rows)))$
  58
  59 /*
  60   Construct a Jordan block of size N x N and eigenvalue EIVAL.
  61 */
  62 JF (eival, n) :=
  63   if not (numberp (n)) then
  64     'JF (eival, n)
  65   elseif not (integerp (n)) then
  66     error (n, "is not an integer, so not a valid Jordan block size")
  67   elseif is (n <= 0) then
  68     error ("Cannot construct matrix with negative size", n)
  69   else
  70     genmatrix (
  71       lambda ([i,j],
  72               if is (i = j) then eival
  73               elseif is (j = i + 1) then 1
  74               else 0),
  75       n, n)$
  76
  77 /*
  78   Calculate the correct partition of multiplicity for the matrix A at
  79   the eigenvalue eival.
  80 */
  81 diag_calculate_mult_partition (A, multiplicity, eival) :=
  82   if is (multiplicity = 1) then [1] else
  83   block ([Abar: A - eival * ident (length (A)),
  84           n: length (A),
  85           nullity],
  86     nullity: n - rank (Abar),
  87     if is (nullity = 1) then [multiplicity]
  88     elseif is (nullity = multiplicity) then makelist (1, nullity)
  89     else block ([blocks_left: nullity, mults: [], Abarpow: Abar, dnull],
  90       for k: 1 do (
  91         Abarpow: Abarpow . Abar,
  92         dnull: n - rank (Abarpow) - nullity,
  93         nullity: nullity + dnull,
  94         if (dnull < blocks_left) then (
  95           mults: append (makelist (k, blocks_left - dnull), mults),
  96           multiplicity: multiplicity - k * (blocks_left - dnull),
  97           blocks_left: dnull),
  98         if is (blocks_left * (k + 1) = multiplicity) then
  99           return (append (makelist (k+1, blocks_left), mults)),
 100         if is (blocks_left * (k + 1) > multiplicity) then
 101           error ("Unexpected blocks left over!"))))$
 102
 103 /*
 104   Calculate the JCF of A. Returns a list of eigenvalues and their
 105   multiplicities.
 106 */
 107 jordan (A) :=
 108   if not (matrixp (A)) then
 109     'jordan(A)
 110   else block([eigenlist: eigenvectors (A)],
 111     map (lambda ([eival, mult],
 112                  cons (eival, diag_calculate_mult_partition (A, mult, eival))),
 113          eigenlist[1][1], eigenlist[1][2]))$
 114
 115 /*
 116   A simple sanity check that arguments are of the form returned by
 117   jordan().
 118 */
 119 diag_jordan_info_check (lst) :=
 120   if not (listp (lst)) then false
 121   else block([ret: true],
 122     for pair in lst do
 123       if not (listp (pair)) then (ret: false, return (false))
 124       elseif is (length (pair) < 2) then
 125         error ("Found a surprisingly short jordan_info list:",
 126                pair),
 127     ret)$
 128
 129 /*
 130   Calculate the minimal polynomial, expressed in the symbol "x", of
 131   any matrix with the given list of eigenvalues and their
 132   multiplicities.
 133
 134   NOTE: This assumes that the first multiplicity for an eigenvalue is
 135         the largest and that the given eigenvalues are distinct. This
 136         is the case for the output of jordan().
 137 */
 138 minimalPoly (jordan_info) :=
 139   if not (diag_jordan_info_check (jordan_info)) then
 140     'minimalPoly (jordan_info)
 141   else
 142   lreduce ("*",
 143     map (lambda ([eival_lst], ('x - first (eival_lst))^(second (eival_lst))),
 144          jordan_info))$
 145
 146 /*
 147   Take a list of eigenvalues and their multiplicities and build the
 148   corresponding Jordan matrix.
 149 */
 150 dispJordan (jordan_info) := block ([blocks: []],
 151   if not (diag_jordan_info_check (jordan_info)) then
 152     'dispJordan (jordan_info)
 153   else
 154   for eival_lst in jordan_info do
 155     for size in rest (eival_lst) do
 156       blocks: cons (JF (first (eival_lst), size), blocks),
 157   diag (reverse (blocks)))$
 158
 159 /*
 160   Takes a "sorted list" and produces a histogram of elements and their
 161   frequencies:
 162
 163      [3,3,2,2,2,1,1,1] => [[3,2],[2,3],[1,3]]
 164
 165   The "sorted" requirement is merely that any object only occurs in
 166   one block. Thus, we will do the right thing on [1,1,0,0,2,2], but
 167   not on [1,1,0,0,1,1].
 168 */
 169 diag_sorted_list_histogram (list) :=
 170   if is (emptyp (list)) then [] else
 171   block([pairs: [], current: first (list), count: 1],
 172     for lst: rest (list) next rest (lst) while not (emptyp (lst))
 173     do
 174       if is (current = first (lst)) then
 175         count: count + 1
 176       else (
 177         pairs: cons ([current, count], pairs),
 178         current: first (lst),
 179         count: 1),
 180     reverse (cons ([current, count], pairs)))$
 181
 182 /*
 183   Find a formula for a general element of the kernel of A.
 184
 185   Returns the solution as a row vector of (linear) expressions in the
 186   variables of %rnum_list.
 187 */
 188 diag_kernel_element (A) :=
 189   block ([vars: makelist (gensym (), length (first (A))), eqns, soln],
 190     if length (vars) = 1 then
 191       eqns: [first(vars) * first (first (A))]
 192     else
 193       eqns: map (first, args (A . apply (matrix, map("[", vars)))),
 194     soln: algsys (eqns, vars),
 195     if emptyp (soln) then error ("Could not solve for the kernel of ", A),
 196     map (rhs, first (soln)))$
 197
 198 /*
 199   Calculate the general form for a chain of generalized eigenvectors
 200   for A, of the form [v1, v2, ..., vn] such that (A-eival*I).v1 = 0 and
 201   (A-eival*I).vk = v(k-1) where n = degree.
 202
 203   If degree is not correct, throws an error.
 204
 205   Otherwise returns [chain, freevars] where chain is a list of the
 206   eigenvectors, each of which is represented as a list, and freevars
 207   is the list of free variables.
 208 */
 209 diag_general_jordan_chain (A, eival, degree) :=
 210   block ([n: length (A), Abar, ret],
 211     if (is (n < 1) or not (length (A[1]) = n)) then
 212     error ("Invalid input matrix: ", A),
 213     Abar: A - eival * ident (n),
 214     ret: [diag_kernel_element (Abar^^degree)],
 215     for i : 2 thru degree do
 216       ret: cons (args (map (first, Abar . first(ret))), ret),
 217     [ret, %rnum_list])$
 218
 219 /*
 220   LI_ROWS should be a list of lists representing row vectors that are
 221   linearly independent. Then CHAIN_EXPR is a list of row vectors in
 222   VARIABLES representing a Jordan chain. The function tries to find a
 223   choice of variables such that CHAIN_EXPR is linearly independent of
 224   LI_ROWS.
 225
 226   (Note: Because chain_expr represents a Jordan chain, it suffices to
 227          check linear independence for the first and last elements of
 228          CHAIN_EXPR)
 229 */
 230 diag_find_li_chain (li_rows, chain_expr, variables) :=
 231   block ([dict: false, extra_rows],
 232     for vals in ident (length (variables)) do (
 233       dict: map ("=", variables, vals),
 234       extra_rows: subst (dict,
 235                          if is (length (chain_expr) > 1) then
 236                          [first (chain_expr), last (chain_expr)]
 237                          else [first (chain_expr)]),
 238       if is (rank (apply (matrix, append (li_rows, extra_rows))) =
 239              length (li_rows) + length (extra_rows))
 240       then return (true)
 241       else dict: false),
 242     if is (dict = false) then
 243       error ("Could not find a linearly independent vector"),
 244     /* Rectform since complex eigenvalues result in horrible
 245        expressions here otherwise */
 246     rectform (subst (dict, chain_expr)))$
 247
 248 /*
 249    Calculates the mode matrix of A, which is a matrix T such that
 250       T^^(-1) . A . T = JordanForm (A)
 251
 252    F should be a list of the eigenvalues of A together with their
 253    multiplicities, as returned by jordan(A).
 254
 255    Notes on algorithm:
 256
 257    rest(ev_lst) is a list of chain lengths, which we convert to
 258    multiplicities. To ensure we get the correct answer, when we hunt
 259    for chains we need to start with those of maximal
 260    length. Fortunately, jordan() returns lists where rest(ev_lst) is
 261    (weakly) decreasing. As such, we can just work down the lengths in
 262    the given order and we'll get the right answer.
 263
 264    For each chain length, we find a general form for a Jordan chain
 265    that length then we evaluate the formula, varying the parameters to
 266    make sure we end up with a linearly independent chain. Note that it
 267    suffices to check that the first term of the chain list (the actual
 268    eigenvector) is linearly independent from the rows we've got so
 269    far.
 270
 271    Since generalised eigenvectors for different eigenvalues are
 272    linearly independent, we don't bother checking there.
 273 */
 274 diag_mode_matrix (a, F) :=
 275   if not (diag_jordan_info_check (F)) or not (matrixp (a)) then
 276     'ModeMatrix (a, F)
 277   else block([msize: length(a), all_rows: []],
 278     for ev_pair in F do
 279       block ([eival: first (ev_pair),
 280               multiplist: diag_sorted_list_histogram (rest (ev_pair)),
 281               mat_rows: [], mat_rank: 0],
 282         for degree_pair in multiplist do
 283           block ([mindeg: first (degree_pair),
 284                   genev_pair, genevs, free_vars],
 285             genev_pair: diag_general_jordan_chain (a, eival,
 286                                                    first (degree_pair)),
 287             for k : 1 thru second (degree_pair) do
 288               mat_rows:
 289                 append (mat_rows,
 290                         diag_find_li_chain (mat_rows,
 291                                             first (genev_pair),
 292                                             second (genev_pair)))),
 293         all_rows: append (all_rows, mat_rows)),
 294     transpose (apply (matrix, all_rows)))$
 295
 296 /*
 297   Finds a matrix T such that T^^(-1) . A . T is the Jordan form of A.
 298 */
 299 ModeMatrix (A, [jordan_info]) :=
 300   if is (length (jordan_info) > 1) then
 301     error ("Too many arguments for ModeMatrix. (Expects 1 or 2)")
 302   elseif not (matrixp (A)) then
 303     if is (length (jordan_info) = 1) then
 304       'ModeMatrix (A, first (jordan_info))
 305     else
 306       'ModeMatrix (A)
 307   else
 308     diag_mode_matrix (A,
 309                       if is (length (jordan_info) = 1)
 310                       then first (jordan_info) else jordan (A))$
 311
 312 /*
 313   Return a list of taylor coefficients for the function f expanded
 314   around some arbitrary point VAR where f(var) = EXPR. The first
 315   element of the list is f(var) and the last is
 316   diff(f(var),var,maxpow).
 317 */
 318 diag_taylor_coefficients (expr, var, maxpow) :=
 319   block ([coeffs: [expr]],
 320     for k:1 thru maxpow do
 321       (expr: diff (expr, var), coeffs: cons (expr/k!, coeffs)),
 322     reverse (coeffs))$
 323
 324 /*
 325   Build a matrix of the form:
 326
 327         [ f(0)  f'(0)  f''(0) ]
 328         [  0    f(0)   f'(0)  ]
 329         [  0      0     f(0)  ]
 330
 331   except with coefficients taken from the COEFFS list, substituting
 332   EIGENVALUE for VAR, where SIZE is the size of the matrix. This is
 333   the result of Taylor expanding a f(A), where A is a Jordan block,
 334   around its eigenvalue.
 335 */
 336 diag_taylor_expand_block (coeffs, var, eigenvalue, size) := (
 337   coeffs: makelist (subst (eigenvalue, var, coeffs[i]), i, 1, size),
 338   apply (matrix,
 339          makelist (makelist (if is (i<k) then 0 else coeffs[i-k+1], i, 1, size),
 340                    k, 1, size)))$
 341
 342 /*
 343   Calculate the value of an analytic function on the matrix
 344   represented by the Jordan list JORDAN. The function is given by EXPR
 345   in VAR.
 346 */
 347 diag_mat_function_jordan (jordan, expr, var) :=
 348   block ([coeffs, blocks: [],
 349           max_degree: lmax (map (second, jordan)) - 1],
 350     /*
 351       Expand f(var) about some arbitrary point as a Taylor series. The
 352       Jordan matrix is diagonal plus a nilpotent matrix order one less
 353       than the largest block. We need coefficients the same order as
 354       that maximum block.
 355     */
 356     coeffs: diag_taylor_coefficients (expr, var, max_degree),
 357     /*
 358       We calculate the value of EXPR on each Jordan block using
 359       DIAG_TAYLOR_EXPAND_BLOCK. The degrees in JORDAN_LST are known to
 360       be decreasing, so we can be slightly clever about not computing
 361       things repeatedly.
 362     */
 363     for jordan_lst in jordan do
 364       block ([cached_block:
 365                 diag_taylor_expand_block (coeffs, var,
 366                                           first (jordan_lst),
 367                                           second (jordan_lst)),
 368               cached_size: second(jordan_lst)],
 369         for size in rest (jordan_lst) do
 370           (if is (size # cached_size) then
 371              (cached_size: size,
 372               cached_block:
 373                 apply (matrix,
 374                        makelist (
 375                          makelist (cached_block[i,j], j, 1, size),
 376                          i, 1, size))),
 377            blocks: cons (cached_block, blocks))),
 378     diag (reverse (blocks)))$
 379
 380 /*
 381   Take an analytic function, f, and a matrix, A, and calculate f(A) by
 382   means of the associated Taylor series.
 383 */
 384 mat_function (f, A) :=
 385   if not (matrixp (A)) then
 386     'mat_function (f, A)
 387   else
 388   block ([jj: jordan (A), var: gensym(), modemat],
 389     modemat: ModeMatrix (A, jj),
 390     modemat . diag_mat_function_jordan (jj, apply (f, [var]), var)
 391             . (modemat)^^(-1))$