Maxima doesn't easily work with derivatives similar to
diff(f(x^2),x) or diff(f(x,x^2),x); Macsyma had an ndiff
package that supported positional derivatives; it made working with
derivatives like these easy. I used the ndiff package for
perturbational calculations and for changing variables in partial
differential equations; with Maxima, these calculations
are more lengthy and error prone.
The other day, I tried to add an ndiff-like package to Maxima.
I pasted a few lines into sdiffgrad and hacked dimension
(a displa related function). Here is an example of what I was
able to do:
(C1) load("pdiff.lisp")$
Loading pdiff.lisp
Finished loading pdiff.lisp
/* Declare that derivatives of f are to be given positionally. */
(C2) qput(f,true,pderiv)$
/* Derivatives of f are indicated with subscripts surrounded with
parenthesis; f_(i) is the i-th derivative of f with respect to the
first argument of f; f_(i,j) is the j-th derivative with respect to
the second argument of f_(i); and etc. */
(C3) diff(f(x),x);
(D3) f (x)
(1)
(C4) diff(f(x^2),x);
2
(D4) 2 x f (x )
(1)
(C5) diff(%,x);
2 2 2
(D5) 4 x f (x ) + 2 f (x )
(2) (1)
(C6) diff(f(x,y),x,1,y,1);
(D6) f (x, y)
(1, 1)
(C7) diff(f(x,x^2),x);
2 2
(D7) f (x, x ) + 2 x f (x, x )
(1, 0) (0, 1)
(C8) diff(f(t),t);
(D8) f (t)
(1)
/* Derivatives work okay; but there are many problems related to evaluation. */
(C9) ev(%);
Improper name or value in functional position:
f
(1)
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
/* I think this error comes from BADFUNCHK that is called by
mapply1; the real problem could be much further up. What do you
all think? */
/* Internally, my positional derivatives look like: */
(C10) diff(f(x),x);
(D10) f (x)
(1)
(C11) ?print(%);
((MQAPPLY SIMP) ((%PDERIVOP SIMP) |$f| 1) |$x|)
What would it take to integrate something like my positional derivative
code into Maxima? (Fix the ev problem and other bugs.)
Regards,
Barton, TOKS (The Old Kind of Scientist)
Here's my code.
--------------------------------------------------
(defun pderivop (f x n)
(append (list '(mqapply simp) (append (list '(%pderivop simp) f) n)) x))
;; Return the list (k1,k2,k3,...,kn), where ki = 1 and all other k's are zero.
(defun i-list (i n)
(let ((s nil))
(dotimes (k n (reverse s))
(cond ((eq k i)
(setq s (cons 1 s)))
(t
(setq s (cons 0 s)))))))
;; Increment the ith element of the list e by one.
(defun incf-ith (i e)
(let ((k (nth i e))
(q (copy-list e)))
(setf (nth i q) (add 1 k))
q))
;; Positional derivatives with symbolic orders do not "work." Maybe a
;; function further up needs to be hacked.
(defun sdiffgrad (e x)
(let ((fun (caar e)) grad args)
(cond ((and (eq fun 'mqapply) (oldget (caaadr e) 'grad))
(sdiffgrad (cons (cons (caaadr e) nil) (append (cdadr e) (cddr e))) x))
;; stuff I added-----------------------------------
((and (eq fun 'mqapply) (eq (caaadr e) '%PDERIVOP))
(setq args (cddr e))
(setq args (mapcar #'(lambda (s) ($diff s x)) args))
(setq fun (cadadr e))
(let ((de 0)
(n (length args))
(d-order (cddadr e)))
(dotimes (i n de)
(setq de (add de (mul (nth i args)
(pderivop fun (cddr e)
(incf-ith i d-order))))))))
(($get (caar e) '|$pderiv|)
(setq args (cdr e))
(setq args (mapcar #'(lambda (s) ($diff s x)) args))
(let ((de 0)
(n (length args)))
(dotimes (i n de)
(setq de (add de (mul (nth i args)
(pderivop (caar e) (cdr e) (i-list i n))))))))
;; end of stuff I added-----------------------------
((or (eq fun 'mqapply) (null (setq grad (oldget fun 'grad))))
(if (not (depends e x)) 0 (diff%deriv (list e x 1))))
((not (= (length (cdr e)) (length (car grad))))
(merror "wrong number of arguments for ~:m" fun))
(t (setq args (sdiffmap (cdr e) x))
(addn (mapcar
#'mul2
(cdr (substitutel
(cdr e) (car grad)
(do ((l1 (cdr grad) (cdr l1))
(args args (cdr args)) (l2))
((null l1) (cons '(mlist) (nreverse l2)))
(setq l2 (cons (cond ((equal (car args) 0) 0)
(t (car l1)))
l2)))))
args)
t)))))
;; Hacked code for displaying positional derivatives. Powers of positoinal
;; derivatives display incorrectly.
(defun dimension (form result lop rop w right)
(let ((level (f1+ level))
(break (if (and w break) (f+ w break))))
(setq form (nformat-check form))
(cond ((atom form)
(dimension-atom form result))
((and (atom (car form)) (setq form (cons '(mprogn) form)) nil))
((or (<= (lbp (caar form)) (rbp lop)) (> (lbp rop) (rbp (caar form))))
(dimension-paren form result))
((memq 'array (car form))
(dimension-array form result))
;; Hokie hack for display of positional derivatives.
((memq '%pderivop (car form))
(setq form (cdr form))
(setq form (list (list (car form) 'simp 'array)
(cons '("") (cdr form))))
(dimension-array form result))
((safe-get (caar form) 'dimension)
(funcall (get (caar form) 'dimension) form result))
(t
(dimension-function form result)))))
;; Taylor doesn't work with my positional derivatives, but fake_taylor
;; does.
; fake_taylor(f,v,pt,ord) := block([ps, i, s],
; ps : 0,
; s : 1,
; ord : 1 + ord,
; for i : 1 thru ord do (
; ps : ps + at(f, v = pt) * s,
; s : s * (v - pt) / i,
; f : diff(f,v)),
; ps);
--------------------------------------------------------