share/vector/rtest_vect.mac

   1 /*******************************************************************************
   2  *
   3  *          Examples for the package vect
   4  *
   5  ******************************************************************************/
   6
   7 kill(all);
   8 done;
   9
  10 (load("vect.mac"),done);
  11 done;
  12
  13 declare([a,b,c,d],nonscalar);
  14 done;
  15
  16 /* ---------- Check some rules for vectors  ----------------------------------*/
  17
  18 /* This does not simplify to zero, because we have removed the property
  19  * commutative for the operator "." */
  20 a . b - b . a;
  21 0;
  22
  23 a ~ b + b ~ a;
  24 0;
  25
  26 vectorsimp(a.(b+c)-(a.b+a.c)),expandall:true;
  27 0;
  28
  29 vectorsimp(a~(b+c)-(a~b+a~c)),expandall:true;
  30 0;
  31
  32 vectorsimp(a~(b~c)-(b*(a.c)-c*(a.b))),expandall:true;
  33 0;
  34
  35 /* Lagrange's identity does not simplify as expected */
  36 vectorsimp((a~b).(c~d)-((a.c)*(b.d)-(b.c)*(a.d))),expandall:true;
  37 0;
  38
  39 /* Does not simplify to the expected result. */
  40 vectorsimp((3*a-2*b)~(2*a+b)),expandall:true;
  41 7*(a~b);
  42
  43 /* ---------- Check some rules for the gradient ------------------------------*/
  44
  45 declare([u,v],scalar, c,constant);
  46 done;
  47 depends([u,v],[x,y,z]);
  48 [u(x,y,z), v(x,y,z)];
  49
  50 grad c;
  51 0;
  52
  53 grad (c*u);
  54 c*grad u;
  55
  56 vectorsimp(grad (u+v)),expandall:true;
  57 grad u + grad v;
  58
  59 /* ---------- Check some rules for the divergence ----------------------------*/
  60
  61 div c;
  62 0;
  63
  64 div (c*u);
  65 c*div u;
  66
  67 vectorsimp(div (u+v)),expandall:true;
  68 div u + div v;
  69
  70 /* The error is a . grad u -> a grad u */
  71 vectorsimp(div (u*a)) - (u*div a + (grad u) . a),expandall:true;
  72 0;
  73
  74 /* There is no rule do expand the following */
  75 vectorsimp(div (a ~ b)),expandall:true;
  76 b . curl a - a . curl b;
  77
  78 /* ---------- Check some rules for the curl ----------------------------------*/
  79
  80 curl c;
  81 0;
  82
  83 curl (c*u);
  84 c* curl u;
  85
  86 vectorsimp(curl (a+b)),expandall:true;
  87 curl a + curl b;
  88
  89 vectorsimp(curl (a*u)),expandall:true;
  90 u*curl a + grad u ~ a;
  91
  92 vectorsimp(curl (a ~ b)),expandall:true;
  93 b . grad a - a . grad b + a * div b - b * div a;
  94
  95 /* ---------- Combinations of the operators ----------------------------------*/
  96
  97 div curl a;
  98 0;
  99 curl grad u;
 100 [0, 0, 0];
 101 vectorsimp(curl curl a),expandall:true;
 102 grad div a - laplacian a;
 103
 104 /* ---------- Examples from bug reports --------------------------------------*/
 105
 106 /* Bud ID: 838301 - vect negate cross product simplification
 107  * a and b are declare to be nonscalar, c to be constant
 108  */
 109
 110 c*(a ~ b);
 111 c * (a ~ b);
 112 c*(b ~ a);
 113 -c * (a ~ b);
 114
 115 /* Bug ID: 1212598 - bug in the VECT.MAK - VECTORSIMP cross product
 116  * a,b,c,d are declared to be nonscalar
 117  */
 118
 119 vectorsimp( (b+2*a) ~ (c+d) ),expandall:true;
 120 b ~ d+b ~ c+2*a ~ d+2*a ~ c;
 121
 122
 123 /* Bug ID: 2806446 - ev_diff in vect.mac */
 124
 125 potential([y*z,x*z,x*y],[x,y,z]);
 126 x*y*z;
 127
 128 /* Bug ID: 2011228 -vect redefines "." as commutative, was:Matrix multiplication
 129  * Two examples which does not work when "." is declared commutative
 130  */
 131
 132 invert(matrix([-2, -2, 1], [3, 1, 1], [3, 3, 1]));
 133 matrix([-1/5,1/2,-3/10],[0,-1/2,1/2],[3/5,0,2/5]);
 134
 135 (a : matrix([1,2],[6,7]), b : matrix([1,2],[1,6]), done);
 136 done;
 137 a.b-b.a;
 138 matrix([-10,-2],[-24,10]);
 139
 140 /*
 141
 142 Some syntax errors can be detected by comparing the declared part of
 143 speech to an actual expression. (see documentation for `infix'.). The
 144 tests are not automatic because a syntax error is not captured by
 145 `errcatch'.
 146
 147 errcatch(if a~b then 1 else 0);
 148 []$
 149
 150 errcatch(if grad a then 1 else 0);
 151 []$
 152
 153 */
 154
 155 /* SF bug #3337: "Wrong scalefactor for cartesian2d" */
 156
 157 (scalefactors (cartesian2d), sf);
 158 [1, 1];
 159
 160 (scalefactors (cartesian3d), sf);
 161 [1, 1, 1];
 162
 163 /* Additional tests for #3337. These are weak in the sense that
 164  * they aren't testing for specific values, but we would have to
 165  * work out what is the correct result in each case.
 166  */
 167
 168 (scalefactors (polar), is (length(sf) = length(first(polar))));
 169 true;
 170
 171 (scalefactors (polarcylindrical), is (length(sf) = length(first(polarcylindrical))));
 172 true;
 173
 174 (scalefactors (spherical), is (length(sf) = length(first(spherical))));
 175 true;
 176
 177 (scalefactors (elliptic), is (length(sf) = length(first(elliptic))));
 178 true;
 179
 180 (scalefactors (ellipticcylindrical), is (length(sf) = length(first(ellipticcylindrical))));
 181 true;
 182
 183 (scalefactors (confocalelliptic), is (length(sf) = length(first(confocalelliptic))));
 184 true;
 185
 186 (scalefactors (prolatespheroidalsqrt), is (length(sf) = length(first(prolatespheroidalsqrt))));
 187 true;
 188
 189 (scalefactors (oblatespheroidalsqrt), is (length(sf) = length(first(oblatespheroidalsqrt))));
 190 true;
 191
 192 (scalefactors (parabolic), is (length(sf) = length(first(parabolic))));
 193 true;
 194
 195 (scalefactors (paraboliccylindrical), is (length(sf) = length(first(paraboliccylindrical))));
 196 true;
 197
 198 (scalefactors (paraboloidal), is (length(sf) = length(first(paraboloidal))));
 199 true;
 200
 201 (scalefactors (prolatespheroidal), is (length(sf) = length(first(prolatespheroidal))));
 202 true;
 203
 204 (scalefactors (oblatespheroidal), is (length(sf) = length(first(oblatespheroidal))));
 205 true;
 206
 207 (scalefactors (bipolar), is (length(sf) = length(first(bipolar))));
 208 true;
 209
 210 (scalefactors (bipolarcylindrical), is (length(sf) = length(first(bipolarcylindrical))));
 211 true;
 212
 213 (scalefactors (toroidal), is (length(sf) = length(first(toroidal))));
 214 true;
 215
 216 (scalefactors (conical), is (length(sf) = length(first(conical))));
 217 true;
 218
 219 (scalefactors (confocalellipsoidal), is (length(sf) = length(first(confocalellipsoidal))));
 220 true;
 221
 222 (reset(tr_bound_function_applyp,vect_cross,dotexptsimp,dotassoc,dotscrules),0);
 223 0$
 224
 225 /* SF bug #3287: "cross product scalar zero versus vector zero result" */
 226
 227 kill (f, g);
 228 done;
 229
 230 0 ~ f(a, b);
 231 [0, 0, 0];
 232
 233 g(c, d) ~ 0;
 234 [0, 0, 0];
 235
 236 f(a, g(b, c)) ~ f(a, g(b, c));
 237 [0, 0, 0];
 238
 239 curl (grad (f (a, b, c)));
 240 [0, 0, 0];
 241
 242 /* examples from the bug report */
 243
 244 express([0,1,1]~[0,1,2]);
 245 [1, 0, 0];
 246
 247 express([0,1,1]~[0,1,1]);
 248 [0, 0, 0];
 249
 250 (e:ident(3),
 251  a_:[a_1,a_2,a_3],
 252  b_:[b_1,b_2,b_3],
 253  /* Levi-Civita Symbol */
 254  ε(i,j,k):=e[i].express(e[j]~e[k]),
 255 0);
 256 0;
 257
 258 ε(1,2,3);
 259 1;
 260
 261 ε(1,3,2);
 262 -1;
 263
 264 ε(1,2,2);
 265 0;
 266
 267 express(ε(1,2,3)*a_[2]*b_[3]);
 268 a_2*b_3;
 269
 270 express(ε(1,3,2)*a_[2]*b_[3]);
 271 -a_2*b_3;
 272
 273 express(ε(1,2,2)*a_[2]*b_[2]);
 274 0;
 275
 276 makelist(sum(sum(ε(i,j,k)*a_[j]*b_[k],j,1,3),k,1,3),i,3);
 277 [a_2*b_3 - a_3*b_2, a_3*b_1 - a_1*b_3, a_1*b_2 - a_2*b_1];
 278
 279 /* bug reported to mailing list 2022-01-04: "Bad vect package?" */
 280
 281 /* example from email -- grad */
 282
 283 (scalefactors(polarcylindrical),
 284  kill (u, E, rho, I, D),
 285  depends(u,r),
 286  E:[-rho*I,0,0],
 287  0);
 288 0;
 289
 290 foo: D*(grad(u)+E*u);
 291 D*([- I*rho*u, 0, 0] + grad(u));
 292
 293 ev (express (foo), 'diff);
 294 [D*('diff(u, r) - I*rho*u), 0, 0];
 295
 296 /* other examples -- cross product, curl */
 297
 298 (kill (a, b, x, y, z),
 299  (a ~ b) + [x, y, z]);
 300 (a ~ b) + [x, y, z];
 301
 302 curl(a) + [x, y, z];
 303 curl(a) + [x, y, z];