To retrieve the eight elements around a given entry in some matrix, using torus topology, so to speak, the following function might be useful:

getEntries := proc(M::Matrix,x::posint,y::posint)
   local
      rowDim := LinearAlgebra:-RowDimension(M),
      colDim := LinearAlgebra:-ColumnDimension(M),
      rowMod := (x) -> ((x-1) mod rowDim) + 1,
      colMod := (y) -> ((y-1) mod colDim) + 1;
      [
         M[rowMod(x  ),colMod(y+1)],
         M[rowMod(x-1),colMod(y+1)],
         M[rowMod(x-1),colMod(y  )],
         M[rowMod(x-1),colMod(y-1)],
         M[rowMod(x  ),colMod(y-1)],
         M[rowMod(x+1),colMod(y-1)],
         M[rowMod(x+1),colMod(y  )],
         M[rowMod(x+1),colMod(y+1)]
      ]
end proc:

The eight entries returned are ordered going counter-clockwise around the 'center'-entry. An example of use:

M := Matrix(3,5,symbol = 'm');
getEntries(M,2,3);
getEntries(M,1,1);

It is a good question. The problem seems to boil down to the wedge product of two anticommuting scalars being erroneously equal irrespective of their ordering:

alpha &wedge beta,
beta &wedge alpha;
alpha*beta,
beta*alpha;

Looking at the code of &wedge, using

showstat(`&wedge`);

this must, it seems, be related to line 8 (at least in Maple 17, which I am using): return FB*FA. Related to that there are two questions:

- Why actually the reverse ordering FB*FA, rather than FA*FB?

- Is Physics:-* active or not (I guess it should be, the Physics package having been loaded), i.e., is the multiplication being performed using the more general multiplication of the Physics package or not? Surprisingly, it seems not to be.

Perhaps the Physics package would be useful: First let

with(Physics):
Setup(noncommutativeprefix = {X,Y,Z});

The Setup line defines any variables that starts with X, Y, or Z to be noncommutative, so that Y*X, say, is not simplified to X*Y, but kept as it is. The expression you mention can then be given as

Commutator(Commutator(a*X+b*Y,c*Z),Commutator(c*Y,b*X+c*Z));

or expanded and simplified, if needed, as

simplify(expand(Commutator(Commutator(a*X+b*Y,c*Z),Commutator(c*Y,b*X+c*Z))));

PS: As the following calculation shows, the Jacobi identity is, unsurprisingly, automatically taken care of:

expand(
   +Commutator(Commutator(X,Y),Z)
   +Commutator(Commutator(Y,Z),X)
   +Commutator(Commutator(Z,X),Y)
);

0

From your finite summation ansatz, it seems to me that it would be useful to have a function that could expand sin(k*x/2) and cos(k*x/2), for any positive integer k, entirely in terms of only (products of) sin(x/2) and cos(x/2) raised to various powers. For in that way it might be possible to pair of terms with identical powers, and thus find the coefficients of the sum.

As Maple does not readily expand even sin(3/2*x), say, I have written the following function which recursively, lowering k by 1/2 in each step, performs the above-mentioned full expansion:

FullExpand := proc(t::algebraic)
   local a,b,r,T := expand(t);
   if T::{`+`,`*`,`^`} then expand(map(FullExpand,T))
   elif T::Or(
      identical(sin(x)),
      identical(cos(x)),
      specfunc(`&*`(realcons,identical(x)),sin),
      specfunc(`&*`(realcons,identical(x)),cos)
   ) then
      r := `if`(op(1,T)::identical(x),1,select(type,op(T),realcons));
      `if`(r > 1/2,FullExpand(eval(expand(op(0,T)(a*x + b*x)),{
         a = r - 1/2,
         b =   + 1/2
      })),T)
   else T
   end if;
end proc:

An example:

FullExpand(sin(5/2*x)^2 + 2*cos(3*x));

Using the Physics package, the summation over i and j can be performed automatically as follows:

with(Physics):
Setup(
   mathematicalnotation = true,
   dimension = 3,
   spacetimeindices = lowercaselatin
):
Define(A,B):
Define(Sigma[k,l] = A[k,i]*A[l,j]*B[~i,~j]):

Note that the last line defines Sigma for any k,l. Just execute

Sigma[1,1];

for the component that you want.

Using Maple 17, I get the following:

dsolve(deqv);

,

If there is no difference between superscripts and subscripts, then perhaps you could do as follows:

Theta := (m,theta) -> piecewise(m >= 0,cos(m*theta),-sin(m*theta)):
Z := (m,n) -> (theta,rho) -> N[n,m]*R[m,n](rho)*Theta(m,theta):

Examples of use:

Z(+2,3)(theta,rho);
Z(-2,3)(theta,rho);

PS: You would still have to specify, I guess, the quantities N and R, the latter being functions of rho.

The two solutions by Mathematica are both included in the solution by Maple, the first by the choice _C1 = +1/2*C[1], the second by the choice _C1 = -1/2*C[1]. Remember that _C1 may be an arbitrary constant.

I don't know if it counts for anything, but when I want to do an analogous thing in a large document, not necessarily in Maple, I just type 'xxx' at the two places that I want to toggle between, and then simply use Ctrl+f to search for them (jumping first to the top of the document using Ctrl+Home if searching for the first one).

It is perhaps not that pretty, and it cannot be used to extract the constant factor of the polynomial, but maybe the following is useful:

getCoeff := (f::algebraic,inds::set,expr::algebraic) ->
   remove(has,coeff(algsubs(expr = 'x',expand(f)),'x'),inds):

Some examples of its use:

getCoeff(f,{a,b,c},a^2*b^2*c^2),
getCoeff(f,{a,b,c},a^2*b^2),
getCoeff(f,{a,b,c},a*b);

It seems that igcd returns zero if and only if all its arguments are zero. For your code, this leads to 61 instances of division by zero:

n := 0:
for a from -10 to 10 do
for b from -10 to 10 do
for c from -10 to 10 do
   if igcd(a*b,b*c,c*a,a*b*c) = 0 then n := n + 1 end if
end do
end do
end do;
n;

check := (M::Matrix) -> evalb(nops((`+`@op)~({
   convert(M   ,listlist)[],
   convert(M^%T,listlist)[]
})) = 1):
check(A);

Perhaps the following is useful:

createDegree := proc(L::list,n::posint)
   local powers;
   powers := foldl(
      (x,y) -> [seq(seq([i[],j[]],i in x),j in y)],
      map(x -> `[]`~([$0..floor(n/degree(x))]),L)[]
   );
   select(x -> 0 < degree(x) and degree(x) <= n,
      map(p -> `*`(zip((x,y) -> x^y,L,p)[]),powers)
   )
end proc:

A suggestion:

createDegree := (L::{list,set},n::posint) -> select(x -> degree(x) <= n,L):

I would do something like the following:

with(Physics):
Setup(mathematicalnotation = true,dimension = 3):
Coordinates(X = (r,theta,phi)):
Setup(metric = Matrix(3,3,(i,j) -> g||i||j(X),shape = symmetric)):
g_[];

where I have given some arbitrary metric. Note that the coordinates are defined outside the Setup environment. Note also that in the realm of tensor calculus of Riemannian geometry, there is no need to specify explicitly that some coordinate system is spherical, say; the coordinates and the metric in conjunction is all there is. But I guess you already know that.

I have chosen to use phi instead of varphi, because the latter does not render well visually.