share/tensor/ex_calc.mac

   1 /* Copyright (C) 2003 Valerij Pipin <pip@iszf.irk.ru>
   2  *
   3  * This program is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU General Public License as
   5  * published by the Free Software Foundation; either version 2 of
   6  * the License, or (at your option) any later version.
   7  *
   8  * This program is distributed in the hope that it will be
   9  * useful, but WITHOUT ANY WARRANTY; without even the implied
  10  * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11  * PURPOSE.  See the GNU General Public License for more details.
  12  *
  13  * Commentary:
  14  * Its purpose is to extend the maxima's itensor ability to manage the
  15  * exterior forms. We introduce the "~" operator for the exterior product,
  16  * and "|" for the inner product of a form with a vector.
  17 */
  18
  19 infix("~");
  20 declare("~",additive);
  21 declare("~",antisymmetric);
  22
  23 igeowedge_flag:false;
  24 __iextdiff_flag:false;
  25
  26 "~"(any,body):=block
  27 (
  28   [i,l1,l2,l11,l22,l1s,lk1,div,rp],
  29   rp:false,
  30   if length(first(indices(any)))>0 then rp:true,
  31   if not rp then
  32   (
  33   if listp(cartan_basis) and not member(body,cartan_basis) and
  34                              member(any,cartan_basis) then -body~any
  35   else if not atom(any) and op(any)="+" then
  36     apply(op(any),[part(any,1)~body,rest(any)~body])
  37   else if not atom(body) and op(body)="+" then
  38     apply(op(body),[any~part(body,1),any~rest(body)])
  39   else if not atom(any) and op(any)="-" then -((-any)~body)
  40   else if not atom(body) and op(body)="-" then -(any~(-body))
  41   else if not atom(any) and op(any)="*" then
  42   (
  43     if atom(first(any)) then first(any)*(rest(any)~body)
  44     else if atom(last(any)) then last(any)*(rest(any,-1)~body)
  45     else (first(any)~body)*rest(any)
  46   )
  47   else if not atom(body) and op(body)="*" then
  48   (
  49     if atom(first(body)) then first(body)*(any~rest(body))
  50     else if atom(last(body)) then last(body)*(any~rest(body,-1))
  51     else (any~first(body))*rest(body)
  52   )
  53   else if atom(any) or atom(body) then
  54   (
  55     if any=body then 0
  56     else nounify("~")(any,body)
  57   )
  58   else if not atom(any) and nounify(op(any))=nounify("~") and
  59           not atom(body) and nounify(op(body))=nounify("~") and
  60           (intersect)(setify(args(any)),setify(args(body)))#[] then 0
  61   else if not atom(any) and nounify(op(any))=nounify("~") then
  62   (
  63      if member(body, any) then 0
  64      else nounify("~")(any,body)
  65   )
  66   else if not atom(body) and nounify(op(body))=nounify("~") then
  67   (
  68      if member(any, body) then 0
  69      else nounify("~")(any,body)
  70   )
  71   else if atom(any) and nounify(op(body))=nounify("~") and
  72           member(any,body) then 0
  73   else if atom(body) and nounify(op(any))=nounify("~") and
  74           member(body,any) then 0
  75   else rp:true
  76   ),
  77   if rp then
  78   (
  79     l1:first(indices(any)),
  80     l2:first(indices(body)),
  81     if l1=[] then return(any*body)
  82     else if l2=[] then return(any*body)
  83     else
  84     (
  85       l11 : makelist(idummy(), i,1, length(l1)),
  86       l22 : makelist(idummy(), i,1, length(l2)),
  87       l1s:map("=",l1,l11), l2s:map("=",l2,l22),
  88       lk1:append(l1,l2),lk2:append(l11,l22),
  89       any:sublis(l1s,any), body:sublis(l2s,body),
  90       div:if igeowedge_flag then length(l1)!*length(l2)! else length(lk1)!,
  91       factor(canform(contract(expand(kdelta(lk1,lk2)*any*body)))/div)
  92     )
  93   )
  94 );
  95
  96 contind(t):=if mapatom(t) then []
  97         else if op(t)=nounify(determinant) then []
  98         else if op(t)=sqrt then []
  99         else if ?rpobj(t) then conti(t)
 100         else if op(t)="+" then
 101         (
 102                 indices(t),
 103                 contind(part(t,1))
 104         )
 105         else if op(t)="-" then contind(part(t,1))
 106         else if op(t)="/" then append(contind(part(t,1)),covind(part(t,2)))
 107         else if op(t)="*" then
 108         block
 109         (
 110                 [l:[]],
 111                 map(lambda([x],l:append(l,contind(x))),args(t)),
 112                 l
 113         )
 114         else error("Cannot determine list of contravariant indices.");
 115
 116 covind(t):=if mapatom(t) then []
 117         else if op(t)=nounify(determinant) then []
 118         else if op(t)=sqrt then []
 119         else if ?rpobj(t) then append(covi(t),deri(t))
 120         else if op(t)="+" then
 121         (
 122                 indices(t),
 123                 covind(part(t,1))
 124         )
 125         else if op(t)="-" then covind(part(t,1))
 126         else if op(t)="/" then append(covind(part(t,1)),contind(part(t,2)))
 127         else if op(t)="*" then
 128         block
 129         (
 130                 [l:[]],
 131                 map(lambda([x],l:append(l,covind(x))),args(t)),
 132                 l
 133         )
 134         else error("Cannot determine list of covariant indices.");
 135
 136 extdiff(forma,[optind]):=if not mapatom(forma) and op(forma)=extdiff then 0
 137 /* A misguided attempt to check for valid forms... DO NOT UNCOMMENT
 138  *      else if (?rpobj)(forma)#false and length(conti(forma))>0
 139  *      then error("A differential form cannot have a contravariant index")
 140  *      else if (?rpobj)(forma)#false and length(covi(forma))>1 and
 141  *              dispsym(forma,length(covi(forma)),0)#
 142  *              [[anti,[makelist(i,i,1,length(covi(forma)))],[]]]
 143  *      then error("Tensor not declared anti(all) in its covariant indices")
 144  */
 145         else if not mapatom(forma) and op(forma)="+" then
 146           sum(apply(extdiff,cons(part(forma,%%i),optind)),%%i,1,length(forma))
 147 /* Another misguided attempt to process products. Fails because it
 148  * separates contractable products, yielding invalid results. DO NOT UNCOMMENT
 149  *      else if not mapatom(forma) and op(forma)="*" then block
 150  *        (
 151  *          [d:(if length(optind)>0 then optind[1] else idummy())],
 152  *          extdiff(first(forma),d)*rest(forma)+
 153  *          first(forma)*extdiff(rest(forma),d)
 154  *        )
 155  */
 156         else block
 157 (
 158   [l1,l11,l1s,lk1,lk2,indd,res,l0:setify(first(indices(forma))),
 159 /* l1:listify(intersect(setify(covi(forma)),setify(first(indices(forma))))),*/
 160    div,ind1:if length(optind)>0 then optind[1] else idummy()],
 161    l1:listify(intersect(setify(covind(forma)),l0)),
 162   if l1=[] then
 163   (
 164     if itorsion_flag = true and length(l0)>0 then return (covdiff(forma,ind1))
 165     else return(idiff(forma,ind1))
 166   )
 167   else if length(l1) >= dim then return(0)
 168   else
 169   (
 170     l11:makelist(idummy(),i,1,length(l1)),
 171     l1s:map("=",l1,l11),
 172     lk1:append([ind1],l1),
 173     indd:idummy(),
 174     lk2:append(l11,[indd]),
 175     res:sublis(l1s,forma),
 176     if itorsion_flag=true then
 177       res:contract(canform(ratexpand(kdelta(lk1,lk2)*covdiff(res,indd))))
 178     else res:contract(canform(ratexpand(kdelta(lk1,lk2)*idiff(res,indd)))),
 179     /* If a frame base is used, the contribution of ifb must be added */
 180     if iframe_flag=true then
 181       for i thru length(l1) do
 182     (
 183       l1s:map("=",[l1[i]],[indd]),
 184       res:res+(if evenp(length(l1)) then 1 else -1)*
 185               'ifb([l1[i],ind1],[indd])*sublis(l1s,forma)
 186     ),
 187     div:if igeowedge_flag then (length(lk1)-1)! else length(lk1)!,
 188     factor(res/div)
 189   )
 190 );
 191
 192 /* VTT: To ensure consistent application of the "first" index in | */
 193 permutator(l):=block
 194 (
 195     [result:1,temp],
 196         for i thru length(l)-1 do if orderlessp(l[i+1],l[i]) then
 197         (temp:l[i+1],l[i+1]:l[i],l[i]:temp,i:0,result:-result),
 198         [result,l]
 199 );
 200
 201 infix("|");
 202 declare("|",additive);
 203 /* VTT: ADDITIVE doesn't always do the trick so we break up sums manually */
 204 "|"(forma,vec):=block
 205 (
 206   [ind1,perm],
 207   forma:distrib(forma),
 208   if not(atom(forma)) and (op(forma)="+" or op(forma)="=") then
 209     apply(op(forma),[part(forma,1)|vec,rest(forma)|vec])
 210 /* When I thought that I'd implement MACSYMA's Cartan package features
 211  * as part of this package, I thought I needed this. But I don't.
 212  *  else if listp(forma) and listp(cartan_basis) then
 213  *  (
 214  *    for i thru length(cartan_basis) do
 215  *      vec:subst(forma[i],cartan_basis[i],-vec),
 216  *    vec
 217  *  )
 218  */
 219   else
 220   (
 221     /* We must not attempt to permute derivative indices. Non-derivative
 222        indices are found as the intersection of free indices, and free
 223        indices of the UNDIFF'd form of forma. However, we still do the
 224        contraction with respect to the first free covariant index which
 225        may be a derivative index. Urgh. */
 226 /*    perm:permutator((?intersect)(first(indices(forma)),
 227                                  first(indices(undiff(forma))))),
 228     if perm[2]=[] then perm:[1,first(indices(forma))],
 229     if perm[2]=[] then 0
 230     else perm[1]*(contract(canform(expand(forma*
 231                  vec([],[first(sort(indices(forma)[1]))])))))*/
 232     ind1:first(sort(first(indices(forma)))),
 233     contract(canform(expand(forma*vec([],[ind1]))))
 234   )
 235 );
 236