Update of /cvsroot/maxima/maxima/share/tensor In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv18256 Modified Files: canten.lisp car_iden.dem ctensor.mac ex_calc.dem ex_calc.mac itensor.lisp kaluza.dem manual.txt swartz.dem symtry.lisp ten_alg.dem tensor-doc.txt Added Files: ctensor1.dem ctensor2.dem ctensor3.dem ctensor4.dem petrov.dem reissner.dem tetrad.dem Removed Files: ajpff.mac ctensor_1.dem ctensor_2.dem ctensor_3.dem ctensor_4.dem jpff.mac kaluza.mac motion.mac petrov.mac test1.mac Log Message: Completed the great renaming: all names now conform to commercial Maxima. Support added for moving frames. Many bugs fixed (and no doubt many new bugs introduced!) --- NEW FILE: ctensor1.dem --- showtime:all$ /* The transformation from Cartesian to spherical coordinates */ /* currently the .mc version is available in DOE MACSYMA */ load("ctensor")$ /* dimension of the manifold */ dim:3$ /* the coordinate labels */ ct_coords:[x,y,z]$ /* the basis metric */ lg:ident(3); (" the required coordinate transformation is: # 1 = X*SIN(Y)*SIN(Z) # 2 = X*SIN(Y)*COS(Z) # 3 = X*COS(Y) PLEASE PROVIDE these right hand sides as INPUT now!!")$ ctransform(lg)$ /* a substitution which reduces the transformed matrix */ ev(%,cos(y) = sqrt(1-sin(y)^2),sin(z) = sqrt(1-cos(z)^2),ratsimp); --- NEW FILE: ctensor2.dem --- showtime:all$ /* This file proves that the Schwarzschild line element satisfies the Einstein vacuum equations as well as computing other tensors and invariants */ /* .mc version is presently available in DOE MACSYMA */ load("ctensor")$ /* the following allows the batch program to run by presetting flags */ setflags()$ /* this calls for the rational simplification of geometrical objects */ ratfac:true$ /* the dimension of the manifold */ dim:4$ /* the coordinate labels */ ct_coords:[x,y,z,t]$ /* THIS IS THE SCHWARZSCHILD METRIC IN STANDARD COORDINATES */ lg: matrix([1/(1-2*m/x),0,0,0],[0,x^2,0,0],[0,0,x^2*sin(y)^2,0],[0,0,0,2*m/x-1]); /* computes metric inverse and determines diagonality */ cmetric()$ /* computes and displays mixed Christoffel symbols */ christof(mcs)$ /* computes and ratsimps Ricci tensor */ uricci(true)$ /* computes scalar curvature */ scurvature(); /* computes Riemann tensor */ lriemann(true)$ /* computes contravariant Riemann tensor */ uriemann(false)$ /* computes the Kretchmann invariant Rijkl^2 */ rinvariant(); --- NEW FILE: ctensor3.dem --- showtime:all$ /* This file proves a simple metric is conformally flat */ /* .mc version currently available in DOE MACSYMA */ load("ctensor")$ setflags()$ ratfac:true$ dim:4$ ct_coords:[x,y,z,t]$ /* THIS IS A CONFORMALLY FLAT METRIC */ lg:a*ident(4)$ dependencies(a(t))$ cmetric()$ christof(mcs)$ uricci(false)$ riemann(false)$ weyl(true)$ --- NEW FILE: ctensor4.dem --- showtime:all$ /* This file finds the Schwarzschild solution of the Einstein vacuum equations */ /* .mc version currently available in DOE MACSYMA */ load("ctensor")$ /* the following allows the batch program to run by presetting flags */ setflags()$ /* this calls for the rational simplification of geometrical objects */ ratfac:true$ /* the dimension of the manifold */ dim:4$ /* the coordinate labels */ ct_coords:[x,y,z,t]$ /* This is the standard spherically symmetric metric */ lg:matrix([a,0,0,0],[0,x^2,0,0],[0,0,x^2*sin(y)^2,0],[0,0,0,-d]); /* functional dependencies */ depends([a,d],x); /* computes inverse metric and specifies diagonality */ cmetric()$ /* computes the mixed Christoffel symbols */ christof(false)$ /* computes and ratsimps Ricci tensor */ uricci(false)$ /* computes and displays the Einstein tensor */ einstein(true); /* makes a list of the non-zero components of the Einstein tensor (EIN) where the 2 indicates the order of the array EIN */ exp:findde(ein,2); /* now begins to solve the field equations */ ode2(last(exp),a,x); /* a kludge to get the solution (the 1,1 component) explicitly */ solve(%,x); resultlist:solve(%,a)$ h:ev(part(resultlist,1),eval); /* to cast the solution into standard form */ h1:h,exp(%c) = 1/(2*m),factor; /* now to find the 4,4 component */ ev(first(exp),h1,diff,factor); ode2(num(%),d,x); expand(radcan(%)); h2:ev(%,%c = 1); /* H1 and H2 should be the solution and to check */ sol:[h1,h2]; exp,sol,diff,ratsimp; --- NEW FILE: petrov.dem --- ("Attempt to compute the Petrov classification of the Kerr metric.")$ ("First, we need to load modules and define the metric:")$ load(ctensor); init_ctensor(); gcd:spmod; ct_coords:[t,r,h,p]$ declare([a,m],constant)$ s(r,h):=r^2+(a*cos(h))^2$ d(r):=a^2-2*m*r+r^2$ dim:4$ fri:matrix( [sqrt(d(r)/s(r,h)), 0, 0, -sqrt(d(r)/s(r,h))*a*sin(h)^2], [0, sqrt(s(r,h)/d(r)), 0, 0], [0, 0, sqrt(s(r,h)), 0 ], [-a/sqrt(s(r,h))*sin(h), 0, 0, (s(r,h)+(a*sin(h))^2)*sin(h)/sqrt(s(r,h))]); ("We need aggressive simplification to ensure that the result is correct.")$ ctrgsimp:true; ratfac:true; cframe_flag:true; cmetric(false); ug:invert(lg)$ nptetrad(false); christof(false); lriemann(false); ricci(false); weyl(false); /* ("Since we use a tetrad frame, the Weyl tensor was computed in the tetrad base. We need the Weyl tensor in the coordinate base to compute the Newman-Penrose coefficients, so we must do a conversion first. Please be patient; this is not very efficient and will take a while:")$ for i thru dim do for j thru dim do for k thru dim do for l thru dim do w2[i,j,k,l]:weyl[i,j,k,l]; for i thru dim do for j thru dim do for k thru dim do for l thru dim do weyl[i,j,k,l]:sum(sum(sum(sum(w2[ii,jj,kk,ll]*fri[ii,i]*fri[jj,j]* fri[kk,k]*fri[ll,l],ii,1,dim),jj,1,dim),kk,1,dim),ll,1,dim); */ ("Now we're ready to compute the Newman-Penrose coefficients:")$ psi(true); petrov(); --- NEW FILE: reissner.dem --- ("Verifying that the Reissner-Nordstrom metric is that of a point charge.")$ ("First, we load ctensor...")$ load(ctensor); ("We also need NCHRPL that contains the MATTRACE function:")$ load(nchrpl); ("Now we set up the metric with (+,-,-,-) signature:")$ dim:4; ct_coords:[t,r,o,f]; lg:matrix([(r^2+q^2-2*m*r)/r^2,0,0,0],[0,r^2/(2*m*r-r^2-q^2),0,0],[0,0,-r^2,0],[0,0,0,-r^2*sin(o)^2]); ug:ratsimp(invert(lg)); ("Next, we compute the Christoffel-symbols:")$ christof(false); ("Before computing the Ricci and Einstein tensors, we create them as matrices; this will make life easier later on:")$ ric:lein:ident(4)$ ricci(false); leinstein(false); ("Part II: we set up the electromagnetic 4-potential:")$ a:[q/r,0,0,0]; ("From which we compute the covariant electromagnetic field tensor:")$ fl:ident(4)$ for i thru dim do for j thru dim do fl[i,j]:diff(a[i],ct_coords[j])-diff(a[j],ct_coords[i]); ("With the help of the metric, we can obtain the mixed and contravariant electromagnetic field tensors:")$ fm:ug.fl; fu:fm.ug; ("And now we can compute the stress-energy tensor (times 8*%PI).")$ t:ratsimp(2*(fl.fm-lg*mattrace(fl.fu)/4)); ("This and the covariant Einstein tensor should be identical:")$ ratsimp(lein-t); --- NEW FILE: tetrad.dem --- load("ctensor")$ ("First, we enter the initial parameters, starting with the coordinate base:")$ init_ctensor(); ct_coords:[t,r,d,c]$ ("We use dimension 4. The four base vectors arranged in matrix form:")$ fri:matrix( [sqrt(f(r)), 0, 0, 0], [0, sqrt(h(r)), 0, 0], [0, 0, r, 0], [0, 0, 0,r*sin(d)]); ("The problem is simple enough for automatic simplification:")$ ctrgsimp:true; ratfac:true; cframe_flag:true; cmetric(); ("Compute the Ricci rotation coefficients")$ christof(true); ("And the Riemann tensor in tetrad base")$ lriemann(true); ("Now compute the Ricci tensor")$ ricci(true); ("For a vacuum metric, the components of the Ricci tensor must be 0")$ ("We begin by adding ric1,1 and ric2,2 which yields a simple equation:")$ ric[2,2]+ric[1,1]; ("Some algebraic manipulation yields a solution for f(r) in terms of h(r):")$ factor(%th(2)); num(%); %/f(r)/h(r); distrib(%); first(%)=-last(%); integrate(%,r); solve(%,f(r)); ("We can choose the constant to be 1 by rescaling the metric:")$ solf:subst(1,num(last(%th(2)[1])),last(%th(2)[1])); ("We substitute the solution for f(r) into ric3,3:")$ ric[3,3]; ("Further algebraic manipulation helps get a solution for h(r):")$ ratsimp(%th(2)); subst(solf,f(r),%); ratsimp(%); num(%); %,diff; ode2(%,h(r),r); exp(first(%))=exp(last(%)); solve(%,h(r)); ("We rename the constant to get the result in the usual form:")$ solh:ratsimp(subst(1/2/m,exp(%c),last(%th(2)[1]))); ("Now substitute the solutions back into our base vectors:")$ subst(solf,f(r),fri); fri:subst(solh,h(r),%); ("And compute the metric, which is the standard Schwarzschild metric:")$ cmetric(true); ("We're done... now let's verify our solution.")$ ("We begin with recomputing the Ricci rotation coefficients:")$ christof(true); ("And the Riemann tensor:")$ lriemann(true); ("And last, the Ricci tensor:")$ ricci(true); Index: canten.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/canten.lisp,v retrieving revision 1.6 retrieving revision 1.7 diff -u -d -r1.6 -r1.7 --- canten.lisp 18 Nov 2004 03:53:28 -0000 1.6 +++ canten.lisp 25 Nov 2004 01:43:29 -0000 1.7 @@ -4,7 +4,7 @@ -(DECLARE-TOP (SPECIAL FREI BOUNI $CANTERM BREAKLIST SMLIST $DUMMYX)) +(DECLARE-TOP (SPECIAL FREI BOUNI $CANTERM BREAKLIST SMLIST $IDUMMYX)) (SETQ NODOWN '($ICHR2 $ICHR1 %ICHR2 %ICHR1 $KDELTA %KDELTA)) @@ -422,7 +422,7 @@ (DEFUN TFINDSTRUC (FACT) (APPEND (CDADR FACT) (CDADDR FACT) (CDDDR FACT) )) -(DEFUN DUMM (X) (EQUAL (CADR (EXPLODEC X)) $DUMMYX)) +(DEFUN DUMM (X) (EQUAL (CADR (EXPLODEC X)) $IDUMMYX)) (DEFUN FINDPERMUT (I1 I2) Index: car_iden.dem =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/car_iden.dem,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- car_iden.dem 18 Nov 2004 03:53:28 -0000 1.3 +++ car_iden.dem 25 Nov 2004 01:43:29 -0000 1.4 @@ -2,13 +2,11 @@ load("tensor/ex_calc.mac"); dummyx:z; allsym:false; -/*decsym(a,2,0,[ANTI(ALL)],[]); -ishow(liediff(extdiff(a([i1,i2]),i4),[v,i5])-extdiff(liediff(a([i1,i2]),[v,i4]),i5))$*/ +decsym(a,2,0,[anti(all)],[]); +canform(liediff(v,extdiff(a([i1,i2]),i3))-extdiff(liediff(v,a([i1,i2])),i3)); decsym(q,3,0,[anti(all)],[]); -ishow(liediff(extdiff(q([i1,i2,i3]),i4),[v,i5])-extdiff(liediff(q([i1,i2,i3]),[v,i4]),i5))$ +canform(liediff(v,extdiff(q([i1,i2,i3]),i4))-extdiff(liediff(v,q([i1,i2,i3])),i4)); decsym(p,4,0,[anti(all)],[]); -ishow(liediff(extdiff(p([i1,i2,i3,i4]),i5),[v,i6])-extdiff(liediff(p([i1,i2,i3,i4]),[v,i5]),i6))$ -/* +canform(liediff(v,extdiff(p([i1,i2,i3,i4]),i5))-extdiff(liediff(v,p([i1,i2,i3,i4])),i5)); decsym(y,5,0,[anti(all)],[]); -ishow(liediff(extdiff(y([i1,i2,i3,i4,i5]),i6),[v,i7])-extdiff(liediff(y([i1,i2,i3,i4,i5]),[v,i6]),i7))$ -*/ +canform(liediff(v,extdiff(y([i1,i2,i3,i4,i5]),i6))-extdiff(liediff(v,y([i1,i2,i3,i4,i5])),i6)); Index: ctensor.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/ctensor.mac,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- ctensor.mac 18 Nov 2004 03:53:28 -0000 1.2 +++ ctensor.mac 25 Nov 2004 01:43:29 -0000 1.3 @@ -1,9 +1,64 @@ /*-*MACSYMA-*-*/ -/* LOAD OF PACKG REMOVED */ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR [...1741 lines suppressed...] - 2*sum(ug[bb,cc]*( - diff(diff(tracer,ct_coords[aa]),ct_coords[cc]) - -sum(mcs[aa,cc,kk] - *diff(tracer,ct_coords[kk]), - kk,1,dim)),cc,1,dim))))$ - -invariant2():="NOT YET IMPLEMENTED"; -bimetric():="NOT YET IMPLEMENTED"; +/* Covariant divergence. The object gxyz must be a mixed 2nd rank tensor + whose first index is covariant. */ +checkdiv(gxyz):=block +( + modedeclare([i,j],fixnum), + if diagmatrixp(gxyz,dim) then for i thru dim do + print(div[i]:ratsimp(ctaylor(1/sqrt(-gdet)*diff(sqrt(-gdet)*gxyz[i,i], + ct_coords[i])-sum(mcs[i,j,j]*gxyz[j,j],j,1,dim)))) + else for i thru dim do + print(div[i]:ratsimp(ctaylor(1/sqrt(-gdet)*sum(diff(sqrt(-gdet)*gxyz[i,j], + ct_coords[j]),j,1,dim)-sum(sum(mcs[i,j,a]*gxyz[a,j],a,1,dim),j,1,dim)))) +)$ Index: ex_calc.dem =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/ex_calc.dem,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- ex_calc.dem 18 Nov 2004 03:53:28 -0000 1.3 +++ ex_calc.dem 25 Nov 2004 01:43:29 -0000 1.4 @@ -2,57 +2,55 @@ ("the operations are allowed only on the covariant tensors")$ load("ex_calc.mac"); dummyx:z; -("The exterior product is denoted by &. Take it on two 1-forms - a([i])&b([j])")$ +("The exterior product is denoted by ~. Take it on two 1-forms + a([i])~b([j])")$ -ishow(a([i])&b([j]))$ +ishow(a([i])~b([j]))$ ("The exterior product of the three 1-forms - a([i])&b([j])&c([k])")$ + a([i])~b([j])~c([k])")$ -ishow(a([i])&b([j])&c([k]))$ +ishow(a([i])~b([j])~c([k]))$ ("Take a sum of two 3 forms")$ -ishow( factor(a([i])&(bk*b([j]))&c([k]) +(ak*a([i]))&c([j])&b([k])))$ +ishow( factor(a([i])~(bk*b([j]))~c([k]) +(ak*a([i]))~c([j])~b([k])))$ -("declare 2-form (don't forget to put allsym:false) and take & with 1-form")$ +("declare 2-form (don't forget to put allsym:false) and take ~ with 1-form")$ allsym:false; decsym(p,2,0,[anti(all)],[]); -ishow(q([i])&p([j,k]))$ +ishow(q([i])~p([j,k]))$ ("the exterior (anticommutative) derivative is denoted by extdiff(x,ind), where ind denotes the component")$ -("So actually d_k & form means in our notation extdiff(form, k), e.g.")$ +("So actually d_k ~ form means in our notation extdiff(form, k), e.g.")$ ishow(extdiff(p([j,k]),i))$ -ishow(extdiff((a([j])&b([k])),i))$ -ishow(extdiff((a([j])&b([k])),k))$ +ishow(extdiff((a([j])~b([k])),i))$ +ishow(extdiff((a([j])~b([k])),k))$ -("the interior product is denoted by |_, say form |_ vector, where +("the interior product is denoted by |, say form | vector, where form is indexed tensor and vector is vector's name")$ -ishow(a([i])|_a)$ +ishow(a([i])|a)$ -ishow((a([i])&b([j]))|_a)$ +ishow((a([i])~b([j]))|a)$ ("to avoid the mistake please use the literally sorted indices - when you apply the |_ ")$ + when you apply the | ")$ -ishow(factor((a([i2])&b([i1]))|_a+ (a([i1])&b([i2]))|_a))$ +ishow(factor((a([i2])~b([i1]))|a+ (a([i1])~b([i2]))|a))$ -("the Lie derivative is denoted by liediff(x,[vector,ind]), (say liediff(form,[vector,ind])) - where vector is vector name and ind is the index of the component")$ -("it is currently not applicable to functions. - Also you have to use the literally sorted indices")$ +("the Lie derivative is denoted by liediff(vector,x)")$ +("it is currently not applicable to functions.")$ -ishow(liediff(a([i1]),[v,i2]))$ +ishow(liediff(v,a([i1])))$ ("tests the consequence of the Cartan identity for 1-forms. The Lie derivative has to commute with the exterior one ")$ -ishow(liediff(extdiff(a([i1]),i2),[v,i3])-extdiff(liediff(a([i1]),[v,i2]),i3))$ +ishow(canform(liediff(v,extdiff(a([i1]),i2))-extdiff(liediff(v,a([i1])),i2)))$ -("where the liediff(extdiff(a([i1]),i2),[v,i3]) is")$ -ishow(liediff(extdiff(a([i1]),i2),[v,i3]))$ +("where the liediff(v,extdiff(a([i1]),i2)) is")$ +ishow(liediff(v,extdiff(a([i1]),i2)))$ ("the verifications of this consequence of the Cartan identity for the higher order forms are in car_iden.dem")$ Index: ex_calc.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/ex_calc.mac,v retrieving revision 1.8 retrieving revision 1.9 diff -u -d -r1.8 -r1.9 --- ex_calc.mac 18 Nov 2004 12:36:13 -0000 1.8 +++ ex_calc.mac 25 Nov 2004 01:43:29 -0000 1.9 @@ -12,63 +12,138 @@ * * Commentary: * Its purpose is to extend the maxima's itensor ability to manage the - exterior forms. We introduce the "&" operator for the exterior product, - and "|_" for the inner product of a form with a vector. + * exterior forms. We introduce the "~" operator for the exterior product, + * and "|" for the inner product of a form with a vector. */ load("tensor/itensor")$ -infix("&"); -declare("&",additive); -"&"(any,body):=block(local(i), - l1:first(indices(any)), - l2:first(indices(body)), - if l1=[] then return(any*body) - else if l2=[] then return(any*body) - else ( l11 : makelist(idummy(), i,1, length(l1)), - l22 : makelist(idummy(), i,1, length(l2)), - l1s:map("=",l1,l11), l2s:map("=",l2,l22), - lk1:append(l1,l2),lk2:append(l11,l22), - any:sublis(l1s,any), body:sublis(l2s,body), - canform(contract( - expand(kdelta(lk1,lk2)*any*body/(length(lk1)-1)!))))); +infix("~"); +declare("~",additive); +declare("~",antisymmetric); -extdiff(forma,ind):=block(l1:first(indices(forma)), - if l1=[] then return(diff(forma,ind)) - else ( l11 : makelist(idummy(), i,1, length(l1)), +"~"(any,body):=block +( + [i,l1,l2,l11,l22,l1s,lk1], + if listp(cartan_basis) and not member(body,cartan_basis) and + member(any,cartan_basis) then -body~any + else if not atom(any) and op(any)="+" then + apply(op(any),[part(any,1)~body,rest(any)~body]) + else if not atom(body) and op(body)="+" then + apply(op(body),[any~part(body,1),any~rest(body)]) + else if not atom(any) and op(any)="-" then -((-any)~body) + else if not atom(body) and op(body)="-" then -(any~(-body)) + else if not atom(any) and op(any)="*" then + ( + if numberp(first(any)) then first(any)*(second(any)~body) + else (first(any)~body)*second(any) + ) + else if not atom(body) and op(body)="*" then + ( + if numberp(first(body)) then first(body)*(any~second(body)) + else (any~first(body))*second(body) + ) + else if atom(any) and atom(body) then + ( + if any=body then 0 + else nounify("~")(any,body) + ) + else if not atom(any) and nounify(op(any))=nounify("~") and + not atom(body) and nounify(op(body))=nounify("~") and + (?intersect)(args(any),args(body))#[] then 0 + else if not atom(any) and nounify(op(any))=nounify("~") then + ( + if member(body, any) then 0 + else nounify("~")(any,body) + ) + else if not atom(body) and nounify(op(body))=nounify("~") then + ( + if member(any, body) then 0 + else nounify("~")(any,body) + ) + else if atom(any) and nounify(op(body))=nounify("~") and + member(any,body) then 0 + else if atom(body) and nounify(op(any))=nounify("~") and + member(body,any) then 0 + else + ( + l1:first(indices(any)), + l2:first(indices(body)), + if l1=[] then return(any*body) + else if l2=[] then return(any*body) + else + ( + l11 : makelist(idummy(), i,1, length(l1)), + l22 : makelist(idummy(), i,1, length(l2)), + l1s:map("=",l1,l11), l2s:map("=",l2,l22), + lk1:append(l1,l2),lk2:append(l11,l22), + any:sublis(l1s,any), body:sublis(l2s,body), + canform(contract(expand(kdelta(lk1,lk2)*any*body)))/(length(lk1)-1)! + ) + ) +); + +extdiff(forma,[optind]):=block +( + [l1:first(indices(forma)), + ind:if length(optind)>0 then optind[1] else idummy()], + if l1=[] then return(idiff(forma,ind)) + else + ( + l11:makelist(idummy(),i,1,length(l1)), l1s:map("=",l1,l11), lk1:append([ind],l1), - indd:idummy(),lk2:append(l11,[indd]), + indd:idummy(), + lk2:append(l11,[indd]), forma:sublis(l1s,forma), - contract(canform(ratexpand( - kdelta(lk1,lk2)*diff(forma,indd)/(length(lk1)-1)!))))); + contract(canform(ratexpand(kdelta(lk1,lk2)*idiff(forma,indd))))/ + (length(lk1)-1)! + ) +); -/* VTT: To ensure consistent application of the "first" index in |_ */ -permutator(l):=block([result:1,temp], +/* VTT: To ensure consistent application of the "first" index in | */ +permutator(l):=block +( + [result:1,temp], for i thru length(l)-1 do if orderlessp(l[i+1],l[i]) then (temp:l[i+1],l[i+1]:l[i],l[i]:temp,i:0,result:-result), [result,l] ); -infix("|_"); -declare("|_",additive); +infix("|"); +declare("|",additive); /* VTT: ADDITIVE doesn't always do the trick so we break up sums manually */ -"|_"(forma,vec):=block(local(ind,perm), +"|"(forma,vec):=block +( + [ind,perm], if not(atom(forma)) and (op(forma)="+" or op(forma)="=") then - apply(op(forma),[part(forma,1)|_vec,rest(forma)|_vec]) - else (perm:permutator(first(indices(forma))), + apply(op(forma),[part(forma,1)|vec,rest(forma)|vec]) + else if listp(forma) and listp(cartan_basis) then + ( + for i thru length(cartan_basis) do + vec:subst(forma[i],cartan_basis[i],-vec), + vec + ) + else + ( + perm:permutator(first(indices(forma))), if perm[2]=[] then 0 - else perm[1]*(contract(canform(expand(forma*vec([],[first(perm[2])]))))))); + else perm[1]*(contract(canform(expand(forma*vec([],[first(perm[2])]))))) + ) +); + + +/* liediff(forma,l1):=block(local(temp1,temp2,ind1), vec:first(l1),ind1:second(l1), if ind1=0 then (ind1:idummy(),v([],[ind1])*diff(forma,ind1)) else( /* VTT: some black magic needed because MAXIMA doesn't handle ADDITIVE as expected */ temp1:rename(extdiff(forma,ind1)), - temp1:ev(temp1|_vec,noeval), + temp1:ev(temp1|vec,noeval), temp1:ev(temp1,eval), temp2:rename(forma), - temp2:ev(temp2|_vec,noeval), + temp2:ev(temp2|vec,noeval), temp2:extdiff(ev(temp2,eval),ind1), temp1+temp2)); - +*/ Index: itensor.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/itensor.lisp,v retrieving revision 1.18 retrieving revision 1.19 diff -u -d -r1.18 -r1.19 --- itensor.lisp 18 Nov 2004 03:53:28 -0000 1.18 +++ itensor.lisp 25 Nov 2004 01:43:29 -0000 1.19 @@ -57,16 +57,16 @@ #+maclisp ($UUO) ;Restore calls to SDIFF so it can be redefined -(DECLARE-TOP (SPECIAL SMLIST $DUMMYX $COORDINATES $IMETRIC $ICOUNTER $DIM +(DECLARE-TOP (SPECIAL SMLIST $IDUMMYX $VECT_COORDS $IMETRIC $ICOUNTER $DIM $CONTRACTIONS $COORD $ALLSYM $METRICCONVERT) (*LEXPR $RENAME $DIFF $COORD $REMCOORD $LORENTZ_GAUGE)) -(SETQ $DUMMYX '$% ;Prefix for dummy indices +(SETQ $IDUMMYX '$% ;Prefix for dummy indices $ICOUNTER 0. ;Dummy variable numeric indexs SMLIST '(MLIST SIMP) ;Simplified MLIST header - $COORDINATES NIL ;Used when differentiating w.r.t. a number + $VECT_COORDS NIL ;Used when differentiating w.r.t. a number $COORD '((MLIST SIMP)) ;Objects treated liked coordinates in DIFF - $ALLSYM T ;If T then all indexed objects symmetric + $ALLSYM NIL ;If T then all indexed objects symmetric $METRICCONVERT T) ;Flag used by $IC_CONVERT ;(DEFUN IFNOT MACRO (CLAUSE) (CONS 'OR (CDR CLAUSE))) @@ -80,7 +80,7 @@ (DEFMFUN $IDUMMY nil ;Sets arguments to dummy indices (progn (setq $ICOUNTER (1+ $ICOUNTER)) - (concat $DUMMYX $ICOUNTER))) + (concat $IDUMMYX $ICOUNTER))) (DEFPROP $KDELTA ((/ . / )) CONTRACTIONS) @@ -201,7 +201,7 @@ (PROG (A B C) (COND ; ((> NARGS 2) (RETURN (MEVAL (CONS '$COVDIFF (CONS ($ICHR1 (ARG 1) (ARG 2)) (CDDR (LISTIFY NARGS))))))) - ((> NARGS 2) (RETURN (MEVAL (CONS '$DIFF (CONS ($ICHR1 (ARG 1) (ARG 2)) (APPLY #'APPEND (MAPCAR #'(LAMBDA (E) (LIST E 1)) (CDDR (LISTIFY NARGS))))))))) + ((> NARGS 2) (RETURN (MEVAL (CONS '$IDIFF (CONS ($ICHR1 (ARG 1) (ARG 2)) (APPLY #'APPEND (MAPCAR #'(LAMBDA (E) (LIST E 1)) (CDDR (LISTIFY NARGS))))))))) ((> NARGS 1) (AND (EQ 1 (LENGTH (ARG 2))) (RETURN ($ICHR1 (ARG 1)))) (merror "ICHR1 cannot have contravariant indices")) ;; (COND @@ -237,7 +237,7 @@ (DEFMFUN $ICHR2 NARGS (PROG (A B C D) (COND - ((> NARGS 2) (RETURN (MEVAL (CONS '$DIFF (CONS ($ICHR2 (ARG 1) (ARG 2)) (APPLY #'APPEND (MAPCAR #'(LAMBDA (E) (LIST E 1)) (CDDR (LISTIFY NARGS))))))))) + ((> NARGS 2) (RETURN (MEVAL (CONS '$IDIFF (CONS ($ICHR2 (ARG 1) (ARG 2)) (APPLY #'APPEND (MAPCAR #'(LAMBDA (E) (LIST E 1)) (CDDR (LISTIFY NARGS))))))))) (T (SETQ A (CADR (ARG 1)) B (CADDR (ARG 1)) C (CADR (ARG 2))) (return (do ((flag) (l (append (cdr (ARG 1)) (cdr (ARG 2))))) @@ -254,13 +254,13 @@ (setq r ($idummy)) (SETQ I (CADR L1) K (CADDR L1) H (CADDDR L1) J (CADR L2)) (RETURN (LIST '(MPLUS) - (SDIFF (LIST '($ICHR2 SIMP) + (IDIFF (LIST '($ICHR2 SIMP) (LIST SMLIST I K) L2) H) (LIST '(MTIMES) -1. - (SDIFF (LIST '($ICHR2 SIMP) + (IDIFF (LIST '($ICHR2 SIMP) (LIST SMLIST I H) (LIST SMLIST J)) K)) @@ -305,7 +305,7 @@ (DEFUN COVDIFF (E) (setq d ($idummy)) (COND - ((OR (ATOM E) (EQ (CAAR E) 'RAT)) (SDIFF E X)) + ((OR (ATOM E) (EQ (CAAR E) 'RAT)) (IDIFF E X)) ((RPOBJ E) (SETQ TEMP (MAPCAR #'(LAMBDA (V) (LIST '(MTIMES) (LIST '($ICHR2 SIMP) @@ -317,7 +317,7 @@ (CONS '(MPLUS) (CONS - (SDIFF E X) + (IDIFF E X) (COND ((OR (CDADR E) (CDDDR E)) (CONS @@ -650,19 +650,132 @@ (AND E ; (DO E E (CDR E) ; (NULL E) -; (SETQ F (SDIFF F (CAR E)))) +; (SETQ F (IDIFF F (CAR E)))) (DO ((E E (CDR E))) ((NULL E) ) - (SETQ F (SDIFF F (CAR E)))) + (SETQ F (IDIFF F (CAR E)))) ) (RETURN F))) + +(defun liediff (v e n) + (cond + ((not (symbolp v)) (merror "~M is not a symbol" v)) + ((or (atom e) (eq (caar e) 'rat)) ; Scalar field + ; v([],[%1])*idiff(e,%1) + (let ((dummy (implode (nconc (exploden $idummyx) (exploden n))))) + (list '(mtimes) + (list (list v) '((mlist)) (list '(mlist) dummy)) + ($idiff e dummy) + ) + ) + ) + ((rpobj e) ; Tensor field +; Dummy implementation for logic tests +; (list '(%liediff) v e) +; Shall the dummy index be in ICOUNTER sequence? Probably yes. +; (let ((dummy (implode (nconc (exploden $idummyx) (exploden n))))) + (let ((dummy ($idummy))) + (append + (list '(mplus) 0 + (list '(mtimes) ; e([...],[...],%1)*v([],[%1]) + (list (list v) '((mlist)) (list '(mlist) dummy)) + ($idiff e dummy) + ) + ) + (maplist #'(lambda (s) ; e([..%1..],[...])*v([],[%1],k) + (list '(mtimes) + (append (list (car e) + (cons '(mlist) + (append (subseq (cdadr e) 0 + (- (length (cdadr e)) (length s)) + ) + (cons dummy (cdr s)) + ) + ) + (caddr e) + ) + (cdddr e) + ) + (list (list v) '((mlist)) (list '(mlist) dummy) (car s)) + ) + ) + (cdadr e) + ) + (maplist #'(lambda (s) ; +e([...],[...],..%1..)*v([],[%1],k) + (list '(mtimes) + (append (list (car e) + (cadr e) + (caddr e) + ) + (subseq (cdddr e) 0 + (- (length (cdddr e)) (length s)) + ) + (cons dummy (cdr s)) + ) + (list (list v) '((mlist)) (list '(mlist) dummy) (car s)) + ) + ) + (cdddr e) + ) + (maplist #'(lambda (s) ; -e([...],[..%1..])*v([],[k],%1) + (list '(mtimes) -1 + (append (list (car e) + (cadr e) + (cons '(mlist) + (append (subseq (cdaddr e) 0 + (- (length (cdaddr e)) (length s)) + ) + (cons dummy (cdr s)) + ) + ) + ) + (cdddr e) + ) + (list (list v) '((mlist)) (list '(mlist) (car s)) dummy) + ) + ) + (cdaddr e) + ) + ) + ) + ) + ((eq (caar e) 'mtimes) ; Leibnitz rule + ; Lv(cadr e)*(cddr e)+(cadr e)*Lv(cddr e) + (list '(mplus) + (cons '(mtimes) (cons (liediff v (cadr e) n) (cddr e))) + (cons '(mtimes) + (list (cadr e) (liediff v + (cond ((cdddr e) (cons '(mtimes) (cddr e))) + (t (caddr e))) + n))) + ) + ) + ((eq (caar e) 'mplus) ; Linearity +; We prefer mapcar to iteration, but the commented code also works +; (list '(mplus) (liediff v (cadr e) n) +; (liediff v (cond ((cdddr e) (cons '(mplus) (cddr e))) +; (t (caddr e))) +; n) +; ) + (cons '(mplus) (mapcar #'(lambda (u) (liediff v u n)) (cdr e))) + ) + (t (merror "~M is not a tensorial expression liediff can handle" e)) + ) +) + +(defmfun $liediff (v e) (liediff v e 1)) + +(defmfun $rediff (x) (meval '(($ev) x $idiff))) +(defmfun $evundiff (x) ($rediff ($undiff x))) + (DEFMFUN $UNDIFF (X) - (COND ((ATOM X) X) + (COND + ((ATOM X) X) ((RPOBJ X) (COND ((CDDDR X) - (NCONC (LIST '(%DERIVATIVE) + (NCONC (LIST '(%IDIFF) (LIST (CAR X) (CADR X) (CADDR X))) (PUTINONES (CDDDR X)))) (T X))) @@ -818,19 +931,19 @@ COUNT))) LOOP1(SETQ V (CADR Z)) (AND (FIXP V) - $COORDINATES + $VECT_COORDS (> V 0.) (NOT (> V $DIM)) (SETQ V - (COND ((ATOM $COORDINATES) - (MEVAL1 (LIST (LIST $COORDINATES 'SIMP 'ARRAY) + (COND ((ATOM $VECT_COORDS) + (MEVAL1 (LIST (LIST $VECT_COORDS 'SIMP 'ARRAY) V))) - ((EQ (CAAR $COORDINATES) 'MLIST) + ((EQ (CAAR $VECT_COORDS) 'MLIST) (COND ((NOT (< V - (LENGTH $COORDINATES))) + (LENGTH $VECT_COORDS))) (merror "Coordinate list too short for derivative index")) - (T (NTH V $COORDINATES)))) + (T (NTH V $VECT_COORDS)))) (T V)))) (COND ((ZEROP COUNT) (RPLACD Z (CDDDR Z)) (GO LOOP2)) ((ZEROP1 (SETQ EXP (SDIFF EXP V))) (RETURN 0.))) @@ -872,16 +985,16 @@ ((EQ (CAAR E) 'MTIMES) (ADDN (SDIFFTIMES (CDR E) X) T)) ((EQ (CAAR E) 'MEXPT) (DIFFEXPT E X)) - ((RPOBJ E) (DIFFRPOBJ E X)) ;New line added - ((AND (BOUNDP '$IMETRIC) (EQ (CAAR E) '%DETERMINANT);New line added - (EQ (CADR E) $IMETRIC)) - ((LAMBDA (DUMMY) - (setq dummy ($idummy)) - (COND ((EQ DUMMY X) (setq dummy ($idummy)))) - (LIST '(MTIMES SIMP) 2. E - (LIST '($ICHR2 SIMP) (CONS SMLIST (LIST DUMMY X)) - (CONS SMLIST (NCONS DUMMY))))) - NIL)) +;; ((RPOBJ E) (DIFFRPOBJ E X)) ;New line added +;; ((AND (BOUNDP '$IMETRIC) (EQ (CAAR E) '%DETERMINANT);New line added +;; (EQ (CADR E) $IMETRIC)) +;; ((LAMBDA (DUMMY) +;; (setq dummy ($idummy)) +;; (COND ((EQ DUMMY X) (setq dummy ($idummy)))) +;; (LIST '(MTIMES SIMP) 2. E +;; (LIST '($ICHR2 SIMP) (CONS SMLIST (LIST DUMMY X)) +;; (CONS SMLIST (NCONS DUMMY))))) +;; NIL)) ((NOT (DEPENDS E X)) (COND ((FIXP X) (LIST '(%DERIVATIVE) E X)) ((ATOM X) 0.) @@ -908,6 +1021,202 @@ ((MEMQ (CAAR E) '(%SUM %PRODUCT)) (DIFFSUMPROD E X)) (T (SDIFFGRAD E X)))) +; VTT: several of these functions have been copied verbatim from comm.lisp and +; comm2.lisp, in order to implement indicial differentiation as distinct from +; differentiation with respect to an external variable. + +(defun idiffmap (e x) (mapcar #'(lambda (term) (idiff term x)) e)) + +(defun idifftimes (l x) + (prog (term left out) + loop (setq term (car l) l (cdr l)) + (setq out (cons (muln (cons (idiff term x) (append left l)) t) out)) + (if (null l) (return out)) + (setq left (cons term left)) + (go loop))) + +(defun idiffexpt (e x) + (if (mnump (caddr e)) + (mul3 (caddr e) (power (cadr e) (addk (caddr e) -1)) (idiff (cadr e) x)) + (mul2 e (add2 (mul3 (power (cadr e) -1) (caddr e) (idiff (cadr e) x)) + (mul2 (simplifya (list '(%log) (cadr e)) t) + (idiff (caddr e) x)))))) + +(defmfun idiffint (e x) + (let (a) + (cond ((null (cdddr e)) + (cond ((alike1 x (caddr e)) (cadr e)) + ((and (not (atom (caddr e))) (atom x) (not (free (caddr e) x))) + (mul2 (cadr e) (idiff (caddr e) x))) + ((or ($constantp (setq a (idiff (cadr e) x))) + (and (atom (caddr e)) (free a (caddr e)))) + (mul2 a (caddr e))) + (t (simplifya (list '(%integrate) a (caddr e)) t)))) + ((alike1 x (caddr e)) (addn (idiffint1 (cdr e) x x) t)) + (t (addn (cons (if (equal (setq a (idiff (cadr e) x)) 0) + 0 + (simplifya (list '(%integrate) a (caddr e) + (cadddr e) (car (cddddr e))) + t)) + (idiffint1 (cdr e) x (caddr e))) + t))))) + +(defun idiffint1 (e x y) + (let ((u (idiff (cadddr e) x)) (v (idiff (caddr e) x))) + (list (if (pzerop u) 0 (mul2 u (maxima-substitute (cadddr e) y (car e)))) + (if (pzerop v) 0 (mul3 v (maxima-substitute (caddr e) y (car e)) -1))))) + +(defun idiff%deriv (e) (let (derivflag) (simplifya (cons '(%idiff) e) t))) + +(defun ideriv (e) + (prog (exp z count) + (cond ((null e) (wna-err '$idiff)) + ((null (cdr e)) (return (itotaldiff (car e)))) + ((null (cddr e)) (nconc e '(1)))) + (setq exp (car e) z (setq e (copy-top-level e))) + loop (if (or (null derivlist) (zl-member (cadr z) derivlist)) (go doit)) + ; DERIVLIST is set by $EV + (setq z (cdr z)) + loop2(cond ((cdr z) (go loop)) + ((null (cdr e)) (return exp)) + (t (go noun))) + doit (cond ((nonvarcheck (cadr z) '$idiff)) + ((null (cddr z)) (wna-err '$idiff)) + ((not (eq (ml-typep (caddr z)) 'fixnum)) (go noun)) + ((minusp (setq count (caddr z))) + (merror "Improper count to IDIFF:~%~M" count))) + loop1(cond ((zerop count) (rplacd z (cdddr z)) (go loop2)) + ((equal (setq exp (idiff exp (cadr z))) 0) (return 0))) + (setq count (f1- count)) + (go loop1) + noun (return (idiff%deriv (cons exp (cdr e)))))) + + +(defmfun idiffncexpt (e x) + ((lambda (base* pow) + (cond ((and (mnump pow) (or (not (eq (ml-typep pow) 'fixnum)) (< pow 0))) ; POW cannot be 0 + (idiff%deriv (list e x 1))) + ((and (atom base*) (eq base* x) (free pow base*)) + (mul2* pow (list '(mncexpt) base* (add2 pow -1)))) + ((ml-typep pow 'fixnum) + ((lambda (deriv ans) + (do ((i 0 (f1+ i))) ((= i pow)) + (setq ans (cons (list '(mnctimes) (list '(mncexpt) base* i) + (list '(mnctimes) deriv + (list '(mncexpt) base* (f- pow 1 i)))) + ans))) + (addn ans nil)) + (idiff base* x) nil)) + ((and (not (depends pow x)) (or (atom pow) (and (atom base*) (free pow base*)))) + ((lambda (deriv index) + (simplifya + (list '(%sum) + (list '(mnctimes) (list '(mncexpt) base* index) + (list '(mnctimes) deriv + (list '(mncexpt) base* + (list '(mplus) pow -1 (list '(mtimes) -1 index))))) + index 0 (list '(mplus) pow -1)) nil)) + (idiff base* x) (gensumindex))) + (t (idiff%deriv (list e x 1))))) + (cadr e) (caddr e))) + +(defmfun idiffsumprod (e x) + (cond ((or (not (atom x)) (not (free (cadddr e) x)) (not (free (car (cddddr e)) x))) + (idiff%deriv (list e x 1))) + ((eq (caddr e) x) 0) + (t (let ((u (idiff (cadr e) x))) + (setq u (simplifya (list '(%sum) + (if (eq (caar e) '%sum) u (div u (cadr e))) + (caddr e) (cadddr e) (car (cddddr e))) + t)) + (if (eq (caar e) '%sum) u (mul2 e u)))))) + +(defun idiffgrad (e x) + (let ((fun (caar e)) grad args) + (cond ((and (eq fun 'mqapply) (oldget (caaadr e) 'grad)) + (idiffgrad (cons (cons (caaadr e) nil) (append (cdadr e) (cddr e))) + x)) + ((or (eq fun 'mqapply) (null (setq grad (oldget fun 'grad)))) + (if (not (depends e x)) 0 (idiff%deriv (list e x 1)))) + ((not (= (length (cdr e)) (length (car grad)))) + (merror "Wrong number of arguments for ~:M" fun)) + (t (setq args (idiffmap (cdr e) x)) + (addn (mapcar + #'mul2 + (cdr (substitutel + (cdr e) (car grad) + (do ((l1 (cdr grad) (cdr l1)) + (args args (cdr args)) (l2)) + ((null l1) (cons '(mlist) (nreverse l2))) + (setq l2 (cons (cond ((equal (car args) 0) 0) + (t (car l1))) + l2))))) + args) + t))))) + +(defmfun $idiff n (let (derivlist) (ideriv (listify n)))) + +(DEFMFUN IDIFF (E X) + (COND + (($constantp E) 0.) + ((ALIKE1 E X) 1.) + ((OR (ATOM E) (MEMQ 'ARRAY (CDAR E))) +;; (ICHAINRULE E X)) +;; (idiff%deriv (list e x 1))) + 0) + ((MGET (CAAR E) '$CONSTANT) 0.) ;New line added + ((EQ (CAAR E) 'MRAT) (RATDX E X)) + ((EQ (CAAR E) 'MPLUS) + (SIMPLUS (CONS '(MPLUS) (IDIFFMAP (CDR E) X)) + 1. + T)) + ((EQ (CAAR E) 'MEQUAL) + (LIST (CAR E) ($IDIFF (CADR E) X) ($IDIFF (CADDR E) X))) + ((EQ (CAAR E) '$MATRIX) + (CONS (CAR E) + (MAPCAR + (FUNCTION (LAMBDA (Y) + (CONS (CAR Y) + (IDIFFMAP (CDR Y) X)))) + (CDR E)))) + ((EQ (CAAR E) 'MTIMES) + (ADDN (IDIFFTIMES (CDR E) X) T)) + ((EQ (CAAR E) 'MEXPT) (IDIFFEXPT E X)) + ((RPOBJ E) (DIFFRPOBJ E X)) + ((AND (BOUNDP '$IMETRIC) (EQ (CAAR E) '%DETERMINANT) + (EQ (CADR E) $IMETRIC)) + ((LAMBDA (DUMMY) + (setq dummy ($idummy)) + (COND ((EQ DUMMY X) (setq dummy ($idummy)))) + (LIST '(MTIMES SIMP) 2. E + (LIST '($ICHR2 SIMP) (CONS SMLIST (LIST DUMMY X)) + (CONS SMLIST (NCONS DUMMY))))) + NIL)) + ((EQ (CAAR E) 'MNCTIMES) + (SIMPLUS (LIST '(MPLUS) + (LIST '(MNCTIMES) + ($IDIFF (CADR E) X) + (CADDR E)) + (LIST '(MNCTIMES) + (CADR E) + ($IDIFF (CADDR E) X))) + 1. + NIL)) + ((EQ (CAAR E) 'MNCEXPT) (IDIFFNCEXPT E X)) + ((EQ (CAAR E) '%INTEGRATE) (IDIFFINT E X)) + ((EQ (CAAR E) '%DERIVATIVE) + (COND ((OR (ATOM (CADR E)) + (MEMQ 'ARRAY (CDAADR E))) +;; (ICHAINRULE E X)) +;; (idiff%deriv (list e x 1))) + 0) + ((FREEL (CDR E) X) 0.) + (T (DIFF%DERIV (LIST E X 1.))))) + ((MEMQ (CAAR E) '(%SUM %PRODUCT)) (IDIFFSUMPROD E X)) + (T (IDIFFGRAD E X)) + ) +) + (defun DIFFRPOBJ (e x) ;Derivative of an indexed object (cond ((and (memq (caar e) $COORD) (null (cdadr e)) (equal (length (cdaddr e)) 1) (null (cdddr e))) @@ -1037,7 +1346,8 @@ ((MEMQ (CAAR E) '($SUM %SUM)) (SETQ TOP (LIST (CADDR E)) BOTTOM (LIST (CADDR E))) ) - ((MEMQ (CAAR E) '(%DERIVATIVE $DIFF)) +;; ((MEMQ (CAAR E) '(%DERIVATIVE $DIFF)) + ((MEMQ (CAAR E) '(%IDIFF $IDIFF)) (DO ((I 1 (1+ I))) ((> I (COND ((CADDDR E) (CADDDR E)) (T 1)))) (SETQ BOTTOM (CONS (CADDR E) BOTTOM))) ) @@ -1160,7 +1470,7 @@ (DEFUN itensor-CLEANUP (A N)((LAMBDA (DUMX)(CLEANUP1 A)) NIL)) ;Sets DUMX to NIL (DEFUN CLEANUP1 (A) - (AND A (SETQ DUMX (IMPLODE (NCONC (EXPLODEN $DUMMYX) ;Keep proper order of + (AND A (SETQ DUMX (IMPLODE (NCONC (EXPLODEN $IDUMMYX) ;Keep proper order of (EXPLODEN N))) N (1+ N)) ;indices (COND ((EQ DUMX (CAR A)) (CLEANUP1 (CDR A))) (T (CONS (CONS (CAR A) DUMX) (CLEANUP1 (CDR A))))))) @@ -1289,14 +1599,20 @@ ((SETQ PROP (ZL-ASSOC INDEX (ZL-GET TENSOR 'TEXPRS))) ;;;VTT (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) (CADR PROP))) ;;; (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) ($RENAME (CADR PROP) (COND ((BOUNDP 'N) N) (T 1))))) - (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) (CAR (CONS ($RENAME (CADR PROP) (1+ $ICOUNTER)) (SETQ $ICOUNTER (1- (COND ((BOUNDP 'N) N) (T 1)))))))) +;; VTT: What is this business with N instead of $ICOUNTER anyway? +;;; (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) (CAR (CONS ($RENAME (CADR PROP) (1+ $ICOUNTER)) (COND ((BOUNDP 'N) (SETQ $ICOUNTER (1- N))) ))))) +;;; (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) (CAR (CONS ($RENAME (CADR PROP) (1+ $ICOUNTER)) (COND ((BOUNDP 'N) (SETQ $ICOUNTER N)) ))))) +;;; (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) (CAR (CONS ($RENAME (CADR PROP) (1+ $COUNTER)) (SETQ $COUNTER (1- (COND ((BOUNDP 'N) N) (T 1)))))))) + (SUBLIS (MAPCAR (FUNCTION CONS)(CDDR PROP) SUBS) ($RENAME (CADR PROP) (COND ((BOUNDP 'N) N) (T 1))))) + + ((SETQ PROP (ZL-GET TENSOR 'TSUBR)) ;; (APPLY PROP (LIST (CONS SMLIST (INCONSTANT L1))(CONS SMLIST (INCONSTANT L2))(CONS SMLIST L3)))) ;; ((NOT (EQ L3 NIL)) (APPLY '$DIFF (SELECT TENSOR (INCONSTANT L1) (INCONSTANT L2) (CDR L3)) (LIST (CAR L3)))) ;; (T (APPEND (LIST (LIST TENSOR 'SIMP)(CONS SMLIST (INCONSTANT L1))(CONS SMLIST (INCONSTANT L2))) L3)))) ;; NIL (APPEND (INCONSTANT L1) (INCONSTANT L2) L3)(LIST (LENGTH (INCONSTANT L1))(LENGTH (INCONSTANT L2))(LENGTH L3)))) (APPLY PROP (LIST (CONS SMLIST L1)(CONS SMLIST L2)(CONS SMLIST L3)))) - ((NOT (EQ L3 NIL)) (APPLY '$DIFF (SELECT TENSOR L1 L2 (CDR L3)) (LIST (CAR L3)))) + ((NOT (EQ L3 NIL)) (APPLY '$IDIFF (SELECT TENSOR L1 L2 (CDR L3)) (LIST (CAR L3)))) (T (APPEND (LIST (LIST TENSOR 'SIMP)(CONS SMLIST L1)(CONS SMLIST L2)) L3)))) NIL (APPEND L1 L2 L3)(LIST (LENGTH L1)(LENGTH L2)(LENGTH L3)))) Index: kaluza.dem =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/kaluza.dem,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- kaluza.dem 18 Nov 2004 03:53:28 -0000 1.4 +++ kaluza.dem 25 Nov 2004 01:43:29 -0000 1.5 @@ -1,41 +1,105 @@ -("Deriving the Kaluza-Klein equation of motion.")$ -("For reference, see http://www.vttoth.com/KK/kk.htm")$ +("Deriving the Kaluza-Klein equation of motion. +For reference, see http://www.vttoth.com/KK/kk.htm")$ + derivabbrev:true; -("KALUZA.MAC contains definitions and contraction properties for G4 and G5")$ -load("kaluza.mac"); -("Predeclaring some 4d indices")$ +("We first load ITENSOR and set up the 5-dimensional metric. +We also set up contraction properties for both the 4-dimensional +and the 5-dimensional metric tensors.")$ + +load("itensor"); +dim:5; +imetric:g5; +defcon(g4); +defcon(g5); +defcon(g4,g4,kdelta); +defcon(g5,g5,kdelta); + +("To set up the metric components, we need some helper functions. +The function predval() determines if a predicate can be evaluated. +It returns false if the predicate would return an error. The +function difflist() applies the differential operator to elements +in a list.")$ + +predval(prd):=block([retval,saved_prederror:prederror], + prederror:false, + retval:ev(prd,pred)=true or ev(prd,pred)=false, + prederror:saved_prederror, + retval +)$ +difflist(exp,lst):=if length(lst)=0 then exp + else difflist(diff(exp,lst[1]),rest(lst))$ + +("Metric components are defined conditionally, allowing us to treat +the fifth index in a unique way.")$ + +a(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('a,append([l1,l2],l3))$ +g4(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('g4,append([l1,l2],l3))$ +g5(l1,l2,[l3]):= + if member(5,l3) then 0 + else if l1#[] then + ( + if not (predval(l1[1]<=4) and predval(l1[2]<=4)) then + funmake('g5,append([l1,l2],l3)) + else if l1[1]<=4 and l1[2]<=4 then + apply('g4,append([l1,l2],l3))+ + g55*difflist(a([l1[1]],[])*a([l1[2]],[]),l3) + else if l1[1]<=4 then g55*apply('a,append([[l1[1]],[]],l3)) + else if l1[2]<=4 then g55*apply('a,append([[l1[2]],[]],l3)) + else if l3#[] then 0 else g55 + ) + else if l2#[] then + ( + if not (predval(l2[1]<=4) and predval(l2[2]<=4)) then + funmake('g5,append([l1,l2],l3)) + else if l2[1]<=4 and l2[2]<=4 then apply('g4,append([l1,l2],l3)) + else if l2[1]<=4 then -apply('a,append([[],[l2[1]]],l3)) + else if l2[2]<=4 then -apply('a,append([[],[l2[2]]],l3)) + else if l3#[] then sum(difflist(a([i],[])*a([],[i]),l3),i,1,4) + else 1/g55+sum(a([i],[])*a([],[i]),i,1,4) + ) + else funmake('g5,append([l1,l2],l3))$ + +("We also need some simplification rules:")$ +matchdeclare([dummy1,dummy2],true); +defrule(evpot,a([dummy1],[],dummy2)-a([dummy2],[],dummy1), + -f([dummy1,dummy2],[]))$ + +("Now we're ready to begin the analysis. First, we predeclare +some 4-dimensional indices:")$ assume(k<=4,l<=4,m<=4)$ -("Equation of motion in empty 5-space")$ -ishow('diff(x([],[a]),t,2)+'ichr2([b,c],[a])*'diff(x([],[b]),t)*'diff(x([],[c]),t)=0)$ +("The equation of motion in empty 5-space:")$ +depends(x,t); +ishow('diff(x([],[a]),t,2)+ + 'ichr2([b,c],[a])*'diff(x([],[b]),t)*'diff(x([],[c]),t)=0)$ ishow(part(first(%),1))$ ishow(subst(m,c,%)+subst(5,c,%))$ ishow(subst(l,b,%)+subst(5,b,%)+part(first(%th(3)),2)=last(%th(3)))$ -("We are only interested in the case where A is a 4D index")$ +("We are only interested in the case where A is a 4D index:")$ ishow(subst(k,a,%th(2)))$ -("We protect one of the Christoffel-symbols from expansion")$ +("We protect one of the Christoffel-symbols from expansion:")$ ishow(subst(chr2klm,'ichr2([l,m],[k]),%th(2)))$ %,ichr2$ ishow(rename(%))$ -("Now we break this up into two parts depending on whether %1=5")$ -map(lambda([u],block(if freeof(%1,u) then u else u+subst(5,%1,u))),first(%th(2)))=last(%th(2))$ +("Now we break this up into two parts depending on whether %1=5:")$ +map(lambda([u],block(if freeof(%1,u) then u else u+subst(5,%1,u))), + first(%th(2)))=last(%th(2))$ assume(%1<=4)$ %th(2),g5$ %,nouns$ ishow(%)$ -("Now we're ready to isolate the electromagnetic field tensor")$ +("Now we're ready to isolate the electromagnetic field tensor:")$ map(lambda([u],factorout(u,g55)),%th(2))$ ishow(apply1(%,evpot))$ -("Contracting and rearranging yields the equation in the usual form")$ +("Contracting and rearranging yields the equation in the usual form:")$ contract(%th(2))$ %,nouns$ ishow(rename(%))$ - %-part(first(%),1)$ ishow(subst('ichr2([l,m],[k]),chr2klm,%))$ @@ -49,12 +113,13 @@ assume(%1<=4)$ %th(6),g5$ -scanmap(lambda([u],apply1( ratsimp(u,g55,g4([],[%1,m]),a([k],[])) ,evpot)),%)$ +scanmap(lambda([u],apply1(ratsimp(u,g55,g4([],[%1,m]),a([k],[])),evpot)),%)$ ishow(%)$ %+%th(5)$ %,nouns$ -ratsubst(chr42([k,l],[m]),g4([],[%1,m])*(g4([l,%1],[],k)+g4([k,%1],[],l)-g4([k,l],[],%1))/2,%)$ +ratsubst(chr42([k,l],[m]), + g4([],[%1,m])*(g4([l,%1],[],k)+g4([k,%1],[],l)-g4([k,l],[],%1))/2,%)$ ishow(%)$ contract(%)$ Index: manual.txt =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/manual.txt,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- manual.txt 18 Nov 2004 03:53:28 -0000 1.3 +++ manual.txt 25 Nov 2004 01:43:29 -0000 1.4 @@ -357,9 +357,9 @@ .end .endfunction -.function(MOTION,dis) +.function(CGEODESIC,dis) computes the covariant form of the geodesic equations of motion for a -given metric. They are stored in the array EM[i]. If the argument +given metric. They are stored in the array GEOD[i]. If the argument 2dis* is TRUE then these equations are displayed. .endfunction @@ -406,7 +406,7 @@ .function(WEYL,dis) computes the covariant Weyl conformal tensor. If the argument 2dis* -is TRUE, the non-zero components W[i,j,k,l] will be displayed. +is TRUE, the non-zero components WEYL[i,j,k,l] will be displayed. Otherwise, these components will be computed and stored. If the switch $var<RATWEYL/TRUE> is set to TRUE, then the components will be rationally simplified. If RATFAC is TRUE then the results will be @@ -887,12 +887,12 @@ .skip 2) the 2vi* may be integers from 1 up to the value of the variable DIM[4]. This will cause the differentiation to be carried out with -respect to the 2vi*th member of the list COORDINATES which should be +respect to the 2vi*th member of the list VECT_COORDS which should be set to a list of the names of the coordinates, e.g., [x,y,z,t] . If -COORDINATES is bound to an atomic variable, then that variable +VECT_COORDS is bound to an atomic variable, then that variable subscripted by 2vi* will be used for the variable of differentiation. This permits an array of coordinate names or -subscripted names like X[1], X[2],... to be used. If COORDINATES has +subscripted names like X[1], X[2],... to be used. If VECT_COORDS has not been assigned a value, then the variables will be treated as in 1) above. .break continue @@ -1570,7 +1570,7 @@ already in use (see the example under INDICES). .endfunction -.function2 DUMMYX +.function2 IDUMMYX is the prefix for dummy indices (see the example under INDICES). .endfunction @@ -1621,12 +1621,12 @@ .end .skip .begin group -(C8) DUMMYX; +(C8) IDUMMYX; (D8) % .end .skip .begin group -(C9) DUMMYX:&$ +(C9) IDUMMYX:&$ .end .skip .begin group Index: swartz.dem =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/swartz.dem,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- swartz.dem 18 Nov 2004 03:53:28 -0000 1.5 +++ swartz.dem 25 Nov 2004 01:43:29 -0000 1.6 @@ -68,19 +68,19 @@ ("computes scalar curvature")$ scurvature(); ("computes Riemann tensor")$ -riemann(true)$ +lriemann(true)$ ("computes contravariant Riemann tensor")$ uriemann(false)$ ("computes the Kretchmann invariant Rijkl^2")$ rinvariant(); diagmetric:true; ("Compute the covariant form of geodesic equations")$ -motion(true)$ +cgeodesic(true)$ -("Compute the contravariant form geodesic equations"); +/*("Compute the contravariant form geodesic equations"); block( for i thru dim do - emc[i]: factorfacsum(ratexpand(sum(ug[i,a]*em[a],a,1,dim))), - for i thru dim do ldisplay(emc[i])); + emc[i]: factorfacsum(ratexpand(sum(ug[i,a]*geod[a],a,1,dim))), + for i thru dim do ldisplay(emc[i]));*/ block([title: "Schwarzschild Potential for Mass M=2",m:2.], Index: symtry.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/symtry.lisp,v retrieving revision 1.11 retrieving revision 1.12 diff -u -d -r1.11 -r1.12 --- symtry.lisp 18 Nov 2004 03:53:28 -0000 1.11 +++ symtry.lisp 25 Nov 2004 01:43:29 -0000 1.12 @@ -4,7 +4,7 @@ ; ** (c) Copyright 1979 Massachusetts Institute of Technology ** -(declare-top (special symtypes $symmetries $allsym csign smlist $dummyx)) +(declare-top (special symtypes $symmetries $allsym csign smlist $idummyx)) (setq symtypes '($SYM $ANTI $CYC) $symmetries '((MLIST SIMP))) @@ -124,7 +124,7 @@ (defun CANTEN (e nfprpobjs) ;CANonical TENsor (prog (cov contr deriv tensor) ((lambda (dummy) (and nfprpobjs dummy (setq e (rename1 e dummy)))) - (NONUMBER (cdaddr ($indices2 e)))) ;NFPRPOBJS is Not From Product + (NONUMBER (cdaddr ($indices e)))) ;NFPRPOBJS is Not From Product (setq cov (copy (cdadr e)) ;of RP (indexed) OBJects contr (copy (cdaddr e)) deriv (copy (cdddr e)) @@ -165,7 +165,7 @@ (defun CLEANUP0 (a) (do ((b a (cdr b)) (n 1 (1+ n)) (l) (dumx)) ((null b) l) - (setq dumx (concat $dummyx n)) + (setq dumx (concat $idummyx n)) (cond ((not (eq dumx (car b))) (setq l (cons (cons (car b) dumx) l)))))) @@ -241,13 +241,13 @@ ((lambda (q) (rename1 q (NONUMBER (cdaddr - ($indices2 + ($indices (cons '(MTIMES) q)))))) (mapcar 'cdr (sortcar (progn (setq free-indices - (NONUMBER (cdadr ($indices2 e)))) + (NONUMBER (cdadr ($indices e)))) (mapcar 'describe-tensor indexed)) 'tensorpred)))))))))) Index: ten_alg.dem =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/ten_alg.dem,v retrieving revision 1.8 retrieving revision 1.9 diff -u -d -r1.8 -r1.9 --- ten_alg.dem 18 Nov 2004 03:53:28 -0000 1.8 +++ ten_alg.dem 25 Nov 2004 01:43:29 -0000 1.9 @@ -28,8 +28,8 @@ g([],[l2[1],j1])*(diff(p([],[]),j1)/rho([])+diff(phi([],[]),j1)) )$ ishow(Eu([i],[]))$ ishow(Eu([],[i]))$ -("control the dummy indices with dummyx: something")$ -dummyx:n$ +("control the dummy indices with idummyx: something")$ +idummyx:n$ ishow(Eu([i],[]))$ @@ -107,17 +107,17 @@ ishow(canform(contract(%)))$ ("the exterior algebra for the antisymmetric covariant objects")$ -("the exterior product is defined through the infix operator & - the exterior derivative through extdiff(), inner product with |_ vec, - and the Lie derivative through liediff(x, [vec, ind])")$ +("the exterior product is defined through the infix operator ~ + the exterior derivative through extdiff(), inner product with | vec, + and the Lie derivative through liediff(vec,x)")$ ("load the package in the system")$ ("now look at the exterior product of two 1-forms and 1-form with 2-form")$ allsym:false$ decsym(p,2,0,[anti(all)],[])$ kill(a)$ -ishow(a([i])&b([j]))$ -ishow(a([i])&p([j,k])+a([i])&b([j])&a([k]))$ +ishow(a([i])~b([j]))$ +ishow(a([i])~p([j,k])+a([i])~b([j])~a([k]))$ ("dualization")$ ("the itensor's old-style Levi-Civita symbol operates only on integers")$ @@ -174,10 +174,10 @@ ("exterior derivatives")$ (allsym:false,decsym(p,2,0,[anti(all)],[]),ishow(extdiff(p([i,j]),k)))$ ("the Lie derivative of a 1-form")$ -ishow(liediff(a([i]),[v,k]))$ +ishow(liediff(v,a([i])))$ ("the Lie derivative has to commute with the exterior one.")$ -ishow(extdiff(liediff(a([i]),[v,j]),k)-liediff(extdiff(a([i]),j),[v,k]))$ -ishow(extdiff(liediff(p([i1,i2]),[v,i3]),i4)-liediff(extdiff(p([i1,i2]),i3),[v,i4]))$ +ishow(canform(extdiff(liediff(v,a([i])),j)-liediff(v,extdiff(a([i]),j))))$ +ishow(canform(extdiff(liediff(v,p([i1,i2])),i3)-liediff(v,extdiff(p([i1,i2]),i3))))$ ("more examples about the Cartan identity are in car_iden.dem")$ @@ -244,7 +244,7 @@ ishow(canform(exp))$ ("transition to the coordinates")$ -dummyx:di$ +idummyx:di$ remcomps(g); imetric:g$ dim:3$ @@ -269,14 +269,14 @@ kill(ct_coords)$ ("the conservation of the helicity in hydrodynamics")$ (if allsym=true then allsym:false)$ -dummyx:j$ +idummyx:j$ ("the covariant formulation of the Euler equation is")$ euler: d_t*v([i1])+ -'liediff(v([i0]),[v,i1])=-extdiff(p([]),i1)/rho([])+extdiff((v([i0])|_v/2-phi([])),i1)$ +'liediff(v,v([i0]))=-extdiff(p([]),i1)/rho([])+extdiff((v([i0])|v/2-phi([])),i1)$ ishow(euler)$ ("take the exterior derivative of both sides")$ -euler1:d_t*ct_coords([i1,i2])+'liediff(ct_coords([i0,i1]),[v,i2])= -extdiff(-extdiff(p([]),i1)/rho([])+extdiff(v([i0])|_v/2-phi([]),i1),i2)$ +euler1:d_t*ct_coords([i1,i2])+'liediff(v,ct_coords([i0,i1]))= +extdiff(-extdiff(p([]),i1)/rho([])+extdiff(v([i0])|v/2-phi([]),i1),i2)$ ishow(euler1)$ ishow(ct_coords([i1,i2])=extdiff(v([i1]),i2))$ ("declare the function's dependencies")$ @@ -294,7 +294,7 @@ ishow(euler2)$ ("so the rhs will vanish after forming the exterior product with the differential of the entropy")$ -rhs(euler1)&extdiff(s([]),i3)$ +rhs(euler1)~extdiff(s([]),i3)$ ishow(factor(expand(apply1(%,p_sr))))$ ("the magnetic field lines are embedded in the moving plasma")$ @@ -302,17 +302,17 @@ ("introduce the 2-form describing the magnetic field")$ decsym(b,2,0,[anti(all)],[])$ ("due to the infinite conductivity the electric field is expressed through")$ -e_field:b([i1,i2])|_v$ +e_field:b([i1,i2])|v$ equat: d_t*b([i2,i3])=-extdiff(e_field,i3)$ ishow(equat)$ ("express B through the vector potential")$ ("i.e., re-compute the electric field")$ -e_field1:extdiff(a([i1]),i2)|_v$ +e_field1:extdiff(a([i1]),i2)|v$ ("substitute it in Equat")$ equat1: d_t*b([i2,i3])=-extdiff(e_field1,i3)$ ("compute the Lie derivative of B")$ -lie_derb: 'liediff(b([i1,i2]),[v,i3])=extdiff(extdiff(a([i1]),i2),L[v,i3])$ +lie_derb: 'liediff(v,b([i1,i2]))=liediff(v,extdiff(a([i1]),i2))$ ("summing up the lhs and rhs parts")$ emded:lhs(equat1)+lhs(lie_derb)=canform(rhs(equat1)+rhs(lie_derb))$ ishow(emded)$ Index: tensor-doc.txt =================================================================== RCS file: /cvsroot/maxima/maxima/share/tensor/tensor-doc.txt,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- tensor-doc.txt 18 Nov 2004 03:53:28 -0000 1.5 +++ tensor-doc.txt 25 Nov 2004 01:43:29 -0000 1.6 @@ -401,8 +401,8 @@ 4 A D X -MOTION(dis) computes the covariant form of the geodesic equations of - motion for a given metric. They are stored in the array EM[i]. If +CGEODESIC(dis) computes the covariant form of the geodesic equations of + motion for a given metric. They are stored in the array GEOD[i]. If the argument dis is TRUE then these equations are displayed. @@ -446,7 +446,7 @@ WEYL(dis) computes the covariant Weyl conformal tensor. If the argument - dis is TRUE, the non-zero components W[i,j,k,l] will be displayed. + dis is TRUE, the non-zero components WEYL[i,j,k,l] will be displayed. Otherwise, these components will be computed and stored. If the switch RATWEYL[TRUE] is set to TRUE, then the components will be rationally simplified. If RATFAC is TRUE then the results will be @@ -844,12 +844,12 @@ 2) the vi may be integers from 1 up to the value of the variable DIM[4]. This will cause the differentiation to be carried out with - respect to the vith member of the list COORDINATES which should be + respect to the vith member of the list VECT_COORDS which should be set to a list of the names of the coordinates, e.g., [x,y,z,t] . If - COORDINATES is bound to an atomic variable, then that variable + VECT_COORDS is bound to an atomic variable, then that variable subscripted by vi will be used for the variable of differentiation. This permits an array of coordinate names or subscripted names like - X[1], X[2],... to be used. If COORDINATES has not been assigned a + X[1], X[2],... to be used. If VECT_COORDS has not been assigned a value, then the variables will be treated as in 1) above. @@ -1410,7 +1410,7 @@ with indices already in use (see the example under INDICES). -DUMMYX is the prefix for dummy indices (see the example under INDICES). +IDUMMYX is the prefix for dummy indices (see the example under INDICES). INDICES(exp) returns a list of two elements. The first is a list of the @@ -1443,10 +1443,10 @@ (- ICHR2 - ICHR2 ICHR2 + ICHR2 + ICHR2 ICHR2 ) I K,J %12 J I K I J,K %12 K I J -(C8) DUMMYX; +(C8) IDUMMYX... [truncated message content] |