Index notation and sympy: How do we use sympy package effectively for equations are written using index notation in python?

Here is my experience to learn how to use python for a mathematical formulations using index notation and obstacles. A record can be found here at Expand index notation equation using sympy

from sympy import simplify, tensorcontraction, tensorproduct
from sympy import *
lam, mu, x, y, z = symbols("lambda mu x y z")
u, v, w = symbols("u v w", cls=Function)
du = derive_by_array([u(x, y, z), v(x, y, z), w(x, y, z)], [x, y, z])
sig = lam * tensorcontraction(du, (0, 1)) * eye(3) + mu * (du + transpose(du))
for i in range(3):
    for j in range(i, 3):
        Eq(Symbol("sigma_{}{}".format(i+1, j+1)), sig[i, j])

The error is TypeError: unsupported operand type(s) for +: ‘MutableDenseMatrix’ and ‘ImmutableDenseNDimArray’

Thank You for Your Attention!

Answer

I like to see actual sympy variables, and the traceback – or at least a clear statement of where the error occurred. I’ll show what I expect.

Running your code piece by piece in a isympy session.

In [1]: lam, mu, x, y, z = symbols("lambda mu x y z")
   ...: u, v, w = symbols("u v w", cls=Function)
   ...: du = derive_by_array([u(x, y, z), v(x, y, z), w(x, y, z)], [x, y, z])

In [2]: du
Out[2]: 
⎡∂               ∂               ∂             ⎤
⎢──(u(x, y, z))  ──(v(x, y, z))  ──(w(x, y, z))⎥
⎢∂x              ∂x              ∂x            ⎥
⎢                                              ⎥
⎢∂               ∂               ∂             ⎥
⎢──(u(x, y, z))  ──(v(x, y, z))  ──(w(x, y, z))⎥
⎢∂y              ∂y              ∂y            ⎥
⎢                                              ⎥
⎢∂               ∂               ∂             ⎥
⎢──(u(x, y, z))  ──(v(x, y, z))  ──(w(x, y, z))⎥
⎣∂z              ∂z              ∂z            ⎦

In [3]: sig = lam * tensorcontraction(du, (0, 1)) * eye(3) + mu * (du + transpos
   ...: e(du))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-181d779e2c85> in <module>
----> 1 sig = lam * tensorcontraction(du, (0, 1)) * eye(3) + mu * (du + transpose(du))

TypeError: unsupported operand type(s) for +: 'MutableDenseMatrix' and 'ImmutableDenseNDimArray'

OK, that occurs in a long expression, but I see only one +. What are the terms?

In [4]: type(du)
Out[4]: sympy.tensor.array.dense_ndim_array.ImmutableDenseNDimArray

transpose(du) is also that, and du+transpose(du) runs.

But there’s another plus – with eye(3):

In [12]: type(eye(3))
Out[12]: sympy.matrices.dense.MutableDenseMatrix

So eye(3) is the wrong type to add to the du expression.

What does the docs say about these Matrix types and their compatibility?


Change the type of the eye matrix, and the addition works:

In [19]: E3 = ImmutableDenseNDimArray(eye(3))
In [20]: E3+du;
In [22]: sig = lam * tensorcontraction(du, (0, 1)) * E3 + mu * (du + transpose(du))

sig is a matrix:

In [25]: type(sig)
Out[25]: sympy.tensor.array.dense_ndim_array.ImmutableDenseNDimArray

Your loop works, but doesn’t save or display anything. We can collect values in a list (or matrix):

In [26]: alist = []
    ...: for i in range(3):
    ...:     for j in range(i, 3):
    ...:         alist.append(Eq(Symbol("sigma_{}{}".format(i+1, j+1)), sig[i, j
    ...: ]))
    ...: 
    ...: 

In [27]: alist
Out[27]: 
⎡        ⎛∂                ∂                ∂             ⎞       ∂           
⎢σ₁₁ = λ⋅⎜──(u(x, y, z)) + ──(v(x, y, z)) + ──(w(x, y, z))⎟ + 2⋅μ⋅──(u(x, y, z
⎣        ⎝∂x               ∂y               ∂z            ⎠       ∂x          

            ⎛∂                ∂             ⎞          ⎛∂                ∂    
)), σ₁₂ = μ⋅⎜──(u(x, y, z)) + ──(v(x, y, z))⎟, σ₁₃ = μ⋅⎜──(u(x, y, z)) + ──(w(
            ⎝∂y               ∂x            ⎠          ⎝∂z               ∂x   

 ...