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
  25 "~"(any,body):=block
  26 (
  27   [i,l1,l2,l11,l22,l1s,lk1,div,rp],
  28   rp:false,
  29   if length(first(indices(any)))>0 then rp:true,
  30   if not rp then
  31   (
  32   if listp(cartan_basis) and not member(body,cartan_basis) and
  33                              member(any,cartan_basis) then -body~any
  34   else if not atom(any) and op(any)="+" then
  35     apply(op(any),[part(any,1)~body,rest(any)~body])
  36   else if not atom(body) and op(body)="+" then
  37     apply(op(body),[any~part(body,1),any~rest(body)])
  38   else if not atom(any) and op(any)="-" then -((-any)~body)
  39   else if not atom(body) and op(body)="-" then -(any~(-body))
  40   else if not atom(any) and op(any)="*" then
  41   (
  42     if atom(first(any)) then first(any)*(rest(any)~body)
  43     else if atom(last(any)) then last(any)*(rest(any,-1)~body)
  44     else (first(any)~body)*rest(any)
  45   )
  46   else if not atom(body) and op(body)="*" then
  47   (
  48     if atom(first(body)) then first(body)*(any~rest(body))
  49     else if atom(last(body)) then last(body)*(any~rest(body,-1))
  50     else (any~first(body))*rest(body)
  51   )
  52   else if atom(any) or atom(body) then
  53   (
  54     if any=body then 0
  55     else nounify("~")(any,body)
  56   )
  57   else if not atom(any) and nounify(op(any))=nounify("~") and
  58           not atom(body) and nounify(op(body))=nounify("~") and
  59           (?intersect)(args(any),args(body))#[] then 0
  60   else if not atom(any) and nounify(op(any))=nounify("~") then
  61   (
  62      if member(body, any) then 0
  63      else nounify("~")(any,body)
  64   )
  65   else if not atom(body) and nounify(op(body))=nounify("~") then
  66   (
  67      if member(any, body) then 0
  68      else nounify("~")(any,body)
  69   )
  70   else if atom(any) and nounify(op(body))=nounify("~") and
  71           member(any,body) then 0
  72   else if atom(body) and nounify(op(any))=nounify("~") and
  73           member(body,any) then 0
  74   else rp:true
  75   ),
  76   if rp then
  77   (
  78     l1:first(indices(any)),
  79     l2:first(indices(body)),
  80     if l1=[] then return(any*body)
  81     else if l2=[] then return(any*body)
  82     else
  83     (
  84       l11 : makelist(idummy(), i,1, length(l1)),
  85       l22 : makelist(idummy(), i,1, length(l2)),
  86       l1s:map("=",l1,l11), l2s:map("=",l2,l22),
  87       lk1:append(l1,l2),lk2:append(l11,l22),
  88       any:sublis(l1s,any), body:sublis(l2s,body),
  89       div:if igeowedge_flag then length(l1)!*length(l2)! else length(lk1)!,
  90       factor(canform(contract(expand(kdelta(lk1,lk2)*any*body)))/div)
  91     )
  92   )
  93 );
  94
  95 contind(t):=if mapatom(t) then []
  96         else if ?rpobj(t) then conti(t)
  97         else if op(t)="+" then
  98         (
  99                 indices(t),
 100                 contind(part(t,1))
 101         )
 102         else if op(t)="-" then contind(part(t,1))
 103         else if op(t)="/" then append(contind(part(t,1)),covind(part(t,2)))
 104         else if op(t)="*" then
 105         block
 106         (
 107                 [l:[]],
 108                 map(lambda([x],l:append(l,contind(x))),args(t)),
 109                 l
 110         )
 111         else error("Cannot determine list of contravariant indices.");
 112
 113 covind(t):=if mapatom(t) then []
 114         else if ?rpobj(t) then append(covi(t),deri(t))
 115         else if op(t)="+" then
 116         (
 117                 indices(t),
 118                 covind(part(t,1))
 119         )
 120         else if op(t)="-" then covind(part(t,1))
 121         else if op(t)="/" then append(covind(part(t,1)),contind(part(t,2)))
 122         else if op(t)="*" then
 123         block
 124         (
 125                 [l:[]],
 126                 map(lambda([x],l:append(l,covind(x))),args(t)),
 127                 l
 128         )
 129         else error("Cannot determine list of covariant indices.");
 130
 131 extdiff(forma,[optind]):=if not mapatom(forma) and op(forma)=extdiff then 0
 132 /* A misguided attempt to check for valid forms... DO NOT UNCOMMENT
 133  *      else if (?rpobj)(forma)#false and length(conti(forma))>0
 134  *      then error("A differential form cannot have a contravariant index")
 135  *      else if (?rpobj)(forma)#false and length(covi(forma))>1 and
 136  *              dispsym(forma,length(covi(forma)),0)#
 137  *              [[anti,[makelist(i,i,1,length(covi(forma)))],[]]]
 138  *      then error("Tensor not declared anti(all) in its covariant indices")
 139  */
 140         else if not mapatom(forma) and op(forma)="+" then
 141           sum(apply(extdiff,cons(part(forma,%%i),optind)),%%i,1,length(forma))
 142 /* Another misguided attempt to process products. Fails because it
 143  * separates contractable products, yielding invalid results. DO NOT UNCOMMENT
 144  *      else if not mapatom(forma) and op(forma)="*" then block
 145  *        (
 146  *          [d:(if length(optind)>0 then optind[1] else idummy())],
 147  *          extdiff(first(forma),d)*rest(forma)+
 148  *          first(forma)*extdiff(rest(forma),d)
 149  *        )
 150  */
 151         else block
 152 (
 153   [l11,l1s,lk1,lk2,indd,res,
 154 /* l1:listify(intersect(setify(covi(forma)),setify(first(indices(forma))))),*/
 155    l1:listify(intersect(setify(covind(forma)),setify(first(indices(forma))))),
 156    div,ind1:if length(optind)>0 then optind[1] else idummy()],
 157   if l1=[] then return(idiff(forma,ind1))
 158   else
 159   (
 160     l11:makelist(idummy(),i,1,length(l1)),
 161     l1s:map("=",l1,l11),
 162     lk1:append([ind1],l1),
 163     indd:idummy(),
 164     lk2:append(l11,[indd]),
 165     res:sublis(l1s,forma),
 166     res:contract(canform(ratexpand(kdelta(lk1,lk2)*idiff(res,indd)))),
 167     /* If a frame base is used, the contribution of ifb must be added */
 168     if iframe_flag=true then
 169       for i thru length(l1) do
 170     (
 171       l1s:map("=",[l1[i]],[indd]),
 172       res:res+(if evenp(length(l1)) then 1 else -1)*
 173               'ifb([l1[i],ind1],[indd])*sublis(l1s,forma)
 174     ),
 175     div:if igeowedge_flag then (length(lk1)-1)! else length(lk1)!,
 176     factor(res/div)
 177   )
 178 );
 179
 180 /* VTT: To ensure consistent application of the "first" index in | */
 181 permutator(l):=block
 182 (
 183     [result:1,temp],
 184         for i thru length(l)-1 do if orderlessp(l[i+1],l[i]) then
 185         (temp:l[i+1],l[i+1]:l[i],l[i]:temp,i:0,result:-result),
 186         [result,l]
 187 );
 188
 189 infix("|");
 190 declare("|",additive);
 191 /* VTT: ADDITIVE doesn't always do the trick so we break up sums manually */
 192 "|"(forma,vec):=block
 193 (
 194   [ind1,perm],
 195   forma:distrib(forma),
 196   if not(atom(forma)) and (op(forma)="+" or op(forma)="=") then
 197     apply(op(forma),[part(forma,1)|vec,rest(forma)|vec])
 198 /* When I thought that I'd implement MACSYMA's Cartan package features
 199  * as part of this package, I thought I needed this. But I don't.
 200  *  else if listp(forma) and listp(cartan_basis) then
 201  *  (
 202  *    for i thru length(cartan_basis) do
 203  *      vec:subst(forma[i],cartan_basis[i],-vec),
 204  *    vec
 205  *  )
 206  */
 207   else
 208   (
 209     /* We must not attempt to permute derivative indices. Non-derivative
 210        indices are found as the intersection of free indices, and free
 211        indices of the UNDIFF'd form of forma. However, we still do the
 212        contraction with respect to the first free covariant index which
 213        may be a derivative index. Urgh. */
 214 /*    perm:permutator((?intersect)(first(indices(forma)),
 215                                  first(indices(undiff(forma))))),
 216     if perm[2]=[] then perm:[1,first(indices(forma))],
 217     if perm[2]=[] then 0
 218     else perm[1]*(contract(canform(expand(forma*
 219                  vec([],[first(sort(indices(forma)[1]))])))))*/
 220     ind1:first(sort(first(indices(forma)))),
 221     contract(canform(expand(forma*vec([],[ind1]))))
 222   )
 223 );
 224