Python符号计算之数组

数组

数组是编程语言中用于复杂计算的基本数据结构，而从数学和物理上的应用来看，不同维度的数组，对应着不同类别的代数空间，即一维数组对应向量；二维对应矩阵，高维对应张量。在sympy中，与高维数组相关的大部分函数，都封装在tensor中，但支持从从sympy中直接导入

import sympy
from sympy import print_latex
from sympy import Array
a1 = Array([[1, 2], [3, 4], [5, 6]])
a1.shape    # (3,2)
a1.rank()   # 2阶张量
print_latex(a1)

$$a_1=\left[\begin{array}{cc}1& 2\\ 3& 4\\ 5& 6\end{array}\right]$$

Array类型的元素也可以是符号，示例如下

from sympy.abc import *
a2 = Array([[[x, y], [z, x*z]], [[1, x*y], [1/x, x/y]]])
a2.shape    # (2,2,2)
a2.rank()   # 3阶张量
print_latex(a2)

$$a_2=\left[\begin{array}{cc}\left[\begin{array}{cc}x& y\\ z& xz\end{array}\right]& \left[\begin{array}{cc}1& xy\\ \frac{1}{x}& \frac{x}{y}\end{array}\right]\end{array}\right]$$

张量乘法

相同维度的数组，可以直接进行元素之间的加减乘除，其中加减法无需多言，但乘法却相对复杂。除了元素之间可以进行乘法计算，将数组看成向量和矩阵，也可以进行矩阵乘法的计算，示例如下

$$\left[\begin{array}{c}c\\ d\end{array}\right]\left[\begin{array}{cc}a& b\end{array}\right]=\left[\begin{array}{cc}ca& cb\\ da& db\end{array}\right]\left[\begin{array}{cc}a& b\end{array}\right]\left[\begin{array}{c}c\\ d\end{array}\right]=ac+bd$$

若将向量看作张量，那么前者是张量乘法，后者则是在乘法之后进行了缩并操作。在sympy中，tensorproduct即为张量乘法，tensorcontraction则为缩并

from sympy import tensorproduct, tensorcontraction
A = Array([a,b])
B = Array([c,d])
p = tensorproduct(A, B)
print(sympy.latex(p))
t = tensorcontraction(p, (0,1))
print_latex(t)    # ac + bd

其中，tensorcontraction的第二个参数$(0,1)$表示被缩并的轴。

求导

sympy中实现了derive_by_array函数，这个函数可以作用在普通的标量函数表达式上，示例如下

from sympy import derive_by_array
from sympy import sin, exp
ex1 = derive_by_array(sin(x*y), x)
print_latex(ex1)

ex2 = derive_by_array(sin(x*y), [x, y, z])
print_latex(ex2)

ex1和ex2的输出分别是

$$\begin{gathered} \frac{\partial \sin (xy)}{\partial x}=y\cos \left(xy\right) \\ \nabla \sin (xy)=\left[\begin{array}{ccc}y\cos \left(xy\right)& x\cos \left(xy\right)& 0\end{array}\right] \end{gathered}$$

下面定义一阶张量$A^m=[e^x,\sin (yz),t]$，以及$x^n=[x,y,z]$，则$B^{nm}=\frac{\partial A^m}{\partial x^n}$，可以计算如下

xn = [x, y, z]
Am = [exp(x), sin(y*z), 0]
Bnm = derive_by_array(Am, xn)
print_latex(Bnm)

结果如下

$$\left[\begin{array}{ccc}{e}^x& 0& 0\\ 0& z\cos \left(yz\right)& 0\\ 0& y\cos \left(yz\right)& 0\end{array}\right]$$

如果对运算结果进行缩并，那么就可以得到$A^m$场的散度

div = tensorcontraction(Bnm, (0,1))
print_latex(div)

$$z\cos \left(yz\right)+e^x$$

即

$$\nabla (e^x,\sin (yz),t)=\frac{\partial e^x}{\partial x}+\frac{\partial \sin (yz)}{\partial y}+\frac{0}{\partial z}=z\cos \left(yz\right)+e^x$$

矢量、矩阵和张量

如果只是从数据结构出发，那么数组和张量是等价的，而矢量、矩阵则为其特例。但在具体的物理应用中，矢量、向量、张量则可能存在数据结构之外的含义。所以sympy中额外给出了这三种数据结构。

首先，矢量从编程角度来说是一维数组，但在物理中，则可用于描述指代有方向的物理量。以三维矢量为例，往往用$\hat{i},\hat{j},\hat{k}$或者${\hat{e}}_x,{\hat{e}}_y,{\hat{e}}_z$来代表$x,y,z$这三个轴向的单位矢量。正因如此，sympy中提供了专门的三维矢量坐标coordSys3D

from sympy.vector import CoordSys3D
N = CoordSys3D('N')
v = (x*y + x*z)*N.i + (y**2 + y*z)*N.j
print_latex(v)

$$\left(xy+xz\right){\hat{\mathbf{i}}}_{\mathbf{N}}+\left(y^2+yz\right){\hat{\mathbf{j}}}_{\mathbf{N}}$$

其次，矩阵除了作为一种数据组织形式外，也有丰富的物理内涵，可以表示刚体旋转等，并且在图像处理领域有着广泛应用。sympy也为其单独提供了Matrix类型

from sympy.matrices import Matrix
M = Matrix([[1,0,0], [0,0,0]])
print_latex(M)
print_latex(Matrix([M, (0, 0, -1)]))
print_latex(Matrix([[1, 2, 3]]))
print_latex(Matrix([1, 2, 3]))
m = Matrix(3,4, lambda i,j : i**2+j**2)
print_latex(m)

其中$m$是通过lambda表达式创建的矩阵，其每个元素都遵照给定函数的规则而创建，效果如下

$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right]\left[\begin{array}{ccc}1& 2& 3\end{array}\right]\left[\begin{array}{c}1\\ 2\\ 3\end{array}\right]\left[\begin{array}{cccc}0& 1& 4& 9\\ 1& 2& 5& 10\\ 4& 5& 8& 13\end{array}\right]$$

Array类型提供了tomatrix方法，可以将二维数组输出为Matrix类型，但高阶数组将会报错

a1.tomatrix()       # 返回Matrix数据
a2.tomatrix()       # 报错

最后，张量在广义相对论中也有着频繁的使用，和高维数组相比，除了在形式上多了类似$\mu \nu$之类的上下指标外，还有协变、逆变的区别，并且基于此提供了协变导数、逆变导数等运算，sympy中为此提供了TensorIndexType类型，其初始化参数有

• name : 张量类型的名字
• dummy_name : 张量指标的名字
• dim : 维度，可以是符号、整数或者None
• eps_dim : epsilon张量的维度
• metric_symmetry : 对称性
• metric_name : 矩阵名字

实例如下

from sympy.tensor.tensor import TensorIndexType
Lorentz = TensorIndexType('Lorentz', dummy_name='L')
Lorentz.metric