Update: Calculate Tensor products and traces with Mathematica

Question

My goal is to algebraically simplify expressions involving tensor products and taking the trace.

That is, I would like to compute $\operatorname{Tr}\left( \prod_i(id+a_i \otimes b_i) \right)$,

using the rules of tensor algebra. E.g. $$\operatorname{Tr}(a_i \otimes b_i)=a_i.b_i$$ and $$a_1 \otimes b_1 a_2 \otimes b_2= b_1.a_2 a_1 \otimes b_2.$$

The final expression should be of the form

$$ \sum_i \prod_j a_{i_j}.b_{i_j}.$$

Any implementation to let Mathematica do this for me would be much appreciated.

So far I tried to use the built-in tensor algebra package

$Assumptions = {Element[a1, Vectors[3, Reals]], 
   Element[a2, Vectors[3, Reals]], Element[b1, Vectors[3, Reals]], 
   Element[b2, Vectors[3, Reals]], Element[n, Vectors[3, Reals]], 
   Element[m, Vectors[3, Reals]], Element[id, Matrices[{3, 3}]]};

and define simplification rules such as

SimId[expr_] := 
 expr //. {Dot[id, tensor___] :> tensor, Dot[tensor___, id] :> tensor,
    MatrixPower[id, 2] :> id, MatrixPower[id, 3] :> id, 
   MatrixPower[id, 4] :> id}
SimTP[expr_] := 
 expr //. {Dot[TensorProduct[before1_, before2_], 
     TensorProduct[after1_, after2_]] :> 
    Dot[before2, after1] TensorProduct[before1, after2]}

This works fine for simple expressions such as

    SimId[SimTP[
  TensorProduct[a2 + b1, 
     n1].(id + TensorProduct[a1, a1 + b2]).TensorProduct[b1, b2] // 
   TensorExpand]]

However, once I enter more complex expressions, it seems to go wrong and give me

    SimId[SimTP[
  TensorProduct[a2, n1].TensorProduct[b2, n1].TensorProduct[a2, 
     b1].TensorProduct[a2, b2] // TensorExpand]]

I am still getting invalid expression of the form

a2.b1 (n1.b2 a2\[TensorProduct]n1).a2\[TensorProduct]b2

and they don't really make sense.

Answer 1

You cannot use @ as it is used for assignments etc. So just set up your own tensor product tp and inner product ip and define

tp[id,v_]:=v
tp[v_,id]:=v
Tra[tp[a,b]]:=ip[a,b]
Trans[tp[a,b]]:=tp[b,a]
ip[tp[a,b],tp[c,d]]:=ip[b,c] tp[a,d]

and

tp[a_+b_,c_]:=tp[a,c]+tp[b,c]

It is unclear as to what you mean by "respect scalar multiplication". If you only want this with respect to explicit numbers you can set, for example,

tp[(s : _Integer | _Rational | _Real|_Complex) a_,b_]:=s tp[a,b]

If, however, you want to use variables that are scalar, you first need to tell Mathematica which symbols represent vectors or scalars. For example, you can define

myVectorsQ[exp_]:=MatchQ[exp,a|b|c]

which is supposed to mean that the symbols a, b and c are vectors and everything else a scalar. Now, you can set

tp[s_*a_?myVectorsQ,b_]:=s tp[a,b]
tp[a_,s_*b_?myVectorsQ]:=s tp[a,b]

Note that this approach is independent of the dimension and that one can define other functions on those tensors (like derivatives or integrals).

Answer 2

I managed to to it relatively simply. This introduces the objects:

$Assumptions = {Element[x, Matrices[{3, 1}, Reals]], Element[a, Matrices[{3, 1}, Reals]], Element[n, Matrices[{3, 1}, Reals]], Element[b, Matrices[{3, 1}, Reals]], Element[m, Matrices[{3, 1}, Reals]], Element[id, Matrices[{3, 3}, Reals, Symmetric]]};

and this set of simplification rules does what I wanted it to do.

Simp[expr_] := expr //. {x__.Dot[Transpose[tensor1__], tensor2__].y__ :> tensor1.tensor2*x.y} SimId[expr_] := expr //. {Dot[id, tensor_] :> tensor, Dot[tensor_, id] :> tensor, MatrixPower[id, 2] :> id, MatrixPower[id, 3] :> id, MatrixPower[id, 4] :> id, MatrixPower[id, 5] :> id, MatrixPower[id, 6] :> id, MatrixPower[id, 7] :> id, MatrixPower[id, 8] :> id, MatrixPower[id, 9] :> id} SimTranspose[expr_] := expr //. {Transpose[n_, {2, 1}] :> Transpose[n], Dot[x_, {2, 1}] :> x} SimTr[expr_] := expr //. {Dot[tensor2__, Transpose[tensor1__]] :> tensor2.tensor1} Sim[expr_] := expr //. {Dot[a, n] :> 0, Dot[n, a] :> 0, Dot[b, m] :> 0, Dot[m, b] :> 0, Dot[n, n] :> 1, Dot[m, m] :> 1} SimAll[expr_] := Sim[SimTr[SimTranspose[Simp[SimId[TensorExpand[expr]]]]]]

Then one can just enter any expression of the described from and get the Trace. The function Sim is optional and can contain certain things things that cancel. E.g. $a \perp n \Rightarrow a.n=0$.