14
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
the "|_" for the inner product of a form with a vector and "@L" for the
17
Lie derivative of the form.
15
* exterior forms. We introduce the "~" operator for the exterior product,
16
* and "|" for the inner product of a form with a vector.
19
load("tensor/itensor")$
21
declare("&",additive);
23
"&"(ANY,BODY):=block(local(I),
24
L1:first(indices(ANY)),
25
L2:first(indices(BODY)),
26
if L1=[] then return(ANY*BODY)
27
else if L2=[] then return(ANY*BODY)
28
else ( L11 : MAKELIST(DUMMY(), I,1, LENGTH(L1)),
29
L22 : MAKELIST(DUMMY(), I,1, LENGTH(L2)),
30
L1S:MAP("=",L1,L11), L2S:MAP("=",L2,L22),
31
LK1:append(L1,L2),LK2:append(L11,L22),
32
ANY:sublis(L1S,ANY), BODY:sublis(L2S,BODY),
33
evf:(if evenp(length(LK1)) then 1 else -1),
35
expand(evf*KDELTA(LK1,LK2)*ANY*BODY/(length(LK1)-1)!)))));
38
declare("@",additive);
39
"@"(forma,ind):=block(L1:first(indices(forma)),
40
if L1=[] then return(diff(forma,ind))
41
else ( L11 : MAKELIST(DUMMY(), I,1, LENGTH(L1)),
44
indd:dummy(),LK2:append(L11,[indd]),
45
forma:sublis(L1S,forma),
46
evf:(if evenp(length(LK1)) then 1 else -1),
47
contract(canform(ratexpand(
48
evf*kdelta(LK1,LK2)*diff(forma,indd)/(length(LK1)-1)!)))));
50
/* VTT: To ensure consistent application of the "first" index in |_ */
51
PERMUTATOR(L):=BLOCK([RESULT:1,TEMP],
52
FOR I THRU LENGTH(L)-1 DO IF ORDERLESSP(L[I+1],L[I]) THEN
53
(TEMP:L[I+1],L[I+1]:L[I],L[I]:TEMP,I:0,RESULT:-RESULT),
58
declare("|_",additive);
20
declare("~",additive);
21
declare("~",antisymmetric);
27
[i,l1,l2,l11,l22,l1s,lk1,div],
28
if listp(cartan_basis) and not member(body,cartan_basis) and
29
member(any,cartan_basis) then -body~any
30
else if not atom(any) and op(any)="+" then
31
apply(op(any),[part(any,1)~body,rest(any)~body])
32
else if not atom(body) and op(body)="+" then
33
apply(op(body),[any~part(body,1),any~rest(body)])
34
else if not atom(any) and op(any)="-" then -((-any)~body)
35
else if not atom(body) and op(body)="-" then -(any~(-body))
36
else if not atom(any) and op(any)="*" then
38
if numberp(first(any)) then first(any)*(second(any)~body)
39
else (first(any)~body)*second(any)
41
else if not atom(body) and op(body)="*" then
43
if numberp(first(body)) then first(body)*(any~second(body))
44
else (any~first(body))*second(body)
46
else if atom(any) and atom(body) then
49
else nounify("~")(any,body)
51
else if not atom(any) and nounify(op(any))=nounify("~") and
52
not atom(body) and nounify(op(body))=nounify("~") and
53
(?intersect)(args(any),args(body))#[] then 0
54
else if not atom(any) and nounify(op(any))=nounify("~") then
56
if member(body, any) then 0
57
else nounify("~")(any,body)
59
else if not atom(body) and nounify(op(body))=nounify("~") then
61
if member(any, body) then 0
62
else nounify("~")(any,body)
64
else if atom(any) and nounify(op(body))=nounify("~") and
65
member(any,body) then 0
66
else if atom(body) and nounify(op(any))=nounify("~") and
67
member(body,any) then 0
70
l1:first(indices(any)),
71
l2:first(indices(body)),
72
if l1=[] then return(any*body)
73
else if l2=[] then return(any*body)
76
l11 : makelist(idummy(), i,1, length(l1)),
77
l22 : makelist(idummy(), i,1, length(l2)),
78
l1s:map("=",l1,l11), l2s:map("=",l2,l22),
79
lk1:append(l1,l2),lk2:append(l11,l22),
80
any:sublis(l1s,any), body:sublis(l2s,body),
81
div:if igeowedge_flag then length(l1)!*length(l2)! else length(lk1)!,
82
factor(canform(contract(expand(kdelta(lk1,lk2)*any*body)))/div)
87
extdiff(forma,[optind]):=if not mapatom(forma) and op(forma)=extdiff then 0
88
/* A misguided attempt to check for valid forms... DO NOT UNCOMMENT
89
* else if (?rpobj)(forma)#false and length(conti(forma))>0
90
* then error("A differential form cannot have a contravariant index")
91
* else if (?rpobj)(forma)#false and length(covi(forma))>1 and
92
* dispsym(forma,length(covi(forma)),0)#
93
* [[anti,[makelist(i,i,1,length(covi(forma)))],[]]]
94
* then error("Tensor not declared anti(all) in its covariant indices")
96
else if not mapatom(forma) and op(forma)="+" then
97
sum(apply(extdiff,cons(part(forma,i),optind)),i,1,length(forma))
98
/* Another misguided attempt to process products. Fails because it
99
* separates contractable products, yielding invalid results. DO NOT UNCOMMENT
100
* else if not mapatom(forma) and op(forma)="*" then block
102
* [d:(if length(optind)>0 then optind[1] else idummy())],
103
* extdiff(first(forma),d)*rest(forma)+
104
* first(forma)*extdiff(rest(forma),d)
109
[l11,l1s,lk1,lk2,indd,res,l1:first(indices(forma)),div,
110
ind:if length(optind)>0 then optind[1] else idummy()],
111
if l1=[] then return(idiff(forma,ind))
114
l11:makelist(idummy(),i,1,length(l1)),
116
lk1:append([ind],l1),
118
lk2:append(l11,[indd]),
119
res:sublis(l1s,forma),
120
res:contract(canform(ratexpand(kdelta(lk1,lk2)*idiff(res,indd)))),
121
/* If a frame base is used, the contribution of ifb must be added */
122
if iframe_flag=true then
123
for i thru length(l1) do
125
l1s:map("=",[l1[i]],[indd]),
126
res:res+(if evenp(length(l1)) then 1 else -1)*
127
'ifb([l1[i],ind],[indd])*sublis(l1s,forma)
129
div:if igeowedge_flag then (length(lk1)-1)! else length(lk1)!,
134
/* VTT: To ensure consistent application of the "first" index in | */
138
for i thru length(l)-1 do if orderlessp(l[i+1],l[i]) then
139
(temp:l[i+1],l[i+1]:l[i],l[i]:temp,i:0,result:-result),
144
declare("|",additive);
59
145
/* VTT: ADDITIVE doesn't always do the trick so we break up sums manually */
60
"|_"(forma,vec):=BLOCK(LOCAL(ind,perm),
61
IF NOT(ATOM(forma)) AND (OP(forma)="+" OR OP(forma)="=") THEN
62
APPLY(OP(forma),[PART(forma,1)|_vec,REST(forma)|_vec])
63
ELSE (perm:PERMUTATOR(first(indices(forma))),
65
ELSE perm[1]*(CONTRACT(CANFORM(EXPAND(forma*vec([],[FIRST(perm[2])])))))));
69
declare("@L",additive);
71
"@L"(forma,L1):=block(local(temp1,temp2,ind1),
72
vec:first(L1),ind1:second(L1),
73
if ind1=0 then (ind1:dummy(),V([],[ind1])*diff(forma,ind1))
75
/* VTT: some black magic needed because MAXIMA doesn't handle ADDITIVE as expected */
76
temp1:rename(forma@ind1),
77
temp1:ev(temp1|_vec,noeval),
80
temp2:ev(temp2|_vec,noeval),
81
temp2:ev(temp2,eval)@ind1,
146
"|"(forma,vec):=block
149
forma:distrib(forma),
150
if not(atom(forma)) and (op(forma)="+" or op(forma)="=") then
151
apply(op(forma),[part(forma,1)|vec,rest(forma)|vec])
152
/* When I thought that I'd implement MACSYMA's Cartan package features
153
* as part of this package, I thought I needed this. But I don't.
154
* else if listp(forma) and listp(cartan_basis) then
156
* for i thru length(cartan_basis) do
157
* vec:subst(forma[i],cartan_basis[i],-vec),
163
/* We must not attempt to permute derivative indices. Non-derivative
164
indices are found as the intersection of free indices, and free
165
indices of the UNDIFF'd form of forma. However, we still do the
166
contraction with respect to the first free covariant index which
167
may be a derivative index. Urgh. */
168
/* perm:permutator((?intersect)(first(indices(forma)),
169
first(indices(undiff(forma))))),
170
if perm[2]=[] then perm:[1,first(indices(forma))],
172
else perm[1]*(contract(canform(expand(forma*
173
vec([],[first(sort(indices(forma)[1]))])))))*/
174
ind:first(sort(first(indices(forma)))),
175
contract(canform(expand(forma*vec([],[ind]))))
added
removed
share/tensor/ex_calc.mac
