「SymPy」符号运算 Vector向量及运算
Posted 行吟客
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了「SymPy」符号运算 Vector向量及运算相关的知识，希望对你有一定的参考价值。
目录
导言
在前几篇文章中，我们学习了SymPy的基本使用、方程求解、微积分相关的知识。
传送链接：
「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开
「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限
「SymPy」符号运算(4) 微积分与有限差分
矢量和张量是理论物理、数值分析中极其重要的工具。现在我们来学习如何用SymPy进行向量运算。
我发现了SymPy库里有两类向量模块1：sympy.vector和sympy.physics.vector，后者是sympy.physics物理模块中的，更适用于求解物理经典力学中的动力学与运动学问题。这里先介绍sympy.vector模块。两个模块的参考文档如下：
sympy.vector模块：https://docs.sympy.org/latest/modules/vector/index.html
sympy.phsics.vector模块 ：https://docs.sympy.org/latest/modules/physics/vector/vectors.html
坐标系和向量
笛卡尔坐标系与向量
sympy.vector能够处理笛卡尔坐标系、球坐标或其他曲线坐标系统。本文只介绍三维笛卡尔坐标系统。
所有的向量vector都需要建立在坐标系上，因此首先需要初始化一个3D笛卡尔坐标系：
from sympy.vector import CoordSys3D
N = CoordSys3D('N') # 坐标系统的名字叫N
N
输出：
CoordSys3D
(
N
,
(
[
1
0
0
0
1
0
0
0
1
]
,
0
^
)
)
\\operatornameCoordSys3D\\left(N, \\left( \\left[\\beginmatrix1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\\endmatrix\\right], \\ \\mathbf\\hat0\\right)\\right)
CoordSys3D
N,
100010001
, 0^
N坐标系统下有三个正交基向量，我们可以通过i，j和k三个属性来使用
N.i
输出：
i
^
N
\\mathbf\\hati_N
i^N
当一个基向量BaseVector乘以标量或者其他SymPy表达式，就得到了VectorMul对象（相加的话得到VectorAdd对象，知道就行，都是Vector的子类）
3*N.i
输出：
(
3
)
i
^
N
(3)\\mathbf\\hati_N
(3)i^N
type(3*N.i)
# <class 'sympy.vector.vector.VectorMul'>
零向量为Vector.zero，与坐标系统无关，因此是在类sympy.vector中定义的
from sympy.vector import Vector
# 输出一个零向量
Vector.zero
# 0
type(Vector.zero)
# sympy.vector.vector.VectorZero
# 基向量+零向量
N.i + Vector.zero
# N.i
# 零向量与零向量的数乘等价
Vector.zero == 2*Vector.zero
# True
除了基向量外。每个坐标系统下自动定义了坐标变量（base scalars, or coordinate variables，注意是标量，表示坐标变量），包括X，Y和Z，定义在BaseScalar类中。可以通过坐标系的属性.x，.y或.z来使用，使用方法与一般的SymPy符号类似，如下式定义了一个势函数场
2
x
2
y
2x^2y
2x2y：
from sympy.vector import CoordSys3D
R = CoordSys3D('R')
# 一个标量电势场
electric_potential = 2*R.x**2*R.y
electric_potential
# 2*R.x**2*R.y
注意，当表达式中存在BaseScalar的实例时，暗示了存在一个随空间变化的场（field），而不包含BaseScalar实例的表达式虽然也可以看作一个场，但是一个常数场：
from sympy.vector import CoordSys3D
from sympy import symbols
R = CoordSys3D('R')
# x, y为一般的符号，并不能表示坐标系的坐标
x, y = symbols('x y')
# 一个矢量场，但是任何位置的值都是2xy
expr = 2*x*y*R.x
坐标原点
每一个坐标系统都有一个坐标原点（point）：
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
N.origin
输出：
Point
(
N
.
o
r
i
g
i
n
,
0
^
)
\\operatornamePoint\\left(N.origin, \\mathbf\\hat0\\right)
Point(N.origin,0^)
可以通过locate_new定义一个新原点，参数包括这个新原点的名字，以及新原点相对于父原点的位置向量。注意用这种方法产生的新原点的是原坐标原点的平移，没有旋转。
from sympy.abc import a, b, c
P = N.origin.locate_new('P', a*N.i + b*N.j + c*N.k)
P
输出：
Point
(
P
,
(
a
)
i
^
N
+
(
b
)
j
^
N
+
(
c
)
k
^
N
,
Point
(
N
.
o
r
i
g
i
n
,
0
^
)
)
\\operatornamePoint\\left(P, (a)\\mathbf\\hati_N + (b)\\mathbf\\hatj_N + (c)\\mathbf\\hatk_N, \\operatornamePoint\\left(N.origin, \\mathbf\\hat0\\right)\\right)
Point(P,(a)i^N+(b)j^N+(c)k^N,Point(N.origin,0^))
再定义一个原点Q：
Q = P.locate_new('Q', -b*N.j)
Q
Point ( Q , ( − b ) j ^ N , Point ( P , ( a ) i ^ N + ( b ) j ^ N + ( c ) k ^ N , Point ( N . o r i g i n , 0 ^ ) ) ) \\operatornamePoint\\left(Q, (- b)\\mathbf\\hatj_N, \\operatornamePoint\\left(P, (a)\\mathbf\\hati_N + (b)\\mathbf\\hatj_N + (c)\\mathbf\\hatk_N, \\operatornamePoint\\left(N.origin, \\mathbf\\hat0\\right)\\right)\\right) Point(Q,(−b)j^N,Point(P,(a)i^N+(b)j^N+(c)k^N,Point(N.origin,0^)))
查询两个原点之间的相对位置向量：
P.position_wrt(Q)
# b*N.j
Q.position_wrt(N.origin)
# a*N.i + c*N.k
或者查询一个原点相对另一个原点的坐标：
Q.express_coordinates(N)
# (a, 0, c)
向量运算
四则运算
+, -, *（数乘）, /，栗子：
v = N.i - 2*N.j
v/3
# 1/3*N.i + (-2/3)*N.j
v + N.k
# N.i + (-2)*N.j + N.k
Vector.zero/2
# 0
(v/3)*4
# 4/3*N.i + (-8/3)*N.j
点乘、叉乘
向量点乘函数.cot()，向量叉乘函数.cross()
v1 = 2*N.i + 3*N.j - N.k
v2 = N.i - 4*N.j + N.k
# v1点乘v2，得到标量
v1.dot(v2)
# -11
# v1叉乘v2，得到矢量
v1.cross(v2)
# (-1)*N.i + (-3)*N.j + (-11)*N.k
v2.cross(v1)
# N.i + 3*N.j + 11*N.k
注意，符号&和^被重载为了点乘和叉乘，方便直接使用（不过官方并不推荐，理由是使用函数显得清晰易懂）
v1 & v2
# -11
v1 ^ v2
# (-1)*N.i + (-3)*N.j + (-11)*N.k
并矢/外积(Dyadic)
向量并矢之后成为二阶张量。外积的函数为.outer()，或者用重载的运算符|。
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
N.i.outer(N.j)
# (N.i|N.j) # 两个基向量的并矢
N.i|N.j
# (N.i|N.j)
# 对并矢进行运算
dyad = N.i.outer(N.k)
dyad*3
# 3*(N.i|N.k)
dyad - dyad
# 0
dyad + 2*(N.j|N.i)
# (N.i|N.k) + 2*(N.j|N.i)
点乘和叉乘在并矢之间、并矢和向量之间可以正常使用，如 i j ⋅ j j = i j \\boldsymbol i \\boldsymbol j \\cdot \\boldsymbol j\\boldsymbol j = \\boldsymbol i\\boldsymbol j ij⋅jj=ij：
d = N.i.outer(N.j)
d.dot(N.j|N.j)
# (N.i|N.j)
i j ⋅ i = 0 \\boldsymbol i \\boldsymbol j \\cdot \\boldsymbol i = \\boldsymbol 0 ij⋅i=0：
d = N.i.outer(N.j)
d.dot(N.i)
# 0
k × i j = j j \\boldsymbol k \\times\\boldsymbol i \\boldsymbol j = \\boldsymbol j\\boldsymbol j k×ij=jj：
N.k ^ d
# (N.j|N.j)
∇ \\nabla ∇算子（Hamiltonian算子）
Del （deloperator）算子，又称哈密顿算子、Nabla算子，用来表示计算标量场的梯度、向量场的散度及旋度。
在SymPy中，
∇
\\nabla
∇算子通过Del类实现。
梯度
对标量场求梯度 ∇ ϕ \\nabla \\phi ∇ϕ
from sympy.vector import CoordSys3D, Del
C = CoordSys3D('C')
delo以上是关于「SymPy」符号运算 Vector向量及运算的主要内容，如果未能解决你的问题，请参考以下文章