# How can I automate this tensor computation?

I am doing this work by hand, but it takes a lot of time and I make several calculation errors, so I was thinking to make Mathematica to calculate this for me, but I am stuck at the very beginning.

I am working with tensors like this:

$$XXV_{ijk} = \dfrac{1}{6}(X_iX_jV_k+X_iX_kV_j+X_jX_iV_k+X_kX_jV_i+X_jX_kV_i+X_kX_iV_j)$$ $$-\dfrac{1}{5}(\delta_{ij}\, (X\cdot X) V_k+\delta_{ik} (X\cdot V) X_j+\delta_{jk} (X\cdot V) X_i)$$

where $$X_i$$ and $$V_i$$ are the components of the 3-vectors $$\vec{X}$$ and $$\vec{V}$$. I have to multiply these tensors, for example

$$XXV \times XXV = \dfrac{2}{25}V^2+\dfrac{8}{25} (X\cdot V)^2$$

I think I can obtain the first part with Tuples and Total(?) but I don't know how to obtain the part with the Kroeneker deltas; if I can write these tensors correctly I think I can multiply these tensors with . and Transpose.

As @yarchik has pointd out, I have to add that my tensors have unit length

• Is there a typo in your last term on the first line, $X_kX_iV_k$? There should be a $j$-index somewhere. – Roman Jul 19 at 14:52
• @Roman yes it is a typo – mattiav27 Jul 19 at 15:13
• I think your second term still needs a power correction. The answer should be proportional to $V^2$. – yarchik Jul 19 at 15:40
• @yarchik you are right, but this is a typo not a calculation error... – mattiav27 Jul 19 at 15:44
• Note that in your example the the correct result should be $\frac{8}{25} V^2 + \frac{2}{25} (X\cdot V)^2$. See all four answers below. – Shadowray Jul 20 at 0:15

You can write it directly as you see it

xxv[i_,j_,k_]:= 1/6( x[i]x[j]v[k]+x[i]x[k]v[j]
+x[j]x[i]v[k]+x[k]x[j]v[i]
+x[j]x[k]v[i]+x[k]x[i]v[j] )
-1/5( KroneckerDelta[i,j]Sum[x[l]x[l],{l,3}]v[k]
+KroneckerDelta[i,k]Sum[x[l]v[l],{l,3}]x[j]
+KroneckerDelta[j,k]Sum[x[l]v[l],{l,3}]x[i] )

FullSimplify[ Sum[xxv[i,j,k] xxv[i,j,k],{i,3},{j,3},{k,3}],
Assumptions->Sum[x[i]^2,{i,3}]==1
&&Sum[x[i]v[i],{i,3}]==xv
&&Sum[v[i]v[i],{i,3}]==vv]

Out[1]= 2/25 (4 vv + xv^2)


where I assumed that your vector x is normalized

You can cast this as a symbolic tensor question, and make use of my TensorSimplify package. Install the paclet with:

PacletInstall[
"TensorSimplify",
"Site" -> "http://raw.githubusercontent.com/carlwoll/TensorSimplify/master"
]


Once installed, load the package with:

<<TensorSimplify


Now, define your tensor using TensorProduct:

XXV = 1/3 (TensorProduct[X,X,V] + TensorProduct[X,V,X] + TensorProduct[V,X,X]) -
1/5 (X.X TensorProduct[Inactive[IdentityMatrix][3], V] +
X.V TensorTranspose[TensorProduct[Inactive[IdentityMatrix][3],X],{1,3,2}] +
X.V TensorProduct[X,Inactive[IdentityMatrix][3]]
);


Note the use of Inactive[IdentityMatrix][3] instead of IdentityMatrix[3]. Then:

TensorSimplify[
TensorContract[TensorProduct[XXV, XXV], {{1, 4}, {2, 5}, {3, 6}}],
Assumptions -> (X|V) ∈ Vectors[3]
]


2/25 (V.X)^2 X.X + 8/25 V.V (X.X)^2

Using X.X == 1 reproduces your result.

This is how I'd do it; maybe it's useful for you.

Define $$\vec{X}$$ and $$\vec{V}$$ as vectors:

X = Array[x, 3];
V = Array[v, 3];


useful $$3\times3\times3$$ tensors for assembling:

a = Outer[Times, X, X, V];
b = (X.X) Outer[Times, IdentityMatrix[3], V];
c = (X.V) Outer[Times, IdentityMatrix[3], X];


assemble $$XXV$$:

XXV = (a + Transpose[a, {3, 1, 2}] + Transpose[a, {2, 3, 1}])/3 -
(b + Transpose[c, {3, 1, 2}] + Transpose[c, {2, 3, 1}])/5;


check a formula:

Total[XXV*XXV, 3] == 2/25 (X.X) ((X.V)^2 + 4 (X.X) (V.V)) // FullSimplify
(*    True    *)

• Nice that you verified my answer. I already started to doubt. – yarchik Jul 19 at 15:20
• @yarchik this is one of my calculation errors... I have corrected the formula in my post – mattiav27 Jul 19 at 15:35
• @mattiav27 it's still wrong, the second term should be 8/25 and you need to specify that you're assuming that $\vec{X}$ has unit length. – Roman Jul 19 at 15:57

You can implement Einstein's summation convention using, for example, temporary variables as summation indices.

ClearAll[delta, CenterDot, dummyIndexQ, tensorSimplify];

SetAttributes[delta, Orderless];
SetAttributes[CenterDot, Orderless];

dummyIndexQ[x_Symbol] := MemberQ[Attributes[x], Temporary];

tensorSimplificationRules = {
delta[a_?dummyIndexQ, a_] :> 3,
delta[a_?dummyIndexQ, b_]^2 :> delta[b, b],
delta[a_, b_?dummyIndexQ] delta[b_, c_] :> delta[a, c],
delta[i_?dummyIndexQ, j_]t_[i_] :> t[j],
(t_[i_?dummyIndexQ])^2 :> (t\[CenterDot]t),
t1_[x_?dummyIndexQ] t2_[x_] :> t1\[CenterDot]t2
};

tensorSimplify[expr_] := FixedPoint[(Expand[#]//.tensorSimplificationRules)&, expr];


Let's define $$XXV_{ijk}$$:

xxv[i_, j_, k_] := (1/3 * (v[k] x[i] x[j] + v[j] x[i] x[k] + v[i] x[j] x[k]) - 1/5 * (delta[i,j] (x\[CenterDot]x) v[k] + delta[i,k] (x\[CenterDot]v) x[j] +     delta[j,k] (x\[CenterDot]v) x[i]))


Result for your example $$XXV_{abc} XXV_{abc}$$:

expr = Module[{a,b,c}, xxv[a,b,c] xxv[a,b,c]];
tensorSimplify[expr]


2/25 (v$$\cdot$$x)^2 (x$$\cdot$$x) + 8/25 (v$$\cdot$$v) (x$$\cdot$$x)^2

Result for more complicated input $$XXV_{abc} XXV_{bcd} XXV_{def} XXV_{efa}$$:

expr2 = Module[{a,b,c,d,e,f}, xxv[a,b,c] xxv[b,c,d] xxv[d,e,f] xxv[e,f,a]];
tensorSimplify[expr2]


$$\frac{524 (x\cdot x)^2 (v\cdot x)^4}{50625}+\frac{1454 v\cdot v (x\cdot x)^3 (v\cdot x)^2}{50625}+\frac{1586 (v\cdot v)^2 (x\cdot x)^4}{50625}$$

Note that all summation indices must be listed inside the first argument of Module`.