# How to implement Einstein summation convention with differential operators

Statement of this problem:

In the textbook, the following differential equilibrium equations can be expressed by tensors: Using Einstein's summation convention, the formula in the figure above can be abbreviated as follows: In addition, the strain coordination equations in the figure below can be abbreviated as: I would like to know how to implement the above summation convention with the help of MMA's tensor operator.

My problem is slightly different from this one because I have involved derivation operations and I need to use notation like $$div(σ)+F=0$$ to memorize differential equilibrium equations to reduce the burden of memorizing deformable compatibility equations.

Objectives to be addressed of this question：

I want to find a universal tensor operation function to express the equations expressed by various tensors in elasticity as shown in the figure below (thank you very much for xzczd's answer, which has made a good demonstration meeting my requirements). If I could, I would like to find a way to express the deformation compatibility equation expressed by strain or stress in a similar way as $$div(σ)+F=0$$ represents the equilibrium differential equation, so as to reduce the burden of memory.

That's the main purpose of this question. I hope I can solve this problem with your help.

What needs further explanation in the comments:

The tensor operation I mentioned mainly refers to the tensor with derivative in the textbook. For example, kl after the comma in the lower corner of $$e_{ij,kl}$$ represents the second derivative of $$e_{ij}$$. This is different from the usual tensor description.

Part of my question can also be expressed as "can I have a function that convert $$σji,j+Fi=0$$ to ".

The textbook I used didn't specify the specific meaning of the first two ee, but I saw the relevant answers, I think it should mean LeviCivitaTensor.

• Which "tensor operator" do you mean. As is, this question is very vague. You migh want to have a look at Div; it works also on matrices. Actually, a basis-free notation for the above would be $\operatorname{div}(\sigma) + F = 0$, wo why struggling with indices? Mar 22, 2020 at 10:53
• @henrik I guess OP's question can be rephrased as "Can I have a function that convert $\sigma{ji,j}+F_i=0$ to …". (Still, OP's clarification is necessary of course. ) "Why struggling with indices" Einstein summation convention appear in text books about fluid dynamics, solid mechanics, etc. quite a bit, being able to handle this with Mathematica is useful, I think. Mar 22, 2020 at 11:28
• @henrik The second set is indeed confusing. If I guess it right, there exist two kinds of $e$ here. The first two $e$ i.e. $e$ in $e_{mik}e_{nlj}$ actually denote LeviCivitaTensor, and the last two $e$ denote strain. (Once again, we still need OP's clarification, of course. ) Mar 22, 2020 at 12:11
• @Edmund Sadly it doesn't, it focus on Einstein summation convention used in list manipulation. Actually we already have a few questions related to Einstein summation convention, but none of them seems to be a duplicate of OP's question, AFAIK. Mar 22, 2020 at 12:30
• ……同学你编辑问题的时候能不能走点心。你这越编辑越乱了。 Mar 23, 2020 at 3:53

Aha, simpler than I thought. Assuming all I guessed in the comments is correct:

BeginPackage["einstein"]

allowtensor;
$tensordimension = 3; Begin["private"] expand[func_, {}] := # & expand[func_, var_] := Function[s, func[s, ##] &[Sequence @@ ({#,$tensordimension} & /@ var)], HoldAll]

tensor[index_List] :=
Function[{expr},
With[{count = Count[expr // Unevaluated, #, Infinity, Heads -> True] & /@ index},
expand[Table, Pick[index, OddQ[#] && # > 0 & /@ count]]@
expand[Sum, Pick[index, EvenQ[#] && # > 0 & /@ count]]@expr], HoldAll]

flat[expr_List] := Flatten@expr
flat[expr_] := expr

SetAttributes[allowtensor, HoldFirst]
allowtensor[a_ + b_, index_List] := allowtensor[a, index] + allowtensor[b, index]
allowtensor[c_ a_, index_List] /; FreeQ[Unevaluated@c, Alternatives @@ index] :=
c allowtensor[a, index]
allowtensor[a_ == b_, index_List] :=
flat@allowtensor[a, index] == flat@allowtensor[b, index] // Thread
allowtensor[expr_, index_List] := tensor[index][expr]

End[]
EndPackage[]


The following is not necessary, but will make the output pretty:

rule[var_] := var[i__] :> Subscript[var, Sequence @@ x /@ {i}]

drule = Derivative[id__][f_][args__] :>
HoldForm@D[f, ##] &[
Sequence @@ (DeleteCases[
Transpose[{{args}, {id}}], {_, 0}] /. {x_, 1} :> x)]];


Then let's check. Some preparation:

inde = {x, y, z};

Clear@x; x[i_] := inde[[i]];


Oh, I've used x both for function definition and independent variable, which isn't a good practice, but this is just a toy example and we know what we're doing, so let it be.

Now check the first example:

allowtensor[D[σ[i, j][x, y, z], x[j]] + F[i] == 0, {i, j}] /.
rule /@ {σ, F} /. drule The second:

ϵ = LeviCivitaTensor;

allowtensor[ϵ[[m, i, k]] ϵ[[n, l, j]] D[e[i, j][x, y, z], x[k], x[l]] == 0,
{i, j, k, l, m, n}] /.
e[i_, j_] /; i > j -> e[j, i] /. rule[e] /.
drule // DeleteDuplicates Notice the output is eliminated to 6 equations because of the symmetry, which should have been clarified in the body of question.

• Your answer is very good. I want to know if you can use differential operators and matrix operations to memorize deformable compatible equations(Similar to using $div(σ)+F=0$ to memorize differential equilibrium equation). Mar 23, 2020 at 1:25
• Why is the result of the following code True? allowtensor[(D[e[i, j][x, y, z], x[k], x[l]] + D[e[k, l][x, y, z], x[i], x[j]] - D[e[i, k][x, y, z], x[j], x[l]] - D[e[j, l][x, y, z], x[i], x[k]]) == 0, {i, j, k, l, m, n}] /. e[i_, j_] /; i > j -> e[j, i] /. rule[e] /. drule // DeleteDuplicates. According to the definition of the textbook, it should be the same as the result of $emik \;enlj\;eij, kl = 0$. Mar 23, 2020 at 2:18
• @pl 1. curl = Curl[#, {x, y, z}] &; Flatten[curl /@ Transpose[curl /@ Table[e[i, j][x, y, z], {i, 3}, {j, 3}]]] /. e[i_, j_] /; i > j -> e[j, i] /. rule[e] /. drule // DeleteDuplicates 2. There's a bug in the original implementation, now it's fixed. Thx for testing. Mar 23, 2020 at 3:51
• I'd offer extra brownie points if you can generate those to be inactive operators for FEM in NDSolve ;-) Mar 25, 2020 at 7:33
• @xzczd, this was half meant as a joke - but if you can figure this out, that would be very useful. But, please do not waste your time on it. Mar 25, 2020 at 9:52

Let me try to partially answer. Partially for the following reason: I know how to implement index vector and tensor notations and how to work with them. I also wanted to implement the Einstein convention and failed. However, even without it one can successfully use the index notations.

Let us first introduce the Kronecker, \[Delta] and Levi-Civita, ee tensors:

Subscript[δ, i_, j_] := KroneckerDelta[i, j];
Subscript[ee, i_, j_, k_] := Signature[{i, j, k}];


Let us try them. This looks as on your screen. I mean that on the screen it looks as we traditionally used to denote vectors and tensors in the index notations, but in the StackExchange is is clumsy. Therefore, in the following I include the screenshots.

Subscript[ee, 1, 2, 3]
Subscript[ee, 1, 1, 3]

(* 1

0 *)


This is the contraction of the Levi-Civita with the Kronecker tensor Sum[Subscript[ee, i, j, k]*Subscript[δ, i, k], {i, 1, 3}, {k,
1, 3}] /. j -> 3
(* 0 *)


This is the example of a vector product: Subscript[s, i_] :=
Sum[Subscript[ee, i, j, k] Subscript[a, j] Subscript[b, k], {j, 1,
3}, {k, 1, 3}];
Subscript[s, 1]

-Subscript[a, 3] Subscript[b, 2] + Subscript[a, 2] Subscript[b, 3]


Here is an example of an electrodynamics calculation of the magnetic field as a part of a Fresnel problem within this technique This is a simple example from the elasticity theory (since you seem to be interested precisely in this area): Subscript[ϵ, 1, 1] =
1/Ε*(Subscript[σ, 1,
1] - ν*(Subscript[σ, 2, 2] + Subscript[σ, 3,
3]));
Subscript[ϵ, 2, 2] =
1/Ε*(Subscript[σ, 2,
2] - ν*(Subscript[σ, 1, 1] + Subscript[σ, 3,
3]));
Subscript[ϵ, 3, 3] =
1/Ε*(Subscript[σ, 3,
3] - ν*(Subscript[σ, 1, 1] + Subscript[σ, 2,
2]));
expr = (Sum[Subscript[ϵ, i, i], {i, 1, 3}]) /.
Subscript[σ, 3,
3] -> ν*(Subscript[σ, 1, 1] + Subscript[σ, 2,
2]) // Factor

(* -(((1 + ν) (-1 + 2 ν) (Subscript[σ, 1, 1] +
Subscript[σ, 2, 2]))/Ε)  *)
`

I have more examples from elasticity theory including operating with derivatives and Green functions. However, I feel that this answer is already too long.

Have fun!

The compatibility equations (from George Herrmann's Elasticity Notes at Stanford in 1978). I think he took this course from Ray Mindlin (look him up----excellent)

This shows how to express them in "dyadic" form and Cartesian tensor form. Perhaps this, along with Alexei's nice answer will help you. I might adopt his nice notation for Kronecker and Alternating symbols.

Note as pointed out in the comments, both the alternating (Levi-Civita) and strain occur in these equations and they both have a symbol that resembles e. So make your handwriting better than mine was then.