Basic tensor calculus with index substitution using sympy

Question

I would like to switch from Mathematica to SymPy to perform some basic index substitution in tensor product. I have large expression like $A_{ab}\times B_{bcd}\times C_{cd}$ for instance. This product can be simplified because it involves some projector using Kronecker symbol. In Mathematica, I defined this Kronecker symbol with substitution rules :

SetAttributes[\[Delta], Orderless];
\[Delta] /: \[Delta][k_, k_] = 3; 
\[Delta] /: \[Delta][k_, l_]^2 = 3;
\[Delta] /: \[Delta][k_, l_]*(f_)[m1___, k_, m2___]*x___ := x*f[m1, l, m2]; 
\[Delta] /: \[Delta][l_, k_]*(f_)[m1___, k_, m2___]*x___ := x*f[m1, l, m2];

That allows me to perform a simple index substitution like $v_{ai}\times\delta_{ij} = v_{aj}$. I can then simplify the expression and obtain the complete expression. It is the first step to further calculus.

Is it possible to define something like this in Python using SymPy? I found several Ricci packages to perform tensor calculus, but it seems way too heavy for what I want to do. I also saw some rules to substitute index to values, but I was not able to define what I want.

Answer 1

I'm not sure I fully understand what you are trying to do but sympy comes with some support for tensor expressions which might do what you want more directly:

https://docs.sympy.org/latest/modules/tensor/array_expressions.html

There is also the KroneckerDelta symbol which can be used in summations (although this might be a bit limited for what you want):

In [8]: k = symbols('k')

In [9]: s = Sum(KroneckerDelta(2, k), (k, 1, 3))

In [10]: s
Out[10]: 
  3       
 ___      
 ?        
  ?   ?   
  ?    2,k
 ?        
 ???      
k = 1     

In [11]: s.doit()
Out[11]: 1

I don't know Mathematica so well but from what I understand of the code you've shown a more direct translation of your Mathematica code would look something like this:

Delta = Function('Delta')
a, b = symbols('a, b', cls=Wild)

rules = [
    (Delta(a, a), 3),
    (Delta(a, b)**2, 3),
]

def replace_all(e):
    for r, v in rules:
        e = e.replace(r, v)
    return e

x, y = symbols('x, y')

expr = Delta(x, x) + Delta(x, y)**2

print(replace_all(expr))

This kind of pattern matching doesn't support sequence variables. Instead the usual way to do this in sympy is by using arbitrary Python functions like expr.replace(f, g) where you define f and g as functions in Python e.g.:

In [24]: is_match = lambda e: e.func == Delta and e.args[0] == e.args[1]

In [25]: replacement = lambda e: 3

In [26]: expr.replace(is_match, replacement)
Out[26]: 
 2          
? (x, y) + 3

In [27]: expr.replace(is_match, replacement)
Out[27]: 
 2          
? (x, y) + 3

Here the functions is_match and replacement could be arbitrarily complicated Python functions created with def rather than just lambda.