share/pdiff/pdiff.lisp

   1 ;; A positional derivative package for Maxima.
   2 ;; Copyright (C) 2002, 2008, Barton Willis
   3
   4 ;;  Maxima pdiff is free software; you can redistribute it and/or
   5 ;;  modify it under the terms of the GNU General Public License,
   6
   7 (in-package :maxima)
   8 ($put '$pdiff 1.4 '$version)
   9
  10 ;; Required for wxMaxima display--I switched from ((%pderivop) ...) to (($pderivop) ...).
  11 (setf (get '$pderivop 'wxxml) 'wxxml-pderivop)
  12
  13 ;; When use_pdiff is true, use positional derivatives for unknown
  14 ;; non-subscripted functions. By unknown function, I mean a function that
  15 ;; is not bound to a formula and that has a derivative that is not known
  16 ;; to Maxima. Let's default $use_pdiff to false.
  17
  18 (defmvar $use_pdiff nil)
  19 (setf $use_pdiff t)
  20
  21 (defmvar $commute_partial_derivatives t)
  22
  23 ;; Return the number of arguments (arity) for f; when the arity is
  24 ;; unknown or the arity is variable (for example addition), return false.
  25 ;; To attempt to find the arity, we check if f is a user-defined function
  26 ;; (examine the $functions list), check if f has a gradef property, check
  27 ;; if f is a trig, inverse trig, log or abs function.
  28
  29 ;; Subscripted functions aren't members of $functions; for this and
  30 ;; other reasons, we disallow subscripted functions from being
  31 ;; positionally differentiated.
  32
  33 (defun get-number-args (f)
  34   (let ((z (cdar (member f (cdr $functions) :key #'caar))))
  35     (cond ((and z (every #'symbolp z)) (length z))
  36           ((and (atom f) (or (trigp f) (arcp f) (memq f '(%log %gamma %erf %sqrt mabs)))) 1)
  37           ((memq f '($bessel_i $bessel_j $bessel_k $bessel_y)) 2)
  38           ((or (equal f "^") (equal f ".") (equal f "^^")) 2)
  39           ((and (op-equalp f 'lambda) (every #'symbolp (margs (second f))))
  40            ($length (second f)))
  41           ((and (symbolp f) (get f 'grad)) (length (first (get f 'grad))))
  42           ((eq t ($constantp f)) 'constant-function)
  43           (t nil))))
  44
  45 ;; pderivop(f,n1,n2,...np) returns the function whose formula is
  46 ;; (x1,x2,...,xn) --> diff(f(x1,x2,..xn),x1,n1,x2,n2,...). The
  47 ;; pderivop function allows the user to do things like
  48 ;;
  49 ;;  (c1) tellsimpafter(pderivop(f,1)(a),1);
  50 ;;  (c2) tellsimpafter(pderivop(f,2)(a),1);
  51 ;;  (c3) diff(f(x),x,2) + diff(f(x),x);
  52 ;;
  53 ;;  (d3)                               f   (x) + f   (x)
  54 ;;                                      (2)       (1)
  55 ;;  (c4) subst(a,x,%);
  56 ;;  (d4)      2
  57
  58 ;; And
  59 ;;   (c5) g(x) := x^2;
  60 ;;   (c6) pderivop(g,1);
  61 ;;   (c7) ev(%);
  62 ;;   (d7) lambda([i54924], 2 i54924)
  63 ;;   (c8) apply(%,[z]);
  64 ;;
  65 ;;   (d8)     2 z
  66
  67 (defprop $pderivop simp-pderivop operators)
  68
  69 (defun simp-pderivop (x yy z)
  70   (declare (ignore yy))
  71   (let ((f (simplifya (cadr x) z))
  72         (n (mapcar #'(lambda (s) (simplifya s z)) (cddr x)))
  73         (gen-args nil)
  74         (nbr-args))
  75     (setq nbr-args (get-number-args f))
  76     (cond ((eq nbr-args 'constant-function)
  77            (cond ((some #'(lambda (s) (and (integerp s) (> s 0))) n)
  78                   `((lambda simp) ((mlist) ((mlist) $x)) 0))
  79
  80                  ((every #'(lambda (s) (and (integerp s) (= s 0))) n)
  81                   `((lambda simp) ((mlist) ((mlist) $x)) ,f))
  82
  83                  (t `(($pderivop simp) ,f ,@n))))
  84
  85           ;; When every derivative is 0 (this includes the case that there are no derivatives),
  86           ;; return f (even when f is not a function---for example pderivop(a < b) --> a < b).
  87           ((every #'(lambda (s) (equal 0 s)) n) f)
  88
  89           ((integerp nbr-args)
  90
  91            (if (not (= nbr-args (length n)))
  92                (merror "~:M expected ~:M derivative argument(s), but it received ~:M" f nbr-args (length n))
  93              (progn
  94                (dotimes (i nbr-args)
  95                  (push (gensym) gen-args))
  96                (push '(mlist) gen-args)
  97                (setq f ($apply f gen-args))
  98                (setq gen-args (cdr gen-args))
  99                (setq f (apply '$diff (cons f (mapcan #'list gen-args n))))
 100                `((lambda) ((mlist) ,@gen-args) ,f))))
 101
 102           ;; pderivop(pderivop(m1,m2,...mk),n1,n2,...,nk) --> pderivop(m1+n1,m2+n2,...mk+nk),
 103           ;; provided all the m1..nk are explicit nonnegative integers. I could allow declared
 104           ;; integers, I suppose. Also, I could give an error when the number of positional
 105           ;; derivatives aren't equal--for example pderivop(pderivop(f,1,2),3).
 106
 107           ((and (op-equalp f '$pderivop)
 108                 $commute_partial_derivatives
 109                 (every #'(lambda (s) (and (integerp s) (> s -1))) n)
 110                 (every #'(lambda (s) (and (integerp s) (> s -1))) (cdr (margs f)))
 111                 (= (length n) (- (length (margs f)) 1)))
 112            `(($pderivop simp) ,(cadr f) ,@(mapcar #'(lambda (a b) (add a b)) (cddr f) n)))
 113
 114           (t `(($pderivop simp) ,f ,@n)))))
 115
 116 ;; Extension to mapply1:
 117
 118 (setf (get '$pderivop 'mapply1-extension) 'mapply-pdiff)
 119
 120 (defun mapply-pdiff (fn args fnname form)
 121   (declare (ignore fnname form))
 122   (if (and $use_pdiff (eq (mop fn) '$pderivop))
 123       (cond ((equal (length (cddr fn)) (length args))
 124              `((mqapply simp) ,fn ,@args))
 125             (t
 126              (merror "The function ~:M expected ~:M argument(s), but it received ~:M"
 127                      (cadr fn) (length (cddr fn)) (length args))))
 128     nil))
 129
 130
 131 ;; Extension to sdiffgrad
 132
 133 (defun sdiffgrad-pdiff (e x)
 134   (let ((args) (fun (caar e)) (de) (n) (d-order))
 135     (labels ((pderivop (f x n) (simplify `((mqapply) (($pderivop) ,f ,@n) ,@x)))
 136              (incf-ith (i e)
 137                        (let ((k (nth i e)) (q (copy-list e)))
 138                          (setf (nth i q) (add 1 k))
 139                          q))
 140              (i-list (i n) (incf-ith i (make-list n :initial-element 0))))
 141
 142       (cond ((and $use_pdiff (eq fun 'mqapply) (eq (caaadr e) '$pderivop))
 143              (setq args (cddr e))
 144              (setq fun (cadadr e))
 145              (setq de 0)
 146              (setq n (length args))
 147              (setq d-order (cddadr e))
 148              (dotimes (i n de)
 149                (setq de (add de (mul ($diff (nth i args) x) (pderivop fun args (incf-ith i d-order)))))))
 150
 151             ;; We disallow positional derivatives of subscripted functions and lambda forms.
 152
 153             ((and $use_pdiff (null (zl-get fun 'grad)) (not ($subvarp (cadr e)))
 154                   (not (memq fun '(mlessp mleqp mequal mgeqp mgreaterp mcond lambda))))
 155              (setq args (cdr e))
 156              (setq fun (caar e))
 157              (setq de 0)
 158              (setq n (length args))
 159              (dotimes (i n de)
 160                (setq de (add de (mul ($diff (nth i args) x) (pderivop fun args (i-list i n)))))))
 161             (t nil)))))
 162
 163 ;; Display support for positional derivatives. Indicate the derivative order
 164 ;; with a subscript surrounded by parenthesis.
 165
 166 (setf (get '$pderivop 'dimension) 'dimension-pderiv)
 167
 168 (defun dimension-pderiv (form result)
 169   (setq form (cdr form))
 170   (setq form `(( ,(car form) simp array) (("") ,@(cdr form))))
 171   (dimension-array form result))
 172
 173 ;; Extend tex to handle positional derivatives.  Depending on the values of
 174 ;; the option variables $tex_uses_prime_for_derivatives
 175 ;; and $tex_uses_named_subscripts_for_derivatives, derivatives can
 176 ;; tex as superscripted primes, subscripted variable names, or
 177 ;; parenthesis surrounded subscripts that indicate the derivative order.
 178
 179 (defmvar $tex_uses_prime_for_derivatives nil)
 180 (defmvar $tex_prime_limit 3)
 181 (defmvar $tex_uses_named_subscripts_for_derivatives nil)
 182 (defmvar $tex_diff_var_names (list '(mlist) '$x '$y '$z))
 183
 184 (setf (get '$pderivop 'tex) 'tex-pderivop)
 185
 186 ;; Examples
 187
 188 ;; 1. $tex_uses_prime_for_derivatives is true,
 189
 190 ;;    (c1) tex(diff(f(x),x));
 191 ;;           $$f^{\prime}(x)$$
 192 ;;    (c2) tex(diff(f(x),x,2));
 193 ;;           $$f^{\prime\prime}(x)$$
 194 ;;    (c4) tex(diff(f(x),x,4));
 195 ;;           $$f^{(4)}(x)$$
 196
 197 ;; In the last example, the derivative order exceeds $tex_prime_limit,
 198 ;; so the derivative is indicated as shown in (c4).
 199
 200 ;; 2. $tex_uses_named subscripts is true
 201 ;;  (c1) tex_uses_prime_for_derivatives : false;
 202 ;;  (c2) tex(diff(f(x),x));
 203 ;;       $$f_{x}(x)$$
 204 ;;  (c3) tex(diff(f(x),x,2));
 205 ;;       $$f_{xx}(x)$$
 206 ;;  (c4) tex(diff(f(y),y,2));
 207 ;;       $$f_{xx}(y)$$
 208
 209 ;; Although the function argument in (c4) is y, the derivative with
 210 ;; respect to the first argument is indicated by a subscript that is
 211 ;; the first element of tex_diff_var_names. A further example
 212
 213 ;;  (c5) tex_diff_var_names : [\a,\b,\c]$
 214 ;;  (c6) tex(diff(f(x,y,z),x,1,y,1,z,1));
 215 ;;       $$f_{abc}(x,y,z)$$
 216
 217 ;; When the derivative order exceeds tex_prime_limit, we don't use named
 218 ;; subscripts for derivatives; otherwise, we could get ridiculously long
 219 ;; subscripts.
 220
 221 ;;   (c43) tex_prime_limit : 3;
 222 ;;   (c44) tex(diff(f(x,y),x,1,y,1));
 223 ;;         $$f_{xy}(x,y)$$
 224 ;;   (c45) tex(diff(f(x,y),x,1066,y,1776));
 225 ;;         $$f_{\left(1066,1776\right)}(x,y)$$
 226
 227 ;; Finally, setting tex_uses_named subscripts and tex_uses_prime_for_derivatives
 228 ;; to false, derivatives are indicated with parenthesis surrounded
 229 ;; subscripts.  There is one subscript for each function argument; when
 230 ;; the derivative order is zero, the subscript is zero.
 231
 232 ;;  (c11) tex_uses_prime_for_derivatives : false;
 233 ;;  (c12) tex_uses_named_subscripts_for_derivatives : false;
 234 ;;  (c13) tex(diff(f(a,b),a,2,b,1));
 235 ;;            $$f_{\left(2,1\right)}(a,b)$$
 236 ;;  (c14) tex(diff(f(a,b),a,0,b,1));
 237 ;;            $$f_{\left(0,1\right)}(a,b)$$
 238 ;;  (c15) tex(diff(f(x,y),x,0,y,1));
 239 ;;            $$f_{\left(0,1\right)}(x,y)$$
 240
 241 (defun tex-pderivop (x l r)
 242   ;(print `(lop = ,lop rop = ,rop x = ,x r = ,r l = ,l))
 243   (cond ((and $tex_uses_prime_for_derivatives (eql 3 (length x)))
 244          (let* ((n (car (last x)))
 245                 (p))
 246
 247            (cond ((<= n $tex_prime_limit)
 248                   (setq p (make-list n :initial-element "\\prime")))
 249                  (t
 250                   (setq p (list "(" n ")"))))
 251
 252            ;; We need to avoid double tex superscripts; when rop is mexpt, use parens.
 253
 254            (cond ((eq rop 'mexpt)
 255                   (append l (list "\\left(") (tex (cadr x) nil nil lop rop)
 256                           (list "^{") p (list "}") (list "\\right)") r))
 257                  (t
 258                   (append l  (tex (cadr x) nil nil lop rop)
 259                           (list "^{") p (list "}")  r)))))
 260
 261         ((and $tex_uses_named_subscripts_for_derivatives
 262               (< (apply #'+ (cddr x)) $tex_prime_limit))
 263          (let ((n (cddr x))
 264                (v (mapcar #'stripdollar (cdr $tex_diff_var_names)))
 265                (p))
 266
 267            (cond ((> (length n) (length v))
 268                   (merror "Not enough elements in tex_diff_var_names to tex the expression")))
 269            (dotimes (i (length n))
 270              (setq p (append p (make-list (nth i n) :initial-element (nth i v)))))
 271            (append l (tex (cadr x) nil nil lop rop) (list "_{") p (list "}") r)))
 272
 273         (t
 274          (append l (tex (cadr x) nil nil lop rop)  (list "_{")
 275                  (tex-matchfix (cons '(mprogn)
 276                                      (cddr x)) nil nil) (list "}") r))))
 277
 278 ;; Convert from positional derivatives to non-positional derivatives.
 279
 280 (defun $convert_to_diff (e)
 281   (setq e ($totaldisrep e))
 282   (cond (($mapatom e) e)
 283         ((and (consp e) (eq (caar e) 'mqapply) (eq (caaadr e) '$pderivop))
 284          (let ((f (second (second e)))
 285                (n (rest (rest (second e))))
 286                (v (rest (rest e)))
 287                (g nil)
 288                ($use_pdiff nil)
 289                (vk) (gk)
 290                (at-list nil))
 291            ;; Work around bug in Clozure CL: https://github.com/Clozure/ccl/issues/109
 292            ;; If ever it's fixed, this conditionalization can be removed.
 293            #+ccl (mapcar #'(lambda (vs) (declare (ignore vs)) (push (gensym) g)) v)
 294            #-ccl (dolist (vs v) (declare (ignore vs)) (push (gensym) g))
 295            (setq e (mapply f g nil))
 296            (setq e (apply '$diff (cons e (mapcan #'list g n))))
 297            (while v
 298              (setq vk (pop v))
 299              (setq gk (pop g))
 300              (if (symbolp vk) (setq e ($substitute vk gk e)) (push (take '(mequal) gk vk) at-list)))
 301            (if at-list ($at e (simplify (cons '(mlist) at-list))) e)))
 302         (t (mapply (mop e) (mapcar #'$convert_to_diff (margs e)) nil))))
 303