A $n$-th-rank tensor $T_n$ based on $\mathbb{R}^3$ is referred here to as harmonic if it is symmetric in all indices and the linear map of the identity on vectors $I_2$ vanished, denoted as $T_n[I_2] = 0_{n-2}$, where $O_{n-2}$ denotes the $(n-2)$-rank zero tensor. In index notation this mean, for example, for $n=4$ $$ T_{ijkl} = T_{jikl} = T_{ijlk} = \dots \quad \land \quad \sum_{k,l=1}^3 T_{ijkl} \delta_{kl} = \sum_{k=1}^3 T_{ijkk} = 0 \quad \forall i,j \ . $$ I use the term harmonic for these tensors due to their connection to homogeneous polynomials (which are connected to symmetric $n$-th-rank tensors) with vanishing Laplacian. How would you generate a harmonic tensor of rank $n$ computationally?
I know, I can generate a symmetrized array with symbols, solve the linear map condition $T_n[I_2] = 0_{n-2}$, reinsert the solution into $T_n$ and then replace the symbols with some, e.g., random numbers. See upcoming code. But I would like to avoid the symbolic solution since I only want to generate fast some instance of harmonic tensors for high $n$. Symmetrize
and RandomReal
allow to generate way faster symmetric tensors, but I can not get my head around on how to enforce $T_n[I_2] = 0_{n-2}$ computationally. Any ideas?
(*Generate a tensor of rank n based on symbols s and with index \
symmetries is*)
genT[n_, s_, is_: {}] :=
Normal@SymmetrizedArray[
pos_ :> Subscript[s, StringJoin[ToString /@ pos]],
ConstantArray[3, n], is];
(*Identity on vectors*)
I2 = IdentityMatrix@3;
(*Linear map*)
lm[A_, B_] :=
ArrayReshape[A,
Join[Dimensions[A][[1 ;; -(1 + TensorRank@B)]], {Times @@
Dimensions@B}]].Flatten[B];
(*Generate harmonic tensor of rank n with symbols s*)
genH[n_, s_] := # /. First@Solve[lm[#, I2] == 0*lm[#, I2]] &@
genT[n, s, Symmetric@Range@n];
(*Generate a harmonic tensor of rank n*)
genHR[n_] := Module[
{s},
# /. Table[
Variables[#][[i]] -> RandomReal[], {i, Length@Variables@#}] &@
genH[n, s]
];
(*Generate instance and check properties*)
T = genHR[6]; // AbsoluteTiming
TensorSymmetry@T
lm[T, I2]
vars = Variables@T
(*Generate a symmetric tensor*)
genSR[n_] :=
Normal@Symmetrize[RandomReal[{-1, 1}, ConstantArray[3, n]]];
TensorSymmetry@genSR[6] // AbsoluteTiming