tests/wester_problems/test_programming_dif.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/Programming/dif.macsyma
   3  * circa 2006-10-23.
   4  *
   5  * Released under the terms of the GNU General Public License, version 2,
   6  * per message dated 2007-06-03 from Michael Wester to Robert Dodier
   7  * (contained in the file wester-gpl-permission-message.txt).
   8  *
   9  * See: "A Critique of the Mathematical Abilities of CA Systems"
  10  * by Michael Wester, pp 25--60 in
  11  * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
  12  * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
  13  */
  14 /* ----------[ M a c s y m a ]---------- */
  15 /* ---------- Initialization ---------- */
  16 showtime: all$
  17 prederror: false$
  18 /* ---------- Programming and Miscellaneous ---------- */
  19 /* Define a simple differentiation operator.  This is most easily done with
  20    pattern matching.  Here is an expression that we ultimately hope to be able
  21    to differentiate. */
  22 expr: a + b*u + c*u^2 + d*(e*u + f)^3 + g*u*exp(h*u);
  23 dif(e, x):= block([oper, n],
  24    if not atom(e) then
  25       oper: part(e, 0),
  26    n: length(e),
  27 /* Start by making the derivative of a sum to be the sum of the derivatives. */
  28    if oper = "+" then
  29       map(lambda([z], dif(z, x)), e)
  30 /* Add the product rule. */
  31    else if oper = "*" then
  32       dif(part(e, 1), x) * part(e, 2..n) + part(e, 1) * dif(part(e, 2..n), x)
  33 /* Now, make the derivative of a constant (with respect to x) zero. */
  34    else if freeof(x, e) then
  35        0
  36 /* Define the derivative of x with respect to x to be one. */
  37    else if e = x then
  38        1
  39 /* Enter the generalized power rule. */
  40    else if oper = "^" and not freeof(x, part(e, 1)) then
  41        part(e, 2)*part(e, 1)^(part(e, 2) - 1)*dif(part(e, 1), x)
  42 /* To get that last term, add in the exponential rule. */
  43    else if oper = "^" and part(e, 1) = %e then
  44       e*dif(part(e, 2), x)
  45    );
  46 /* Now, try it out! => b + 2 c u + 3 d e (e u + f)^2 + g exp(h u) (1 + h u) */
  47 dif(expr, u);