使用 sympy 展开索引符号方程

Question

下面我有一个使用索引符号编写的方程式。这个方程可以用图中的六个方程来表示。

第一个方程使用索引符号（爱因斯坦符号：https://en.wikipedia.org/wiki/Einstein_notation）展开。 在 U_k,k 中，逗号是导数的约定。由于我们有重复的索引 (k,k)，我们应用求和约定并得到 (du_1/dx_1 + du_2/dx_2 + du_3/dx_3)。在图中，术语 u_1、u_2 和 u_3 分别写为 u、v 和 w，它们由 x_1、x_2 和 x_3（x、y 和 z）来区分。

我是 python 和 Sympy 的新手，我看到有一个张量模块，但我看不到 Sympy 中是否已经实现了一些东西，我可以在其中编写第一个方程（带索引）并从中获得其他六个关系。

Answer 1

目前，没有直接的方法可以满足您的要求。

我将建议一个应该起作用的技巧，通过使用 IndexedBase（即，由其他符号索引的符号），然后进行替换。

声明你的符号：

U = IndexedBase("U")
l = symbols("lambda")
var("mu u v w x y z i j k")

声明你的表达式（没有爱因斯坦求和，所以明确输入一个求和）：

sij = l*Sum(U[k, k], (k, 0, 2)) * KroneckerDelta(i, j) + mu*(U[i, j] + U[j, i])

定义一个替换函数：将U[0, 0]根据索引匹配到Derivative(u, x)等：

def replace_U(uij):
    i1, i2 = uij.indices
    return Derivative([u, v, w][i1], [x, y, z][i2])

现在，遍历所有 i, j 索引，首先用整数值替换 U[i, j]，然后替换像 U[0, 0], ...

for ii in range(3):
    for ji in range(3):
        if ii < ji:
            continue
        pprint(sij.doit().xreplace({i: ii, j: ji}).replace(lambda v: v.base == U, replace_U))

记住：Sum.doit() 扩展了总和。条件 ii >= ji 用于避免重复打印相同的表达式（您的方程在 i, j 中是对称的）。

注意，上面代码中的U[i, j]只是一个符号，i和j除了作为U 的索引。替换函数为它们分配导数。

我得到的输出是：

  ⎛d       d       d    ⎞       d    
λ⋅⎜──(u) + ──(v) + ──(w)⎟ + 2⋅μ⋅──(u)
  ⎝dx      dy      dz   ⎠       dx   
  ⎛d       d    ⎞
μ⋅⎜──(u) + ──(v)⎟
  ⎝dy      dx   ⎠
  ⎛d       d       d    ⎞       d    
λ⋅⎜──(u) + ──(v) + ──(w)⎟ + 2⋅μ⋅──(v)
  ⎝dx      dy      dz   ⎠       dy   
  ⎛d       d    ⎞
μ⋅⎜──(u) + ──(w)⎟
  ⎝dz      dx   ⎠
  ⎛d       d    ⎞
μ⋅⎜──(v) + ──(w)⎟
  ⎝dz      dy   ⎠
  ⎛d       d       d    ⎞       d    
λ⋅⎜──(u) + ──(v) + ──(w)⎟ + 2⋅μ⋅──(w)
  ⎝dx      dy      dz   ⎠       dz   

我是 python 和 Sympy 的新手，我看到有一个张量模块 但我看不到 Sympy 中是否已经实现了某些东西 我可以在哪里写第一个方程（带索引）并从中获得 其他六个关系。

有 sympy.tensor.tensor，它支持抽象索引符号（一种受限的爱因斯坦求和）。不幸的是，它不支持衍生产品。

如果您想处理组件，最近添加了一个：sympy.tensor.array。它提供多维数组、张量积和收缩，以及数组的导数。

Answer 2

正如Francesco Bonazzi 所说，还有另一种使用sympy.tensor.array 来实现方程的方法。这种方法的优点是易于实现导数。

示例代码（针对jupyter）如下：

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)) * Array(eye(3)) + mu * (du + transpose(du))
for i in range(3):
    for j in range(i, 3):
        display(Eq(Symbol("sigma_{}{}".format(i+1, j+1)), sig[i, j]))

输出：

Answer 3

这是我对这个问题的回答。该代码显示/打印结果，并且由于其按索引表示法计算而不会出错。

欢迎评论进一步 (i) 适用于科学出版物的优化和打印结果。目前看起来很尴尬。

from sympy import Symbol, Derivative
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))*Array(eye(3))  + mu * (du + transpose(du))
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]))
        init_printing()