share/linearalgebra/matrixexp.lisp

   1 ;;  Copyright 2004 by Barton Willis
   2
   3 ;;  This is free software; you can redistribute it and/or
   4 ;;  modify it under the terms of the GNU General Public License,
   5 ;;  http://www.gnu.org/copyleft/gpl.html.
   6
   7 ;; This software has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 ;; Innocent looking problems, such as matrixexp(matrix([x,1,0],[1,1,1],[0,1,1])),
  11 ;; generate huge (and most likely worthless) expressions. These huge
  12 ;; expressions strain Maxima's rational function code; to avoid errors
  13 ;; such as "..quotient by polynomial of higher degree" I found it necessary to
  14 ;;
  15 ;;  (1) set gcd == spmod and algebraic == true,
  16
  17 ;;  (2) always give ratsimp and fullratsimp a value (even if nil) for its
  18 ;;      optional second argument,
  19
  20 ;;  (3) set ratvars to an empty list at the start of most functions.
  21
  22 ;; I didn't try to find the real cause for these bugs.  The function
  23 ;; spectral_rep does check its output. Although these checks are slow,
  24 ;; I recommend that they stay.
  25
  26 ;; Map the CL function f over the elements of a Maxima matrix mat and
  27 ;; return a Maxima matrix. The first argument should be a CL function.
  28
  29 ($put '$matrixexp 1 '$version)
  30
  31 ;; When mat is a square matrix, return exp(mat * x). The second
  32 ;; argument is optional and it defaults to 1.
  33
  34 (defun $matrixexp (mat &optional (x 1))
  35   (let ((sp) (d) (p) (id) (n ($length ($args mat))) (f))
  36     ($ratvars)
  37     ($require_square_matrix mat '$first '$matrixexp)
  38     (setq mat ($spectral_rep mat))
  39     (setq sp ($first mat))
  40     (setq p ($second mat))
  41     (setq sp (cons '(mlist simp)
  42                    (mapcar #'(lambda (s) ($exp (mult s x))) (cdr sp))))
  43     (setq d (mult x ($third mat)))
  44     (setq id ($ident n))
  45     (setq f id)
  46     (setq n (+ n 1))
  47     (dotimes (i n)
  48       (setq f (add id (div (ncmul2 d f) (- n i)))))
  49     ($fullratsimp (ncmul2 (ncmul2 sp p) f) nil)))
  50
  51 ;; Let f(var) = expr.  This function returns f(mat), where 'mat' is a
  52 ;; square matrix.  Here expr is an expression---it isn't a function!
  53
  54 (defun require-lambda (e n pos fun-name)
  55   (let ((var))
  56     (if (and (consp e) (consp (car e)) (eq 'lambda (mop e)) (= 3 (length e))
  57              ($listp (nth 1 e)) (setq var (cdr (nth 1 e))) (= n (length var))
  58              (every #'(lambda (s) (or (symbolp s) ($subvarp s))) var))
  59         (list var (nth 2 e))
  60       (merror "The ~:M argument to `~:M' must be a lambda form with ~:M variable(s)" pos fun-name n))))
  61
  62 (defun $matrixfun (lamexpr mat)
  63   (let ((z (gensym)) (expr) (var) (sp) (d) (p) (di)
  64         (n ($length ($args mat))) (f 0))
  65
  66     ($require_square_matrix mat '$second '$matrixexp)
  67     (setq expr (require-lambda lamexpr 1 '$first '$matrixfun))
  68     (setq var (nth 0 (nth 0 expr)))
  69     (setq expr (nth 1 expr))
  70     ($ratvars)
  71     (setq expr ($substitute z var expr))
  72     (setq mat ($spectral_rep mat))
  73     (setq sp ($first mat))
  74     (setq p ($second mat))
  75     (setq d ($third mat))
  76     (setq di ($ident n))
  77     (setq sp (cdr sp))
  78     (dotimes (i (+ n 1))
  79       (setq f (add
  80                f
  81                (ncmul2 di (ncmul2
  82                            (cons '(mlist simp)
  83                                  (mapcar #'(lambda (s)
  84                                              ($substitute s z expr)) sp)) p))))
  85       (setq di (ncmul2 di d))
  86       (setq expr (div ($diff expr z) (factorial (+ i 1)))))
  87     ($fullratsimp (simplify f) nil)))
  88
  89 ;; Return the residue of the rational expression e with respect to the
  90 ;; variable var at the point pt.  Assumptions:
  91
  92 ;;  (1) the denominator of e divides ker,
  93 ;;  (2) e is a rational expression,
  94 ;;  (3) ker is a polynomial,
  95 ;;  (4) pt is z zero of ker and ord is its order.
  96
  97 (defun rational-residue (e var pt ker ord)
  98   (let (($gcd '$spmod) ($algebraic t) ($ratfac nil) (p) (q) (f (sub var pt)))
  99     ($ratvars)
 100     (setq e ($fullratsimp e var))
 101     (setq p ($num e))
 102     (setq q ($denom e))
 103     (setq p (mult p ($quotient ker q var)))
 104     (setq e ($fullratsimp (div p ($quotient ker (power f ord) var)) var))
 105     ($fullratsimp
 106      ($substitute pt var (div ($diff e var (- ord 1)) (factorial (- ord 1))))
 107      nil)))
 108
 109
 110 (defun $spectral_rep (mat)
 111   ($require_square_matrix mat '$first '$spectral_rep)
 112   (let (($gcd '$spmod) ($algebraic t) ($resultant '$subres) (ord) (zi)
 113         ($ratfac nil) (z (gensym)) (res) (m) (n ($length ($args mat)))
 114         (p) (p1) (p2) (sp) (proj))
 115
 116     ($ratvars)
 117     (setq p ($newdet (sub mat (mult z ($ident n)))))
 118     (if (oddp n) (setq p (mult -1 p)))
 119
 120     ;; p1 = (z - z1)(z - z2) ... (z - zk), where z1 thru zk are
 121     ;; the distinct zeros of p.
 122
 123     (setq p1 ($first ($divide p ($gcd p ($diff p z) z) z)))
 124     (setq p2 ($resultant p1 ($diff p1 z) z))
 125
 126     (cond ((and (not ($constantp p2)) (not (like 0 p2)))
 127            ($ratvars)
 128            (setq p2 ($sqfr p2))
 129            (if (mminusp p2) (setq p2 (mult -1 p2)))
 130            (setq p2 (if (mtimesp p2) (margs p2) (list p2)))
 131            (setq p2 (mapcar #'(lambda (s) (if (mexptp s) (nth 1 s) s)) p2))
 132            (setq p2 `((mtimes simp) ,@p2))
 133         (mtell "Proviso: assuming ~:M ~%" `((mnotequal simp) ,p2 0))))
 134
 135     (setq sp ($solve p z))
 136     (setq sp (mapcar '$rhs (cdr sp)))
 137     (cond ((not (eq n (apply #'+ (cdr $multiplicities))))
 138            (print `(ratvars = ,$ratvars gcd = '$gcd algebraic = ,$algebraic))
 139            (print `(ratfac = ,$ratfac))
 140            (merror "Unable to find the spectrum")))
 141
 142     (setq res ($fullratsimp (ncpower (sub (mult z ($ident n)) mat) -1) z))
 143     (setq m (length sp))
 144     (dotimes (i m)
 145       (setq zi (nth i sp))
 146       (setq ord (nth (+ i 1) $multiplicities))
 147       (push (matrix-map #'(lambda (e) (rational-residue e z zi p ord)) res)
 148             proj))
 149
 150     (setq proj (nreverse proj))
 151     (setq m (length proj))
 152     (dotimes (i m)
 153       (setq mat (sub mat (mult (nth i sp) (nth i proj)))))
 154
 155     (cond ((check-spectral-rep proj n)
 156            (push `(mlist simp) proj)
 157            (push `(mlist simp) sp)
 158            `((mlist simp) ,sp ,proj, ($fullratsimp mat nil)))
 159           (t
 160            (merror "Unable to find the spectral representation")))))
 161
 162 (defun check-spectral-rep (proj n)
 163   (let* ((m (length proj)) (okay) (zip ($zeromatrix n n)) (qi)
 164          ($gcd '$spmod) ($algebraic t) ($ratfac nil))
 165     ($ratvars)
 166
 167     (setq proj (mapcar #'(lambda (s) ($fullratsimp s nil)) proj))
 168     (setq okay (like ($ident n) ($fullratsimp (apply 'add proj) nil)))
 169
 170     (dotimes (i m)
 171       (setq qi (nth i proj))
 172       (setq qi ($fullratsimp (sub qi (ncmul2 qi qi)) nil))
 173       (setq okay (and okay (like zip qi))))
 174
 175     (dotimes (i m)
 176       (setq qi (nth i proj))
 177       (dotimes (j i)
 178         (setq okay (and okay (like zip ($fullratsimp
 179                                         (ncmul2 qi (nth j proj)) nil))))
 180         (setq okay (and okay (like zip ($fullratsimp
 181                                         (ncmul2 (nth j proj) qi) nil))))))
 182     okay))
 183
 184
 185
 186
 187
 188
 189
 190
 191