tests/wester_problems/test_vector_analysis.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/VectorAnalysis/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 /* ---------- Vector Analysis ---------- */
  19 /* Vector norm => sqrt(15) */
  20 mat_norm([1 + %i, -2, 3*%i]);
  21 load(vect)$
  22 /* Cross product: (2, 2, -3) x (1, 3, 1) => (11, -5, 4) */
  23 vect_express([2, 2, -3] ~ [1, 3, 1]);
  24 declare([a, b, c, d, f, g], nonscalar)$
  25 /* (a x b) . (c x d) => (a . c) (b . d) - (a . d) (b . c) */
  26 (a ~ b) . (c ~ d);
  27 vectorsimp(%), expandall;
  28 /* => (2 y z^3 - 2 x^2 y^2 z,   x y,   2 x y^2 z^2 - x z) */
  29 vect_express(curl [x*y*z, x^2*y^2*z^2, y^2*z^3]);
  30 ev(%, diff);
  31 /* DEL . (f x g) => g . (DEL x f) - f . (DEL x g) */
  32 div (f ~ g);
  33 vectorsimp(%), expandall;
  34 remove([a, b, c, d, f, g], nonscalar)$
  35 /* Express DEL . a in spherical coordinates (r, theta, phi) for
  36    a = (a_r(r, theta, phi), a_theta(r, theta, phi), a_phi(r, theta, phi)).
  37    Here, phi is in the x-y plane and theta is the angle with the z-axis.
  38    => 1/r^2 d/dr[r^2 a_r] + 1/[r sin(theta)] d/dtheta[sin(theta) a_theta]
  39       + 1/[r sin(theta)] da_phi/dphi
  40    => da_r/dr + (2 a_r)/r + 1/r da_theta/dtheta + a_theta/[r tan(theta)]
  41       + 1/[r sin(theta)] da_phi/dphi
  42    See Keith R. Symon, _Mechanics_, Third Edition, Addison-Wesley Publishing
  43    Company, 1971, p. 103. */
  44 vect_coordsys(spherical)$
  45 vect_coords;
  46 depends([a_r, a_theta, a_phi], [r, theta, phi]);
  47 vect_express(div [a_r, a_phi, a_theta]);
  48 expand(ev(%, 'diff));
  49 remove([a_r, a_theta, a_phi], dependency)$
  50 /* Express dR/dt in spherical coordinates (r, theta, phi) where R is the
  51    position vector r*Rhat(theta, phi) with Rhat being the unit vector in the
  52    direction of R => (dr/dt, r dtheta/dt, r sin(theta) dphi/dt)
  53    [Symon, p. 98] */
  54 depends([r, theta, phi], t, rhat, [theta, phi])$
  55 vect_express(diff([r*rhat, 0, 0], t));
  56 remove([r, theta, phi, rhat], dependency, ".", commutative)$
  57 vect_coordsys(cartesian3d)$
  58 /* Scalar potential => x^2 y + y + 2 z^3 */
  59 scalar_potential([2*x*y, x^2 + 1, 6*z^2]);
  60 /* Vector potential => (x y z, x^2 y^2 z^2, y^2 z^3) is one possible solution.
  61    See Harry F. Davis and Arthur David Snider, _Introduction to Vector
  62    Analysis_, Third Edition, Allyn and Bacon, Inc., 1975, p. 97. */
  63 vector_potential([2*y*z^3 - 2*x^2*y^2*z, x*y, 2*x*y^2*z^2 - x*z]);
  64 vect_express(curl(%))$
  65 ratsimp(ev(%, diff));
  66 /* Orthogonalize the following vectors (Gram-Schmidt).  See Lee W. Johnson and
  67    R. Dean Riess, _Introduction to Linear Algebra_, Addison-Wesley Publishing
  68    Company, 1981, p. 104 => [[0 1 2 1], [0 -1 1 -1], [2 1 0 -1]]^T */
  69 [[0, 1, 2, 1], [0, 1, 3, 1], [1, 1, 1, 0], [1, 3, 6, 2]];
  70 gramschmidt(%);
  71 /* ---------- Quit ---------- */
  72 quit();