share/contrib/vector3d/vector3d.dem

   1 /* File name: vector3d.dem
   2  * Demonstration and test of the vector3d package
   3  * Walter Eastes: June, 2004
   4  */
   5 load("vector3d.mac");
   6
   7 "Let's define some vectors to use"$
   8 a : [ax,ay,az]; b : [bx,by,bz]; c : [cx,cy,cz];
   9
  10 "These are vector and scalar functions of the coordinates"$
  11 u : [ux,uy,uz]; v : [vx,vy,vz];
  12 depends([ux,uy,uz,vx,vy,vz], [x,y,z]); depends([f], [x,y,z]);
  13
  14 "Now let's prove some elementary identities"$
  15 is(cross(a,b) = -cross(b,a));
  16 is(cross(a,h*b) = cross(h*a,b));
  17 is(cross(a,a) = [0,0,0]);
  18 is( expand(cross(a, cross(b,c))) = expand((a . c)*b - (a . b)*c) );
  19
  20 "Now some identities for the differential operators"$
  21 is( expand(div(cross(u,v))) = expand(curl(u) . v - u . curl(v)) );
  22 is( curl(curl(u)) = grad(div(u)) - laplacian(u) );
  23 is( div(grad(f)) = laplacian(f) );
  24
  25 "Now let's change to spherical polar coordinates"$
  26 coordsys(spherical);
  27
  28 "Define vector functions of these coordinates"$
  29 u : [ur,uth,uph]; v : [vr,vth,vph];
  30 depends([ur,uth,uph,vr,vth,vph], [r,th,ph]); depends([f], [r,th,ph]);
  31
  32 "Now see if the differential operator identities still hold"$
  33 is( expand(div(cross(u,v))) = expand(curl(u) . v - u . curl(v)) );
  34 is( expand(curl(curl(u))) = expand(grad(div(u)) - laplacian(u)) );
  35 is(equal( div(grad(f)), laplacian(f) ));
  36
  37 "Finally prove the identity in the original demo"$
  38 "Need these vectors first"$
  39 a : [ar,ath,aph]; b : [br,bth,bph];
  40 is( expand( cross(a,b) . cross(u,v) + cross(b,u) . cross(a,v) +
  41    cross(u,a) . cross(b,v) ) = 0);