「SymPy」符号运算(5) Vector向量及运算

行吟客

已于 2023-01-21 19:11:21 修改

阅读量1.1k

点赞数 3

分类专栏： SymPy符号运算系列文章标签：算法 python ipython scipy

于 2023-01-21 18:53:22 首次发布

本文链接：https://blog.csdn.net/weixin_45549554/article/details/128747161

版权

SymPy符号运算系列专栏收录该内容

8 篇文章 12 订阅

订阅专栏

导言

在前几篇文章中，我们学习了SymPy的基本使用、方程求解、微积分相关的知识。

传送链接：

「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
 「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开
 「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限
 「SymPy」符号运算(4) 微积分与有限差分

矢量和张量是理论物理、数值分析中极其重要的工具。现在我们来学习如何用SymPy进行向量运算。

我发现了SymPy库里有两类向量模块¹：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

输出：
$\operatorname{CoordSys3D}\left(N, \left( \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right], \ \mathbf{\hat{0}}\right)\right)$
N坐标系统下有三个正交基向量，我们可以通过i，j和k三个属性来使用

N.i

输出：
$\mathbf{\hat{i}_{N}}$
当一个基向量BaseVector乘以标量或者其他SymPy表达式，就得到了VectorMul对象（相加的话得到VectorAdd对象，知道就行，都是Vector的子类）

3*N.i

输出：
$(3)\mathbf{\hat{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符号类似，如下式定义了一个势函数场 $2x^2y$ ：

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

输出：
$\operatorname{Point}\left(N.origin, \mathbf{\hat{0}}\right)$
可以通过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

输出：
$\operatorname{Point}\left(P, (a)\mathbf{\hat{i}_{N}} + (b)\mathbf{\hat{j}_{N}} + (c)\mathbf{\hat{k}_{N}}, \operatorname{Point}\left(N.origin, \mathbf{\hat{0}}\right)\right)$
再定义一个原点Q：

Q = P.locate_new('Q', -b*N.j)
Q

$\operatorname{Point}\left(Q, (- b)\mathbf{\hat{j}_{N}}, \operatorname{Point}\left(P, (a)\mathbf{\hat{i}_{N}} + (b)\mathbf{\hat{j}_{N}} + (c)\mathbf{\hat{k}_{N}}, \operatorname{Point}\left(N.origin, \mathbf{\hat{0}}\right)\right)\right)$

查询两个原点之间的相对位置向量：

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)

点乘和叉乘在并矢之间、并矢和向量之间可以正常使用，如 $\boldsymbol i \boldsymbol j \cdot \boldsymbol j\boldsymbol j = \boldsymbol i\boldsymbol j$ ：

d = N.i.outer(N.j)
d.dot(N.j|N.j)
# (N.i|N.j)

$\boldsymbol i \boldsymbol j \cdot \boldsymbol i = \boldsymbol 0$ ：

d = N.i.outer(N.j)
d.dot(N.i)
# 0

$\boldsymbol k \times\boldsymbol i \boldsymbol j = \boldsymbol j\boldsymbol j$ ：

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')
delop = Del()                           # 获得 Del 类的实例/对象
gradient_field = delop(C.x*C.y*C.z)     # 对标量场求梯度
gradient_field
# (Derivative(C.x*C.y*C.z, C.x))*C.i + (Derivative(C.x*C.y*C.z, C.y))*C.j
# + (Derivative(C.x*C.y*C.z, C.z))*C.k
gradient_field.doit()                   # 通过 .doit() 化简
# C.y*C.z*C.i + C.x*C.z*C.j + C.x*C.y*C.k

# 另一种方法
delop.gradient(C.x*C.y*C.z).doit()
# C.y*C.z*C.i + C.x*C.z*C.j + C.x*C.y*C.k

此外，也可以使用专门用来求梯度的函数gradient()

from sympy.vector import gradient
gradient(C.x*C.y*C.z)
# C.y*C.z*C.i + C.x*C.z*C.j + C.x*C.y*C.k

旋度

对矢量场求旋度 $\nabla \times \boldsymbol F$

from sympy.vector import CoordSys3D, Del
C = CoordSys3D('C')
delop = Del()
delop.cross(C.x*C.y*C.z*C.i).doit()
# C.x*C.y*C.j + (-C.x*C.z)*C.k
(delop ^ C.x*C.y*C.z*C.i).doit()
# C.x*C.y*C.j + (-C.x*C.z)*C.k

此外，也可以使用专门用来求旋度的函数curl()

from sympy.vector import curl
curl(C.x*C.y*C.z*C.i)
# C.x*C.y*C.j + (-C.x*C.z)*C.k

散度

对矢量场求散度 $\nabla \cdot \boldsymbol F$

from sympy.vector import CoordSys3D, Del
C = CoordSys3D('C')
delop = Del()
delop.dot(C.x*C.y*C.z*(C.i + C.j + C.k)).doit()
# C.x*C.y + C.x*C.z + C.y*C.z
(delop & C.x*C.y*C.z*(C.i + C.j + C.k)).doit()
# C.x*C.y + C.x*C.z + C.y*C.z

此外，也可以使用专门用来求散度的函数divergence()

from sympy.vector import divergence
divergence(C.x*C.y*C.z*(C.i + C.j + C.k))
# C.x*C.y + C.x*C.z + C.y*C.z

除了上述三种常见的 $\nabla$ 用途，Del()还可以用来求方向导数（directional derivative）。

By definition, the directional derivative of a field $\mathbf{F}$ along a vector $\vec{v}$ at point $x$ represents the instantaneous rate of change of $\mathbf{F}$ moving through x with the velocity $\vec{v}$ . It is represented mathematically as: $(\vec{v} \cdot \nabla) \mathbf{F}$

from sympy.vector import CoordSys3D, Del
C = CoordSys3D('C')
delop = Del()
vel = C.i + C.j + C.k
scalar_field = C.x*C.y*C.z
vector_field = C.x*C.y*C.z*C.i
(vel.dot(delop))(scalar_field)
# C.x*C.y + C.x*C.z + C.y*C.z
(vel & delop)(vector_field)
# (C.x*C.y + C.x*C.z + C.y*C.z)*C.i

此外，也可以使用专门的函数directional_derivative()

from sympy.vector import directional_derivative
directional_derivative(C.x*C.y*C.z, 3*C.i + 4*C.j + C.k)
# C.x*C.y + 4*C.x*C.z + 3*C.y*C.z

向量操作

查看向量的元素，用.components属性

from sympy.vector import CoordSys3D
C = CoordSys3D('C')
v = 3*C.i + 4*C.j + 5*C.k
v.components
# {C.i: 3, C.j: 4, C.k: 5}

获得向量的模，用.magnitude()函数

v.magnitude()

输出：
$\sqrt{2}$

获得标准化后的向量，用.normalize()函数

v.normalize()

输出：
$(\frac{3 \sqrt{2}}{10})\mathbf{\hat{i}_{C}} + (\frac{2 \sqrt{2}}{5})\mathbf{\hat{j}_{C}} + (\frac{\sqrt{2}}{2})\mathbf{\hat{k}_{C}}$

（对向量内的符号）求导 diff(expr, symbol)，或者Derivative(expr, symbol)

v = (sin(a)**2 + cos(a)**2)*N.i - (2*cos(b)**2 - 1)*N.k>>> diff(v, b)
(4*sin(b)*cos(b))*N.k
# (4*sin(b)*cos(b))*N.k

from sympy import Derivative
Derivative(v, b).doit()
# (4*sin(b)*cos(b))*N.k

（对向量内的符号）求积分Integral(expr, symbol)

from sympy import Integral
v1 = a*N.i + sin(a)*N.j - N.k
Integral(v1, a)
Integral(v1, a).doit()
# a**2/2*N.i + (-cos(a))*N.j + (-a)*N.k

查询某个矢量场是否为保守场（旋度为零、标量场梯度形成的场）

from sympy.vector import CoordSys3D, is_conservative
R = CoordSys3D('R')
field = R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k
is_conservative(field)
# True
curl(field)
# 0

也可以计算某个保守场对应的标量势场（忽略了积分常数）：

from sympy.vector import CoordSys3D, scalar_potential
R = CoordSys3D('R')
conservative_field = 4*R.x*R.y*R.z*R.i + 2*R.x**2*R.z*R.j + 2*R.x**2*R.y*R.k
scalar_potential(conservative_field, R)
# 2*R.x**2*R.y*R.z

化简.simplify()、三角函数化简trigsimp()

v = (sin(a)**2 + cos(a)**2)*N.i - (2*cos(b)**2 - 1)*N.k
trigsimp(v)
# N.i + (-cos(2*b))*N.k

v.simplify()
# N.i + (-cos(2*b))*N.k

.factor()函数可以对基向量的系数进行因式分解

from sympy.abc import a, b, c
from sympy import sin, cos, 
trigsimp, diff
v = (a*b + a*c + b**2 + b*c)*N.i + N.j
v.factor()
# ((a + b)*(b + c))*N.i + N.j

向量转化为矩阵`Matrix`

通过sympy.vector的to_matrix()可以将向量转化为矩阵

from sympy.vector import *
N = CoordSys3D('N')
v = N.i + N.j + N.k
v.to_matrix(N)

输出：
$\left[\begin{matrix}1\\1\\1\end{matrix}\right]$

实践与应用

证明 $\nabla \cdot (f\boldsymbol v) = f(\nabla \cdot \boldsymbol v) + \boldsymbol v\cdot(\nabla f)$

from sympy.vector import CoordSys3D, Del
from sympy import symbols, Function

# 建立坐标系，初始化Del算子
C = CoordSys3D('C')
delop = Del()

# 定义符号与函数，v1,v2,v3为矢量v的三个分量大小，定义为C.x, C.y和C.z的函数
v1, v2, v3, f = symbols('v1 v2 v3 f', cls=Function)
vfield = v1(C.x, C.y, C.z)*C.i + v2(C.x, C.y, C.z)*C.j + v3(C.x, C.y, C.z)*C.k
ffield = f(C.x, C.y, C.z)

# 构造等式左边的表达式
lhs = (delop.dot(ffield * vfield)).doit()
# 构造等式右边的表达式
rhs = ((vfield.dot(delop(ffield))) + (ffield * (delop.dot(vfield)))).doit()

# 判断左右两边是否相等
lhs.expand().simplify() == rhs.expand().doit().simplify()
# True

向量积分的应用

要求标量场或者矢量场在某一个区域的积分，首先需要定义一个区域。SymPy可以通过参数方程、隐函数、几何模型对象来定义区域。

首先导入一些必要的常量和函数：

from sympy import sin, cos, exp, pi, symbols
from sympy.vector import CoordSys3D, ParametricRegion, ImplicitRegion, vector_integrate
from sympy.abc import r, x, y, z, theta, phi
C = CoordSys3D('C')

用到的主要函数：vector_integrate(field, *region) ，field是标量或矢量场，*region是积分区域。

计算周长、表面积、体积

为计算圆的周长，首先需要定义这个圆的边界曲线。可以用参数方程定义：

param_circle = ParametricRegion((4*cos(theta), 4*sin(theta)), (theta, 0, 2*pi))

也可以用用隐函数来定义：

implicit_circle = ImplicitRegion((x, y), x**2 + y**2 - 4)

定义的区域为圆的边界曲线，计算曲线的周长，即单位标量场的曲线积分的绝对值

vector_integrate(1, param_circle)
# 8*pi
vector_integrate(1, implicit_circle)
# 4*pi

定义一个实心球区域，并求体积：

solidsphere = ParametricRegion((r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),
	(phi, 0, pi), (theta, 0, 2*pi), (r, 0, 4))
vector_integrate(1, solidsphere)
# 256*pi/3

计算物体质量

有一个三角形区域，顶点为 $(0, 0)$ ， $(0, 5)$ ， $(5, 0)$ ，密度分布 $\rho(x, y) = xy (\mathrm{kg/m^2})$ ，求总质量：

triangle = ParametricRegion((x, y), (x, 0, 5), (y, 0, 5 - x))
vector_integrate(C.x*C.y, triangle)
# 625/24

有一个圆柱区域，以z为轴，高h，半径a，密度分布 $\rho(x, y) = x^2 + y^2 (\mathrm{kg / m^2})$ ，求总质量：

a, h = symbols('a h', positive=True)
cylinder = ParametricRegion((r*cos(theta), r*sin(theta), z),
			(theta, 0, 2*pi), (z, 0, h), (r, 0, a))
vector_integrate(C.x**2 + C.y**2, cylinder)
# pi*a**4*h/2

计算通量流量

假设有一个区域上存在常数向量场 $E(x,y,z)=a\hat{\mathbf k}$ ，有一个以 $r$ 为半径的半球落在 $x - y$ 平面，求通量：

semisphere = ParametricRegion((r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),
		(phi, 0, pi/2), (theta, 0, 2*pi))
flux = vector_integrate(a*C.k, semisphere)
flux
# pi*a*r**2

假设有一个区域上存在向量场 $E(x,y,z)=x^2\hat{\mathbf k}$ 在 $x - y$ 平面上方，向量场 $y^2 \hat{\mathbf k}$ ，求通过边长 $L$ 、中心在原点的立方体的向量场的通量：

L = symbols('L', positive=True)
top_face = ParametricRegion((x, y, L/2), (x, -L/2, L/2), (y, -L/2, L/2))
bottom_face = ParametricRegion((x, y, -L/2), (x, -L/2, L/2), (y, -L/2, L/2))
flux = vector_integrate(C.x**2*C.k, top_face) + vector_integrate(C.y**2*C.k, bottom_face)
flux
# L**4/6

证明Stokes定理

from sympy.vector import curl
curve = ParametricRegion((cos(theta), sin(theta)), (theta, 0, pi/2))
surface = ParametricRegion((r*cos(theta), r*sin(theta)), (r, 0, 1),(theta, 0, pi/2))
F = C.y*C.i + C.z*C.k + C.x*C.k
vector_integrate(F, curve)
# -pi/4
vector_integrate(curl(F), surface)
# -pi/4
>>> circle = ParametricRegion((cos(theta), sin(theta), 1), (theta, 0,2*pi))
cone = ParametricRegion((r*cos(theta), r*sin(theta), r), (r, 0, 1),(theta, 0, 2*pi))
cone = ParametricRegion((r*cos(theta), r*sin(theta), r), (r, 0, 1),(theta, 0, 2*pi))
f = (-C.y**3/3 + sin(C.x))*C.i + (C.x**3/3 + cos(C.y))*C.j + C.x*C.y*C.z*C.k
>>> vector_integrate(f, circle)
# pi/2
vector_integrate(curl(f), cone)
# pi/2

证明散度定理

from sympy.vector import divergence
sphere = ParametricRegion((4*sin(phi)*cos(theta), 4*sin(phi)*sin(theta), 4*cos(phi)),
 		   				  (phi, 0, pi), (theta, 0, 2*pi))
solidsphere = ParametricRegion((r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),
							   (r, 0, 4),(phi, 0, pi), (theta, 0, 2*pi))
field = C.x**3*C.i + C.y**3*C.j + C.z**3*C.k
vector_integrate(field, sphere)
# 12288*pi/5
vector_integrate(divergence(field), solidsphere)
# 12288*pi/5

Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103 ↩︎