tests/wester_problems/test_programming.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/Programming/problems.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 /* How easy is it to substitute x for a + b in the following expression?
  20    => (x + c)^2 + (d - x)^2 */
  21 expr: (a + b + c)^2 + (d - a - b)^2;
  22 ratsubst(x, a + b, expr);
  23 factorsum(%);
  24 subst(b = x - a, expr);
  25 remvalue(expr)$
  26 /* How easy is it to substitute r for sqrt(x^2 + y^2) in the following
  27    expression? => x/r */
  28 x/sqrt(x^2 + y^2);
  29 subst(sqrt(x^2 + y^2) = r, %);
  30 /* Change variables so that the following transcendental expression is
  31    converted into a rational expression   [Vernor Vinge]
  32    => (r - 1)^4 (u^4 - r u^3 - r^3 u + r u + r^4)/[u^4 (2 r - 1)^2] */
  33 q: (1/r^4 + 1/(r^2 - 2*r*cos(t) + 1)^2
  34           - 2*(r - cos(t))/(r^2 * (r^2 - 2*r*cos(t) + 1)^(3/2))) /
  35    (1/r^4 + 1/(r - 1)^4 - 2*(r - 1)/(r^2 * (r^2 - 2*r + 1)^(3/2)));
  36 assume(r > 1, u >= r - 1, u > 0);
  37 subst([r^2 - 2*r*cos(t) + 1 = u^2, cos(t) = (r^2 - u^2 + 1)/(2*r)], q);
  38 factor(scanmap(factor, %));
  39 forget(r > 1, u >= r - 1, u > 0)$
  40 /* Establish a rule to symmetrize a differential operator:   [Stanly Steinberg]
  41    f g'' + f' g' -> (f g')' */
  42 defrule(symmetrize, f(x)*diff(g(x), x, 2) + diff(f(x), x)*diff(g(x), x),
  43                     'diff(f(x)*'diff(g(x), x), x));
  44 q: f(x)*diff(g(x), x, 2) + diff(f(x), x)*diff(g(x), x);
  45 apply1(q, symmetrize);
  46 /* => 2 (f g')' + f g */
  47 apply1(2*q + f(x)*g(x), symmetrize);
  48 apply1(expand(2*q + f(x)*g(x)), symmetrize);
  49 matchdeclare(x, true)$
  50 defrule(symmetrize, f(x)*diff(g(x), x, 2) + diff(f(x), x)*diff(g(x), x),
  51                     'diff(f(x)*'diff(g(x), x), x));
  52 q: f(y)*diff(g(y), y, 2) + diff(f(y), y)*diff(g(y), y);
  53 apply1(q, symmetrize);
  54 apply1(2*q + f(y)*g(y), symmetrize);
  55 matchdeclare(f, true, g, true)$
  56 defrule(symmetrize, f(x)*diff(g(x), x, 2) + diff(f(x), x)*diff(g(x), x),
  57                     'diff(f(x)*'diff(g(x), x), x));
  58 q: ff(x)*diff(gg(x), x, 2) + diff(ff(x), x)*diff(gg(x), x);
  59 errcatch(apply1(q, symmetrize));
  60 remvalue(q)$
  61 /* Infinite lists: [1 2 3 4 5 ...] * [1 3 5 7 9 ...]
  62    => [1 6 15 28 45 66 91 ...] */
  63 l1: [1, 2, 3, 4, 5]$
  64 l2: [1, 3, 5, 7, 9]$
  65 l1 * l2;
  66 remvalue(l1, l2)$
  67 /* Write a simple program to compute Legendre polynomials */
  68 p[n](x):= ratsimp(1/(2^n*n!) * diff((x^2 - 1)^n, x, n));
  69 /* p[0](x) = 1,   p[1](x) = x,   p[2](x) = (3 x^2 - 1)/2,
  70    p[3](x) = (5 x^3 - 3 x)/2,   p[4](x) = (35 x^4 - 30 x^2 + 3)/8 */
  71 for i:0 thru 4 do display(p[i](x))$
  72 /* p[4](1) = 1 */
  73 p[4](1);
  74 /* Now, perform the same computation using a recursive definition */
  75 pp[0](x):= 1$
  76 pp[1](x):= x$
  77 pp[n](x):= ratsimp(((2*n - 1)*x*pp[n - 1](x) - (n - 1)*pp[n - 2](x))/n);
  78 for i:0 thru 4 do disp('pp[i](x) = pp[i](x))$
  79 pp[4](1);
  80 remarray(p, pp)$
  81 /* Iterative computation of Fibonacci numbers */
  82 myfib(n):= block([i, j, k, f],
  83            if n < 0 then
  84               error('undefined)
  85            else if n < 2 then
  86               n
  87            else
  88              (j: 0,   k: 1,
  89               for i:2 thru n do
  90                 (f: j + k,   j: k,   k: f),
  91               f))$
  92 /* Convert the function into FORTRAN syntax */
  93 errcatch(apply('gentran, [%]));
  94 /* Create a list of the first 11 values of the function. */
  95 makelist(myfib(i), i, 0, 10);
  96 remfunction(myfib)$
  97 /* Define the function p(x) = x^2 - 4 x + 7 such that p(lambda) = 0 for
  98    lambda = 2 +- i sqrt(3) and p(A) = [[0 0], [0 0]] for A = [[1 -2], [2 3]]
  99    (the lambda are the eigenvalues and p(x) is the characteristic polynomial of
 100    A)   [Johnson and Reiss, p. 184] */
 101 p(x):= x^2 - 4*x + 7;
 102 ratsimp(p(2 + %i*sqrt(3)));
 103 p(matrix([1, -2], [2, 3]));
 104 remfunction(p)$
 105 /* Define a function to be the result of a calculation */
 106 -log(x^2 - 2^(1/3)*x + 2^(2/3))/(6 * 2^(2/3))
 107    + atan((2*x - 2^(1/3))/(2^(1/3) * sqrt(3))) / (2^(2/3) * sqrt(3))
 108    + log(x + 2^(1/3))/(3 * 2^(2/3));
 109 define(f(x), %);
 110 expr: f(y);
 111 /* Display the top-level structure of a nasty expression, hiding the
 112    lower-level details. */
 113 reveal(%, 1);
 114 reveal(expr, 2);
 115 remvalue(expr)$
 116 remfunction(f)$
 117 /* Convert the following expression into TeX or LaTeX */
 118 declare(x, complex)$
 119 y = sqrt((exp(x^2) + exp(-x^2))/(sqrt(3)*x - sqrt(2)));
 120 tex(%);
 121 remove(x, complex)$