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         "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 vsimp(expn) := apply1 (vsimp1(expn), rule_vector0_additive_identity);
 260
 261 vsimp1(expn) :=
 262    if mapatom(expn) then expn
 263    else block([pv, qv, rv, sv],
 264       expn: map('vsimp1, expn),
 265       if mapatom(expn) then return(expn),
 266       if inpart(expn,0)="~" then expn:removecrosssc1(expn,pv,rv,sv)
 267       else if piece="grad" then (
 268          if (expandall or expandgrad or expandgradprod or
 269             expandprod) and
 270             not mapatom(pv:inpart(expn,1)) and inpart(pv,0)="*"
 271             then expn:apply("+", maplist(lambda([u],\gradprod(u,pv)), pv)))
 272       else if piece="div" then(
 273          if (expandall or expanddiv or expanddivprod or
 274             expandprod) and not
 275             mapatom(pv:inpart(expn,1)) and inpart(pv,0)="*" then
 276             expn: apply("+", maplist(lambda([u],\divprod(u,pv)), pv)))
 277       else if piece="curl" then (
 278          if (expandall or expandcurl or expandcurlcurl) and not
 279             mapatom(pv:inpart(expn,1)) and inpart(pv,0)="curl"
 280             then expn: grad( div(pv:inpart(pv,1))) - laplacian( pv))
 281       else if piece="laplacian" then
 282          if expandlaplaciantodivgrad then
 283             expn: div (grad( inpart(expn,1)))
 284          else if (expandall or expandlaplacian or
 285             expandlaplacianprod or expandprod) and not mapatom
 286             (pv:inpart(expn,1)) and inpart(pv,0)="*" then(
 287             qv: inpart(pv,1),
 288             rv: delete(qv,pv),
 289             expn: rv*laplacian(qv) + 2*grad (rv) * grad( qv) + qv*
 290                laplacian(rv)),
 291       expn) $
 292
 293 crosssimp(ex,pv,rv,sv) :=
 294    if not mapatom(ex) and inpart(ex,0)="~" then (
 295       if expandall or expandcross or expandcrosscross then
 296          ex: trycrosscross(ex,pv,rv,sv),
 297       if not mapatom(ex) and inpart(ex,0)="~" and(
 298          expandall or expandcross or expandcrossplus or
 299          expandplus) then ex: trycrossplus(ex,pv,rv,sv),
 300       ex)
 301    else ex $
 302
 303 removecrosssc(expn,pv,rv,sv) :=
 304    if not mapatom(expn) and inpart(expn,0)="~" then
 305       removecrosssc1(expn,pv,rv,sv)
 306    else expn $
 307
 308 removecrosssc1(expn,pv,rv,sv) :=  block(
 309       [left, right],
 310       left: partitionsc(inpart(expn,1)),
 311       right: partitionsc(inpart(expn,2)),
 312       if firstcrossscalar and (left[2]=1 or right[2]=1)
 313          then( print("warning: declare vector indeterminants
 314 nonscalar to avoid errors & to get full simplification"),
 315             firstcrossscalar:false,
 316                 return(expn)),
 317       left[1]*right[1]*crosssimp(left[2]~right[2],pv,rv,sv))$
 318
 319 partitionsc(ex) :=
 320    if mapatom(ex) then
 321       if nonscalarp(ex) then [1,ex]
 322       else [ex,1]
 323    else if inpart(ex,0)="*" then block([sc,nonsc],
 324       sc: nonsc: 1,
 325       for fact in ex do
 326          if nonscalarp(fact) then nonsc:nonsc*fact
 327          else sc:sc*fact,
 328       [sc,nonsc])
 329    else [1,ex] $
 330
 331 /*
 332 trycrossplus(expn,pv,rv,sv) :=(
 333    pv:inpart(expn,1), rv:inpart(expn,2),
 334    if not mapatom(pv) and inpart(pv,0)="+" then
 335       if not mapatom(rv) and inpart(rv,0)="+" then
 336          map(lambda([u],trycrossplus(u,pv,rv,sv)),
 337          map(lambda([u],crossrv(u,rv,sv)), pv))
 338       else map(lambda([u],crossrv(u,rv,sv)), pv)
 339    else if not mapatom(rv) and inpart(rv,0)="+" then
 340       map(lambda([u],pvcross(pv,u,sv)), rv)
 341    else expn) $
 342 */
 343
 344 trycrossplus(expn,pv,rv,sv) :=
 345   (
 346     pv:inpart(expn,1),
 347     rv:inpart(expn,2),
 348     if not mapatom(pv) and inpart(pv,0)="+" then
 349        map(lambda([u],crossrv(u,rv,sv)), pv)
 350     else
 351        if not mapatom(rv) and inpart(rv,0)="+" then
 352           map(lambda([u],pvcross(pv,u,sv)), rv)
 353        else
 354           expn)$
 355
 356 trycrosscross(expn,pv,rv,sv) := (
 357       pv:inpart(expn,1),  rv:inpart(expn,2),
 358       if not mapatom(rv) and inpart(rv,0)="~" then (
 359          sv: inpart(rv,2), rv:inpart(rv,1),
 360          rv*pv.sv - sv*pv.rv)
 361       else if not mapatom(pv) and inpart(pv,0)="~" then(
 362          sv:inpart(pv,2), pv:inpart(pv,1),
 363          sv*rv.pv - pv*rv.sv)
 364       else expn) $
 365
 366 pvcross(pv,rv,sv) := removecrosssc(pv~rv,pv,rv,sv) $
 367
 368 crossrv(pv,rv,sv) := removecrosssc(pv~rv,pv,rv,sv) $
 369
 370 \gradprod(uu,pv) := delete(uu,pv)*grad(uu) $
 371
 372 \divprod(uu,pv) := block([dotscrules],
 373    dotscrules: false,
 374
 375    if nonscalarp(uu) then delete(uu,pv)*div(uu)
 376    else delete(uu,pv).grad(uu) )$
 377
 378
 379 potential1(gr) := block(
 380    [origin, grperm, jj, result,%dum],
 381    if not listp(gr) or length(gr)#dimension then error(
 382       "1st arg of potential must be a list of length equal to",
 383       "the dimensionality of the coordinate system"),
 384    origin: zeroloc(),
 385    result: [],
 386    for jj:dimension step -1 thru 1 do
 387       result: cons(sf[jj]*gr[jj], result),
 388    grperm:[],
 389    for eqn in origin do (
 390       jj:1,
 391       while lhs(eqn)#coordinates[jj] do jj:jj+1,
 392       grperm: endcons(result[jj], grperm)),
 393    result:sum(myint(subless(jj,origin,grperm), %dum, rhs(origin[jj]),
 394       lhs(origin[jj])), jj, 1, dimension),
 395    gr: gr-express1(grad( result)),
 396    gr: ev_diff(gr),
 397    gr: trigsimp(radcan(gr)),
 398    origin:1,
 399     /* variable name should not be re-used! */
 400    while origin<=dimension and gr[origin]=0 do
 401          (mode_declare(origin,fixnum),origin:origin+1),
 402    if origin<=dimension then print("unable to prove that the",
 403       "following difference between the input and the gradient",
 404       "of the returned result is zero", gr),
 405    trigsimp(radcan(result))) $
 406
 407 define_variable(potentialzeroloc, 0, any) $
 408
 409 zeroloc() :=
 410    if not listp(potentialzeroloc) then
 411        map(lambda([uu],uu=potentialzeroloc), coordinates)
 412    else if disjunct(coordinates,map('lhs,potentialzeroloc)) # []
 413       then error("potentialzeroloc must be a list of length",
 414       "equaling the dimensionality of the coordinate system",
 415       "containing equations with each coordinate variable",
 416       "on the lhs of exactly 1 equation,",
 417       "or else potentialzeroloc must not be a list")
 418    else potentialzeroloc$
 419
 420
 421
 422 eval_when([translate,batch,demo,load,loadfile],
 423
 424 cyc(ii) ::= buildq([ii], 1 + remainder(ii+shift,3)) )$
 425
 426 vpot(kurl) := block(
 427    [origin, shift],
 428    mode_declare(shift,fixnum),
 429    if not listp(kurl) or length(kurl)#3 then error(
 430       "1st arg of vectorpotential must be a list of length 3"),
 431    origin: zeroloc(),
 432    shift: 1,
 433    while shift<=3 and lhs(origin[1])#coordinates[shift] do
 434       shift:shift+1,
 435    shift: shift+1,
 436    if shift>4 or lhs(origin[2])#coordinates[cyc(2)] or
 437       lhs(origin[3])#coordinates[cyc(3)] then error(
 438       "left sides of potentialzeroloc must be a cyclic",
 439       "permutation of coordinates"),
 440    origin: [(myint(sf[cyc(1)]*sf[cyc(3)]*kurl[cyc(2)],
 441       lhs(origin[3]),rhs(origin[3]),lhs(origin[3])) - myint(
 442       sf[cyc(1)]*sf[cyc(2)]*subst(origin[3],kurl[cyc(3)]),
 443       lhs(origin[2]),rhs(origin[2]),lhs(origin[2])))/sf[cyc(1)],
 444       -myint(sf[cyc(2)]*sf[cyc(3)]*kurl[cyc(1)],
 445       lhs(origin[3]),rhs(origin[3]),lhs(origin[3]))/sf[cyc(2)],
 446       0],
 447    origin: [origin[cyc(cyc(1))], origin[cyc(cyc(2))],
 448       origin[cyc(cyc(3))]],
 449    kurl: kurl-express1(curl (origin)),
 450    kurl: ev_diff(kurl),
 451    kurl: trigsimp(radcan(kurl)),
 452    for jj:1 thru 3 do if kurl[jj]#0 then print(
 453       "unable to prove that the following difference between a",
 454       "component of the input and of the curl output is zero",
 455       kurl[jj]),
 456    origin) $
 457
 458
 459
 460 disjunct(l1,l2) := append(setdiff(l1,l2), setdiff(l2,l1)) $
 461
 462 setdiff(l1,l2) :=
 463    if l1=[] then []
 464    else if member(first(l1),l2) then setdiff(rest(l1),l2)
 465    else cons(first(l1), setdiff(rest(l1),l2)) $
 466
 467 subless(kk,origin,grperm) := (mode_declare(kk,fixnum),block([ans,%dum],
 468    ans: ratsubst(%dum, lhs(origin[kk]), grperm[kk]),
 469    for l1:1 thru kk-1 do
 470       ans: ratsubst(rhs(origin[l1]), lhs(origin[l1]),ans),
 471    ans)) $
 472
 473 myint(fun,var,low,high):=block([result,atlow,athigh],
 474   result:integrate(fun,var),
 475   if freeof(nounify('integrate),result) then (
 476         atlow:evlimit(result,var,low),
 477         if atlow=false then go(nogood),
 478         athigh:evlimit(result,var,high),
 479         if athigh=false then go(nogood),
 480         return(radcan(athigh-atlow))),
 481   nogood, defint(fun,var,low,high))$
 482
 483 evlimit(expr,var,lim):=block([temp],
 484   if lim='minf or lim='inf then go(uselimit),
 485   temp:errcatch(subst(lim,var,expr)),
 486   if temp#[] then return(temp[1]),
 487   uselimit, temp:limit(expr,var,lim),
 488   if member(temp,'[inf,minf,und,ind,infinity]) then return(false),
 489   if freeof(nounify('limit),temp) then temp)$
 490