share/vector/vect.mac

   1 /*-*-MACSYMA-*-*/
   2 /* Or use the BATCHLOAD command to load this with TTYOFF:TRUE */
   3 /*  NOTE:  THE CURRENT VERSION OF VECT IS THE ONE DUE TO STOUTEMYER.
   4 IT WILL BE REPLACED SOON BY AN EXTENDED VERSION WHICH HANDLES BOTH
   5 VECTORS AND DYADICS.
   6         MICHAEL C. WIRTH (MCW)
   7         12/18/78
   8 Style changes made in order to TRANSLATE, 3/1/81 George Carrette (GJC)
   9 */
  10
  11 herald_package(vect)$
  12 put('vect,true,'version);
  13
  14 eval_when([translate,batch,demo,load,loadfile],
  15   /* trgsmp is in the standard Maxima image but check anyway */
  16   if get('sin,'complement_function)=false then load("trgsmp"),
  17   if not ?boundp('coords) then load("vect_transform"),
  18   matchdeclare([etrue,ttrue,vtrue], true, lessp, before, scalarm, vscalarp),
  19   tr_bound_function_applyp:false
  20   /* we do not want FOO(F):=F(1) to mean APPLY(F,[1]) */
  21   )$
  22
  23 infix("~", 134, 133, 'expr, 'expr, 'expr) $
  24 prefix("grad", 142, 'expr, 'expr) $
  25 prefix("div", 142, 'expr, 'expr) $
  26 prefix("curl", 142, 'expr, 'expr) $
  27 prefix("laplacian", 142, 'expr, 'expr) $
  28
  29 declare(order,commutative,
  30         ordern,nary,
  31         ["grad","div","curl","laplacian"],outative,
  32         ["~", "grad", "curl"], nonscalar)$
  33
  34 /* Protect against "circular rule" error in case vect.mac is loaded twice.
  35  * Most of these operators are defined here in vect.mac, so it's unlikely
  36  * to cause trouble by removing their rules. But "." is a built-in operator
  37  * and someone might define a rule for some purpose other than vect.mac,
  38  * so removing rules for "." is a little heavy-handed, but I don't know
  39  * how protect against circular rule errors except by removing rules.
  40  */
  41
  42 /* No easy way to determine if an operator has rules ... Well, this works. */
  43 existing_rules_ops : map (?getop, map (lambda ([x], ?mget (x, '?ruleof)), rules));
  44 has_rules (op) := member (op, existing_rules_ops);
  45
  46 if has_rules (".") then
  47  (print ("vect: warning: removing existing rule or rules for \".\"."),
  48   remove (".", rule));
  49 if has_rules ("~") then remove ("~", rule);
  50 if has_rules ("grad") then remove ("grad", rule);
  51 if has_rules ("div") then remove ("div", rule);
  52 if has_rules ("curl") then remove ("curl", rule);
  53 if has_rules ("laplacian") then remove ("laplacian", rule);
  54 if has_rules (vectorpotential) then remove (vectorpotential, rule);
  55 if has_rules (express) then remove (express, rule);
  56 if has_rules (potential) then remove (potential, rule);
  57
  58 tellsimpafter(0~etrue, zerovector ()) $
  59 tellsimpafter(etrue~0, zerovector ()) $
  60 tellsimpafter(etrue~etrue, zerovector ())$
  61 tellsimpafter(etrue~ttrue.vtrue, etrue.(ttrue~vtrue))$
  62 tellsimp(etrue~lessp, -(lessp~etrue)) $
  63 tellsimpafter(div (curl( etrue)), 0) $
  64 tellsimpafter(curl (grad( etrue)), zerovector ()) $
  65 tellsimpafter(vectorpotential(etrue,ttrue),(scalefactors(ttrue), vpot(etrue)))$
  66 tellsimpafter(vectorpotential(etrue), vpot(etrue)) $
  67 tellsimpafter(express(etrue), express1(etrue)) $
  68 tellsimpafter(express(etrue,ttrue), (scalefactors(ttrue), express(etrue)))$
  69 tellsimpafter(potential(etrue, ttrue),(scalefactors(ttrue),potential1(etrue)))$
  70 tellsimpafter(potential(etrue), potential1(etrue)) $
  71
  72 /* Create a zero vector of the prevailing number of dimensions.
  73  */
  74 zerovector () := makelist (0, dimension);
  75
  76 /* Variables and switches */
  77
  78 define_variable(coordinates, '[x, y, z],any)$
  79 define_variable(dimension,3,fixnum)$
  80 define_variable(dimenimbed,1,fixnum)$
  81 define_variable(trylength,1,fixnum)$
  82 define_variable(bestlength,1,fixnum)$
  83 define_variable(sfprod,1,any)$
  84 define_variable(sf,[1,1,1],list)$
  85
  86 define_variable(expandall,false,boolean)$
  87 define_variable(expanddot,false,boolean)$
  88 define_variable(expanddotplus,false,boolean)$
  89 define_variable(expandgrad,false,boolean)$
  90 define_variable(expandplus,false,boolean)$
  91 define_variable(expandgradplus,false,boolean)$
  92 define_variable(expanddiv,false,boolean)$
  93 define_variable(expanddivplus,false,boolean)$
  94 define_variable(expandcurl,false,boolean)$
  95 define_variable(expandcurlplus,false,boolean)$
  96 define_variable(expandlaplacian,false,boolean)$
  97 define_variable(expandlaplacianplus,false,boolean)$
  98 define_variable(expandprod,false,boolean)$
  99 define_variable(expandgradprod,false,boolean)$
 100 define_variable(expanddivprod,false,boolean)$
 101 define_variable(expandcurlcurl,false,boolean)$
 102 define_variable(expandlaplaciantodivgrad,false,boolean)$
 103 define_variable(expandlaplacianprod,false,boolean)$
 104 define_variable(expandcross,false,boolean)$
 105 define_variable(expandcrosscross,false,boolean)$
 106 define_variable(expandcrossplus,false,boolean)$
 107 define_variable(firstcrossscalar,false,boolean)$
 108
 109 vect_cross: true$
 110
 111
 112 ev_diff(_expr_):=apply('ev,[_expr_,'diff])$
 113
 114
 115 scalefactors(transformation) := block(
 116    if listp(first(transformation)) then (
 117       coordinates: rest(transformation),
 118       transformation: first(transformation))
 119    else coordinates: listofvars(transformation),
 120    dimension: length(coordinates),
 121    dimenimbed: length(transformation),
 122    for row:1 thru dimension do
 123       for col:1 thru dimenimbed do jacobian[row,col]:
 124          trigsimp(ratsimp(diff(transformation[col],
 125             coordinates[row]))),
 126    sfprod:1,
 127    sf : makelist (1, dimension),
 128    for row:1 thru dimension do (
 129       for col:1 thru row-1 do (
 130          sf[row]: gcov(row,col),
 131          if sf[row]#0 then
 132             print("warning: coordinate system is nonorthogonal unless following simplifies to zero:", sf[row])),
 133       sf[row]: radcan(sqrt(gcov(row,row))),
 134       sfprod: sfprod*sf[row])) $
 135
 136 gcov(ii,jj) := trigsimp(ratsimp(sum(
 137    jacobian[ii,kk]*jacobian[jj,kk], kk, 1, dimenimbed))) $
 138
 139 express1(expn) := block(
 140    [ans],
 141    if mapatom(expn) then
 142       if nonscalarp(expn) then (ans:[],
 143          for jj: dimension step -1 thru 1 do
 144             ans: cons(expn[coordinates[jj]], ans),
 145          return(ans))
 146       else return(expn),
 147    expn: map('express1, expn),
 148    if mapatom(expn) or listp(expn) then return(expn),
 149
 150    if inpart(expn,0)="grad" then (ans:[],
 151       expn: inpart(expn,1),
 152       for jj: dimension step -1 thru 1 do ans:
 153          cons('diff(expn,coordinates[jj])/sf[jj], ans),
 154       return(ans)),
 155
 156    if piece="div" then (expn: inpart(expn,1),
 157       if not listp(expn) then error("div called on scalar arg:",
 158          expn),
 159       return(sum('diff(sfprod*expn[jj]/sf[jj],
 160          coordinates[jj]), jj, 1, dimension)/sfprod)),
 161
 162    if piece="laplacian" then return(sum('diff(sfprod*'diff(
 163       inpart(expn,1),coordinates[jj])/sf[jj]**2,
 164       coordinates[jj]), jj, 1, dimension) / sfprod),
 165
 166    if piece="curl" then (expn:inpart(expn,1),
 167       if listp(expn) then (
 168          if length(expn)=2 then return(('diff(sf[2]*expn[2],
 169             coordinates[1])-'diff(sf[1]*expn[1],
 170             coordinates[2]))/ sf[1]/sf[2]),
 171          if dimension=3 then return([
 172              ('diff(sf[3]*expn[3],coordinates[2])-
 173              'diff(sf[2]*expn[2],coordinates[3]))/
 174              sf[2]/sf[3],
 175              ('diff(sf[1]*expn[1],coordinates[3])-
 176               'diff(sf[3]*expn[3],coordinates[1]))/
 177              sf[1]/sf[3],
 178              ('diff(sf[2]*expn[2],coordinates[1]) -
 179               'diff(sf[1]*expn[1],coordinates[2]))/
 180              sf[1]/sf[2]])),
 181       error("curl used in space of wrong dimension")),
 182
 183    if piece="~" then (
 184       ans: inpart(expn,1),  expn:inpart(expn,2),
 185       if listp(ans) and listp(expn) and length(ans)=length(expn)
 186          then (if length(ans)=2 then return(ans[1]*expn[2]
 187              -ans[2]*expn[1]),
 188             if length(ans)=3 then return([ans[2]*expn[3]-
 189                ans[3]*expn[2], ans[3]*expn[1]-ans[1]*expn[3],
 190                ans[1]*expn[2]-ans[2]*expn[1]])),
 191       error("~ used with improper arguments:",ans,expn)),
 192
 193    expn) $
 194
 195
 196 dotassoc: dotexptsimp: false$
 197 dotscrules: true $
 198 define_variable(expandflags, '[
 199    expandall,
 200       expanddot,
 201          expanddotplus,
 202       expandcross,
 203          expandcrossplus,
 204          expandcrosscross,
 205       expandgrad,
 206          expandgradplus,
 207          expandgradprod,
 208       expanddiv,
 209          expanddivplus,
 210          expanddivprod,
 211       expandcurl,
 212          expandcurlplus,
 213          expandcurlcurl,
 214       expandlaplacian,
 215          expandlaplacianplus,
 216          expandlaplacianprod,
 217    expandlaplaciantodivgrad,
 218    expandplus,
 219    expandprod ],list) $
 220
 221
 222 apply('declare, [expandflags, 'evflag]) $
 223
 224 vectorsimp(expn) := block(
 225    [dotdistrib, dotscrules, inflag, firstcrossscalar],
 226    inflag: firstcrossscalar: true,
 227    dotdistrib: expandall or expanddot or expanddotplus
 228       or expandplus,
 229    if expandall or expandgrad or expandgradplus or expandplus
 230       then declare("grad",additive),
 231    if expandall or expanddiv or expanddivplus or expandplus
 232       then declare("div",additive),
 233    if expandall or expandcurl or expandcurlplus or expandplus
 234       then declare("curl",additive),
 235    if expandall or expandlaplacian or expandlaplacianplus or
 236       expandplus then declare("laplacian",additive),
 237    expn: vsimp(expn),
 238    if expandall then expn: ratexpand(expn),
 239    if expandall or expandgrad or expandgradplus or expandplus
 240       then remove("grad",additive),
 241    if expandall or expanddiv or expanddivplus or expandplus
 242       then remove("div",additive),
 243    if expandall or expandcurl or expandcurlplus or expandplus
 244       then remove("curl",additive),
 245    if expandall or expandlaplacian or expandlaplacianplus or
 246       expandplus then remove("laplacian",additive),
 247    expn) $
 248
 249
 250 before(arg) := inpart(order(etrue,arg),1)#etrue$
 251
 252 vscalarp(arg) := not nonscalarp(arg)$
 253
 254 matchdeclare (nn, nonscalarp);
 255 matchdeclare (zz, vector0p);
 256 vector0p(e) := listp(e) and length(e) = dimension and every (lambda ([x], x = 0), e);
 257 defrule (rule_vector0_additive_identity, nn + zz, nn);
 258
 259 matchdeclare (aa, all, bb, ordergreatp (aa));
 260 defrule (rule_dot_commutative, aa . bb, bb . aa);
 261
 262 vsimp(expn) := apply1 (vsimp1(expn), rule_vector0_additive_identity, rule_dot_commutative);
 263
 264 vsimp1(expn) :=
 265    if mapatom(expn) then expn
 266    else block([pv, qv, rv, sv],
 267       expn: map('vsimp1, expn),
 268       if mapatom(expn) then return(expn),
 269       if inpart(expn,0)="~" then expn:removecrosssc1(expn,pv,rv,sv)
 270       else if piece="grad" then (
 271          if (expandall or expandgrad or expandgradprod or
 272             expandprod) and
 273             not mapatom(pv:inpart(expn,1)) and inpart(pv,0)="*"
 274             then expn:apply("+", maplist(lambda([u],\gradprod(u,pv)), pv)))
 275       else if piece="div" then(
 276          if (expandall or expanddiv or expanddivprod or
 277             expandprod) and not
 278             mapatom(pv:inpart(expn,1)) and inpart(pv,0)="*" then
 279             expn: apply("+", maplist(lambda([u],\divprod(u,pv)), pv)))
 280       else if piece="curl" then (
 281          if (expandall or expandcurl or expandcurlcurl) and not
 282             mapatom(pv:inpart(expn,1)) and inpart(pv,0)="curl"
 283             then expn: grad( div(pv:inpart(pv,1))) - laplacian( pv))
 284       else if piece="laplacian" then
 285          if expandlaplaciantodivgrad then
 286             expn: div (grad( inpart(expn,1)))
 287          else if (expandall or expandlaplacian or
 288             expandlaplacianprod or expandprod) and not mapatom
 289             (pv:inpart(expn,1)) and inpart(pv,0)="*" then(
 290             qv: inpart(pv,1),
 291             rv: delete(qv,pv),
 292             expn: rv*laplacian(qv) + 2*grad (rv) * grad( qv) + qv*
 293                laplacian(rv)),
 294       expn) $
 295
 296 crosssimp(ex,pv,rv,sv) :=
 297    if not mapatom(ex) and inpart(ex,0)="~" then (
 298       if expandall or expandcross or expandcrosscross then
 299          ex: trycrosscross(ex,pv,rv,sv),
 300       if not mapatom(ex) and inpart(ex,0)="~" and(
 301          expandall or expandcross or expandcrossplus or
 302          expandplus) then ex: trycrossplus(ex,pv,rv,sv),
 303       ex)
 304    else ex $
 305
 306 removecrosssc(expn,pv,rv,sv) :=
 307    if not mapatom(expn) and inpart(expn,0)="~" then
 308       removecrosssc1(expn,pv,rv,sv)
 309    else expn $
 310
 311 removecrosssc1(expn,pv,rv,sv) :=  block(
 312       [left, right],
 313       left: partitionsc(inpart(expn,1)),
 314       right: partitionsc(inpart(expn,2)),
 315       if firstcrossscalar and (left[2]=1 or right[2]=1)
 316          then( print("warning: declare vector indeterminants
 317 nonscalar to avoid errors & to get full simplification"),
 318             firstcrossscalar:false,
 319                 return(expn)),
 320       left[1]*right[1]*crosssimp(left[2]~right[2],pv,rv,sv))$
 321
 322 partitionsc(ex) :=
 323    if mapatom(ex) then
 324       if nonscalarp(ex) then [1,ex]
 325       else [ex,1]
 326    else if inpart(ex,0)="*" then block([sc,nonsc],
 327       sc: nonsc: 1,
 328       for fact in ex do
 329          if nonscalarp(fact) then nonsc:nonsc*fact
 330          else sc:sc*fact,
 331       [sc,nonsc])
 332    else [1,ex] $
 333
 334 /*
 335 trycrossplus(expn,pv,rv,sv) :=(
 336    pv:inpart(expn,1), rv:inpart(expn,2),
 337    if not mapatom(pv) and inpart(pv,0)="+" then
 338       if not mapatom(rv) and inpart(rv,0)="+" then
 339          map(lambda([u],trycrossplus(u,pv,rv,sv)),
 340          map(lambda([u],crossrv(u,rv,sv)), pv))
 341       else map(lambda([u],crossrv(u,rv,sv)), pv)
 342    else if not mapatom(rv) and inpart(rv,0)="+" then
 343       map(lambda([u],pvcross(pv,u,sv)), rv)
 344    else expn) $
 345 */
 346
 347 trycrossplus(expn,pv,rv,sv) :=
 348   (
 349     pv:inpart(expn,1),
 350     rv:inpart(expn,2),
 351     if not mapatom(pv) and inpart(pv,0)="+" then
 352        map(lambda([u],crossrv(u,rv,sv)), pv)
 353     else
 354        if not mapatom(rv) and inpart(rv,0)="+" then
 355           map(lambda([u],pvcross(pv,u,sv)), rv)
 356        else
 357           expn)$
 358
 359 trycrosscross(expn,pv,rv,sv) := (
 360       pv:inpart(expn,1),  rv:inpart(expn,2),
 361       if not mapatom(rv) and inpart(rv,0)="~" then (
 362          sv: inpart(rv,2), rv:inpart(rv,1),
 363          rv*pv.sv - sv*pv.rv)
 364       else if not mapatom(pv) and inpart(pv,0)="~" then(
 365          sv:inpart(pv,2), pv:inpart(pv,1),
 366          sv*rv.pv - pv*rv.sv)
 367       else expn) $
 368
 369 pvcross(pv,rv,sv) := removecrosssc(pv~rv,pv,rv,sv) $
 370
 371 crossrv(pv,rv,sv) := removecrosssc(pv~rv,pv,rv,sv) $
 372
 373 \gradprod(uu,pv) := delete(uu,pv)*grad(uu) $
 374
 375 \divprod(uu,pv) := block([dotscrules],
 376    dotscrules: false,
 377
 378    if nonscalarp(uu) then delete(uu,pv)*div(uu)
 379    else delete(uu,pv).grad(uu) )$
 380
 381
 382 potential1(gr) := block(
 383    [origin, grperm, jj, result,%dum],
 384    if not listp(gr) or length(gr)#dimension then error(
 385       "1st arg of potential must be a list of length equal to",
 386       "the dimensionality of the coordinate system"),
 387    origin: zeroloc(),
 388    result: [],
 389    for jj:dimension step -1 thru 1 do
 390       result: cons(sf[jj]*gr[jj], result),
 391    grperm:[],
 392    for eqn in origin do (
 393       jj:1,
 394       while lhs(eqn)#coordinates[jj] do jj:jj+1,
 395       grperm: endcons(result[jj], grperm)),
 396    result:sum(myint(subless(jj,origin,grperm), %dum, rhs(origin[jj]),
 397       lhs(origin[jj])), jj, 1, dimension),
 398    gr: gr-express1(grad( result)),
 399    gr: ev_diff(gr),
 400    gr: trigsimp(radcan(gr)),
 401    origin:1,
 402     /* variable name should not be re-used! */
 403    while origin<=dimension and gr[origin]=0 do
 404          (mode_declare(origin,fixnum),origin:origin+1),
 405    if origin<=dimension then print("unable to prove that the",
 406       "following difference between the input and the gradient",
 407       "of the returned result is zero", gr),
 408    trigsimp(radcan(result))) $
 409
 410 define_variable(potentialzeroloc, 0, any) $
 411
 412 zeroloc() :=
 413    if not listp(potentialzeroloc) then
 414        map(lambda([uu],uu=potentialzeroloc), coordinates)
 415    else if disjunct(coordinates,map('lhs,potentialzeroloc)) # []
 416       then error("potentialzeroloc must be a list of length",
 417       "equaling the dimensionality of the coordinate system",
 418       "containing equations with each coordinate variable",
 419       "on the lhs of exactly 1 equation,",
 420       "or else potentialzeroloc must not be a list")
 421    else potentialzeroloc$
 422
 423
 424
 425 eval_when([translate,batch,demo,load,loadfile],
 426
 427 cyc(ii) ::= buildq([ii], 1 + remainder(ii+shift,3)) )$
 428
 429 vpot(kurl) := block(
 430    [origin, shift],
 431    mode_declare(shift,fixnum),
 432    if not listp(kurl) or length(kurl)#3 then error(
 433       "1st arg of vectorpotential must be a list of length 3"),
 434    origin: zeroloc(),
 435    shift: 1,
 436    while shift<=3 and lhs(origin[1])#coordinates[shift] do
 437       shift:shift+1,
 438    shift: shift+1,
 439    if shift>4 or lhs(origin[2])#coordinates[cyc(2)] or
 440       lhs(origin[3])#coordinates[cyc(3)] then error(
 441       "left sides of potentialzeroloc must be a cyclic",
 442       "permutation of coordinates"),
 443    origin: [(myint(sf[cyc(1)]*sf[cyc(3)]*kurl[cyc(2)],
 444       lhs(origin[3]),rhs(origin[3]),lhs(origin[3])) - myint(
 445       sf[cyc(1)]*sf[cyc(2)]*subst(origin[3],kurl[cyc(3)]),
 446       lhs(origin[2]),rhs(origin[2]),lhs(origin[2])))/sf[cyc(1)],
 447       -myint(sf[cyc(2)]*sf[cyc(3)]*kurl[cyc(1)],
 448       lhs(origin[3]),rhs(origin[3]),lhs(origin[3]))/sf[cyc(2)],
 449       0],
 450    origin: [origin[cyc(cyc(1))], origin[cyc(cyc(2))],
 451       origin[cyc(cyc(3))]],
 452    kurl: kurl-express1(curl (origin)),
 453    kurl: ev_diff(kurl),
 454    kurl: trigsimp(radcan(kurl)),
 455    for jj:1 thru 3 do if kurl[jj]#0 then print(
 456       "unable to prove that the following difference between a",
 457       "component of the input and of the curl output is zero",
 458       kurl[jj]),
 459    origin) $
 460
 461
 462
 463 disjunct(l1,l2) := append(setdiff(l1,l2), setdiff(l2,l1)) $
 464
 465 setdiff(l1,l2) :=
 466    if l1=[] then []
 467    else if member(first(l1),l2) then setdiff(rest(l1),l2)
 468    else cons(first(l1), setdiff(rest(l1),l2)) $
 469
 470 subless(kk,origin,grperm) := (mode_declare(kk,fixnum),block([ans,%dum],
 471    ans: ratsubst(%dum, lhs(origin[kk]), grperm[kk]),
 472    for l1:1 thru kk-1 do
 473       ans: ratsubst(rhs(origin[l1]), lhs(origin[l1]),ans),
 474    ans)) $
 475
 476 myint(fun,var,low,high):=block([result,atlow,athigh],
 477   result:integrate(fun,var),
 478   if freeof(nounify('integrate),result) then (
 479         atlow:evlimit(result,var,low),
 480         if atlow=false then go(nogood),
 481         athigh:evlimit(result,var,high),
 482         if athigh=false then go(nogood),
 483         return(radcan(athigh-atlow))),
 484   nogood, defint(fun,var,low,high))$
 485
 486 evlimit(expr,var,lim):=block([temp],
 487   if lim='minf or lim='inf then go(uselimit),
 488   temp:errcatch(subst(lim,var,expr)),
 489   if temp#[] then return(temp[1]),
 490   uselimit, temp:limit(expr,var,lim),
 491   if member(temp,'[inf,minf,und,ind,infinity]) then return(false),
 492   if freeof(nounify('limit),temp) then temp)$
 493