share/contrib/impdiff.mac

   1 /*
   2 This subroutine implicit derivatives of multivariable expressions
   3 Copyright (C) 1999  Dan Stanger
   4
   5 This library is free software; you can redistribute it and/or modify it
   6 under the terms of the GNU Library General Public License as published
   7 by the Free Software Foundation; either version 2 of the License, or (at
   8 your option) any later version.
   9
  10 This library is distributed in the hope that it will be useful, but
  11 WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 Library General Public License for more details.
  14
  15 You should have received a copy of the GNU Library General Public
  16 License along with this library; if not, write to the Free Software
  17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  18
  19 Dan Stanger dan.stanger@internut.com
  20 please contact me for updates to this code or
  21 look on my web page (currently www.diac.com/~dxs)
  22 this code was tested using Macsyma 2.4 from Macsyma, Inc.
  23 */
  24 /* f is an array, the indexes are the derivative degree
  25      in the indvarlist order.
  26    indvarlist independent variable list
  27    depvar the dependent variable
  28    orderlist the order desired
  29 */
  30
  31 /* If makeOrders is not already defined, then load it.
  32  * Following test looks for a Maxima function named makeOrders.
  33  * Probably should also look for a Lisp function,
  34  * in case makeOrders is reimplemented as a Lisp function.
  35  */
  36 if ?mget ('makeOrders, '?mexpr) = false then load (makeOrders);
  37
  38 /* Display an array, either a hashed array or a declared array,
  39  * as a list of equations a[<subscripts>] = <value>.
  40  * Seems generally useful -- might want to move this into core at some point.
  41  */
  42 display_array ('a) := block ([stuff, ainfo],
  43   stuff : errcatch (apply (arrayinfo, [a])),
  44   if stuff # []
  45     then
  46      (ainfo : first (stuff),
  47       if first (ainfo) = declared
  48         then display_array_declared (a, ainfo[3])
  49         else display_array_hashed (a, rest (rest (ainfo)))));
  50
  51 display_array_declared (a, dimensions) :=
  52  (map (lambda ([n], apply (set, makelist (i, i, 1, n))), dimensions),
  53   apply (cartesian_product, %%),
  54   makelist (arraymake (a, i), i, args (%%)),
  55   apply (display, %%));
  56
  57 display_array_hashed (a, indices) :=
  58  (makelist (arraymake (a, i), i, indices),
  59   apply (display, %%));
  60
  61 implicit_derivative(f,indvarlist,orderlist,depvar):=
  62    block([
  63         l:makeOrders(indvarlist,orderlist),
  64         orders,orderslength,diffargs],local(orders,diffargs),
  65         orderslength:length(l)-1,
  66         array([orders,diffargs],orderslength),
  67         fillarray(orders,l),
  68         depends(depvar,indvarlist),
  69         for i:1 thru orderslength do block([d,s],
  70            diffargs[i]:apply('append,maplist(lambda([x,y],[x,y]),
  71                               indvarlist,orders[i])),
  72            d:apply('diff, cons(arrayapply(f,orders[0]), diffargs[i])),
  73            for j:(i-1) step -1 thru 1 do block(
  74                 d:subst(
  75                    arrayapply(f,orders[j]),
  76                    apply('diff,cons(depvar,diffargs[j])), d)),
  77            s:solve(d,apply('diff,cons(depvar,diffargs[i]))),
  78            arraysetapply(f,orders[i],rhs(first(s)))),
  79        f)$
  80
  81 /* the following example will fill the array f with derivatives
  82
  83     load (impdiff);
  84     f [0, 0] : x^2 + y^3 - z^4 = 0;
  85     implicit_derivative (f, [x, y], [2, 3], z);
  86     display_array (f);
  87
  88  */
  89