SymPy 1.13 中文文档(三十)
Wigner 符号
Wigner、Clebsch-Gordan、Racah 和 Gaunt 系数
集合函数用于精确计算 Wigner 3j、6j、9j、Clebsch-Gordan、Racah 以及 Gaunt 系数,所有的计算结果都是有理数乘以有理数的平方根[Rasch03]。
请参阅个别函数的描述以获取更多详细信息和示例。
参考文献
[Regge58] (1,2)
‘Clebsch-Gordan 系数的对称性质’, T. Regge, Nuovo Cimento, 卷 10, pp. 544 (1958)
[Regge59]
‘Racah 系数的对称性质’, T. Regge, Nuovo Cimento, 卷 11, pp. 116 (1959)
[Edmonds74] (1,2,3,4,5,6,7,8,9,10)
A. R. Edmonds. Angular momentum in quantum mechanics. Investigations in physics, 4.; Investigations in physics, no. 4. Princeton, N.J., Princeton University Press, 1957.
[Rasch03] (1,2,3,4,5,6,7,8,9)
J. Rasch 和 A. C. H. Yu, ‘为预先计算的 Wigner 3j、6j 和 Gaunt 系数提供高效存储方案’, SIAM J. Sci. Comput. 卷 25, 期 4, pp. 1416-1428 (2003)
[Liberatodebrito82]
‘FORTRAN 程序计算三个球面谐波的积分’, A. Liberato de Brito, Comput. Phys. Commun., 卷 25, pp. 81-85 (1982)
[Homeier96] (1,2)
‘一些真实球面谐波耦合系数的性质及其与 Gaunt 系数的关系’, H. H. H. Homeier 和 E. O. Steinborn J. Mol. Struct., 卷 368, pp. 31-37 (1996)
Credits and Copyright
这段代码取自 Sage,经过了所有作者的允许:
groups.google.com/forum/#!topic/sage-devel/M4NZdu-7O38
作者
-
Jens Rasch (2009-03-24): Sage 的初始版本
-
Jens Rasch (2009-05-31): updated to sage-4.0
-
Oscar Gerardo Lazo Arjona (2017-06-18): added Wigner D matrices
-
Phil Adam LeMaitre (2022-09-19): 添加了真实的 Gaunt 系数
版权所有 (C) 2008 Jens Rasch jyr2000@gmail.com
sympy.physics.wigner.clebsch_gordan(j_1, j_2, j_3, m_1, m_2, m_3)
计算 Clebsch-Gordan 系数。 (\left\langle j_1 m_1 ; j_2 m_2 | j_3 m_3 \right\rangle)。
这个函数的参考资料是[Edmonds74]。
Parameters:
j_1, j_2, j_3, m_1, m_2, m_3 :
整数或半整数。
返回:
有理数乘以有理数的平方根。
示例
>>> from sympy import S
>>> from sympy.physics.wigner import clebsch_gordan
>>> clebsch_gordan(S(3)/2, S(1)/2, 2, S(3)/2, S(1)/2, 2)
1
>>> clebsch_gordan(S(3)/2, S(1)/2, 1, S(3)/2, -S(1)/2, 1)
sqrt(3)/2
>>> clebsch_gordan(S(3)/2, S(1)/2, 1, -S(1)/2, S(1)/2, 0)
-sqrt(2)/2
注记
Clebsch-Gordan 系数将通过其与 Wigner 3j 符号的关系进行评估:
[\left\langle j_1 m_1 ; j_2 m_2 | j_3 m_3 \right\rangle =(-1)^{j_1-j_2+m_3} \sqrt{2j_3+1} \operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,-m_3)]
查看 Wigner 3j 符号的文档,它展示了比 Clebsch-Gordan 系数更高的对称关系。
Authors
- Jens Rasch(2009-03-24):初始版本
sympy.physics.wigner.dot_rot_grad_Ynm(j, p, l, m, theta, phi)
返回球面谐波的旋转梯度的点积。
解释
此函数返回以下表达式的右侧:
[\vec{R}Y{_j^{p}} \cdot \vec{R}Y{l^{m}} = (-1)^{m+p} \sum\limits{k=|l-j|}^{l+j}Y{k^{m+p}} * \alpha{l,m,j,p,k} * \frac{1}{2} (k²-j²-l²+k-j-l)]
参数
j, p, l, m …. 球面谐波的指数(表达式或整数)theta, phi …. 球面谐波的角度参数
示例
>>> from sympy import symbols
>>> from sympy.physics.wigner import dot_rot_grad_Ynm
>>> theta, phi = symbols("theta phi")
>>> dot_rot_grad_Ynm(3, 2, 2, 0, theta, phi).doit()
3*sqrt(55)*Ynm(5, 2, theta, phi)/(11*sqrt(pi))
sympy.physics.wigner.gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None)
计算高斯系数。
参数:
l_1, l_2, l_3, m_1, m_2, m_3 :
整数。
prec - 精度,默认为 None。
提供精度可以大大加快计算速度。
返回:
有理数乘以有理数的平方根
(如果 prec=None),或者如果指定了精度,则为实数。
解释
高斯系数定义为三个球面谐波的积分:
[\begin{split}\begin{aligned} \operatorname{Gaunt}(l_1,l_2,l_3,m_1,m_2,m_3) &=\int Y_{l_1,m_1}(\Omega) Y_{l_2,m_2}(\Omega) Y_{l_3,m_3}(\Omega) ,d\Omega \ &=\sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} \operatorname{Wigner3j}(l_1,l_2,l_3,0,0,0) \operatorname{Wigner3j}(l_1,l_2,l_3,m_1,m_2,m_3) \end{aligned}\end{split}]
示例
>>> from sympy.physics.wigner import gaunt
>>> gaunt(1,0,1,1,0,-1)
-1/(2*sqrt(pi))
>>> gaunt(1000,1000,1200,9,3,-12).n(64)
0.006895004219221134484332976156744208248842039317638217822322799675
在 (l) 和 (m) 的非整数值上使用是错误的:
sage: gaunt(1.2,0,1.2,0,0,0)
Traceback (most recent call last):
...
ValueError: l values must be integer
sage: gaunt(1,0,1,1.1,0,-1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer
注记
高斯系数遵循以下对称规则:
-
在列的任何排列下不变
[\begin{split}\begin{aligned} Y(l_1,l_2,l_3,m_1,m_2,m_3) &=Y(l_3,l_1,l_2,m_3,m_1,m_2) \ &=Y(l_2,l_3,l_1,m_2,m_3,m_1) \ &=Y(l_3,l_2,l_1,m_3,m_2,m_1) \ &=Y(l_1,l_3,l_2,m_1,m_3,m_2) \ &=Y(l_2,l_1,l_3,m_2,m_1,m_3) \end{aligned}\end{split}]
-
在空间反射下不变,即
[Y(l_1,l_2,l_3,m_1,m_2,m_3) =Y(l_1,l_2,l_3,-m_1,-m_2,-m_3)]
-
对于 (3j) 符号的 72 个 Regge 对称性,关于对称:
-
对于不满足三角关系的 (l_1), (l_2), (l_3) 为零
-
违反条件之一为零:(l_1 \ge |m_1|), (l_2 \ge |m_2|), (l_3 \ge |m_3|)
-
仅对 (l_i) 的偶数和 (L = l_1 + l_2 + l_3 = 2n) 的 (n \in \mathbb{N}) 非零
算法
此函数使用[Liberatodebrito82]的算法精确计算 Gaunt 系数的值。注意,该公式包含大阶乘的交替求和,因此不适合有限精度算术,只适用于计算代数系统[Rasch03]。
作者
Jens Rasch(2009-03-24):Sage 的初始版本。
sympy.physics.wigner.racah(aa, bb, cc, dd, ee, ff, prec=None)
计算 Racah 符号 (W(a,b,c,d;e,f))。
参数:
a, …, f :
整数或半整数。
prec :
精度,默认为
None。提供精度可以大大加快计算速度。
返回:
有理数乘以有理数的平方根
(如果 prec=None),或者如果指定了精度,则为实数。
示例
>>> from sympy.physics.wigner import racah
>>> racah(3,3,3,3,3,3)
-1/14
注记
Racah 符号与 Wigner 6j 符号相关:
[\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6) =(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)]
请参阅 6j 符号,了解其更丰富的对称性和其他属性。
算法
此函数使用[Edmonds74]的算法来精确计算 6j 符号的值。请注意,该公式包含大因子的交错和,因此不适合有限精度算术,只适用于计算代数系统[Rasch03]。
作者
- Jens Rasch(2009-03-24):初始版本
sympy.physics.wigner.real_gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None)
计算实高斯系数。
参数:
l_1, l_2, l_3, m_1, m_2, m_3:
整数。
prec - 精度,默认:None。
提供精度可以显著加快计算速度。
返回:
有理数乘以有理数的平方根。
解释
实高斯系数被定义为对三个实球谐函数的积分:
[\begin{split}\begin{aligned} \operatorname{RealGaunt}(l_1,l_2,l_3,m_1,m_2,m_3) &=\int Z^{m_1}{l_1}(\Omega) Z^{m_2}{l_2}(\Omega) Z^{m_3}_{l_3}(\Omega) ,d\Omega \ \end{aligned}\end{split}]
或者,可以通过将实球谐函数与标准球谐函数通过酉变换 (U) 关联来定义,即 (Z^{m}{l}(\Omega)=\sum{m'}U^{m}{m'}Y^{m'}{l}(\Omega)) [Homeier96]。则实高斯系数被定义为
[\begin{split}\begin{aligned} \operatorname{RealGaunt}(l_1,l_2,l_3,m_1,m_2,m_3) &=\int Z^{m_1}{l_1}(\Omega) Z^{m_2}{l_2}(\Omega) Z^{m_3}{l_3}(\Omega) ,d\Omega \ &=\sum{m'_1 m'_2 m'3} U^{m_1}{m'1}U^{m_2}{m'2}U^{m_3}{m'_3} \operatorname{Gaunt}(l_1,l_2,l_3,m'_1,m'_2,m'_3) \end{aligned}\end{split}]
酉矩阵 (U) 的分量为
[\begin{aligned} U^m_{m'} = \delta_{|m||m'|}(\delta_{m'0}\delta_{m0} + \frac{1}{\sqrt{2}}\big[\Theta(m) \big(\delta_{m'm}+(-1)^{m'}\delta_{m'-m}\big)+i\Theta(-m)\big((-1)^{-m} \delta_{m'-m}-\delta_{m'm}(-1)^{m'-m}\big)\big]) \end{aligned}]
其中 (\delta_{ij}) 是克罗内克 δ 符号,(\Theta) 是定义为的阶跃函数
[\begin{split}\begin{aligned} \Theta(x) = \begin{cases} 1 ,\text{for}, x > 0 \ 0 ,\text{for}, x \leq 0 \end{cases} \end{aligned}\end{split}]
例子
>>> from sympy.physics.wigner import real_gaunt
>>> real_gaunt(2,2,4,-1,-1,0)
-2/(7*sqrt(pi))
>>> real_gaunt(10,10,20,-9,-9,0).n(64)
-0.00002480019791932209313156167176797577821140084216297395518482071448
对于 (l) 和 (m) 的非整数值使用是错误的::
real_gaunt(2.8,0.5,1.3,0,0,0) 追溯到(最近的调用)... ValueError: l values must be integer real_gaunt(2,2,4,0.7,1,-3.4) 追溯到(最近的调用)... ValueError: m values must be integer
注意
实高斯系数继承自标准高斯系数,在对 ((l_i, m_i)) 对的任意置换下不变,并要求 (l_i) 的和为偶数以产生非零值。它还遵循以下对称规则:
-
若 (l_1), (l_2), (l_3) 未满足条件 (l_1 \in {l_{\text{max}}, l_{\text{max}}-2, \ldots, l_{\text{min}}}), 其中 (l_{\text{max}} = l_2+l_3), 则为零
[\begin{split}\begin{aligned} l_{\text{min}} = \begin{cases} \kappa(l_2, l_3, m_2, m_3) & \text{if}, \kappa(l_2, l_3, m_2, m_3) + l_{\text{max}}, \text{is even} \ \kappa(l_2, l_3, m_2, m_3)+1 & \text{if}, \kappa(l_2, l_3, m_2, m_3) + l_{\text{max}}, \text{is odd}\end{cases} \end{aligned}\end{split}]
并且 (\kappa(l_2, l_3, m_2, m_3) = \max{\big(|l_2-l_3|, \min{\big(|m_2+m_3|, |m_2-m_3|\big)}\big)})
-
负 (m_i) 的个数为奇数时为零
算法
该函数精确计算实数 Gaunt 系数的值,使用了 [Homeier96] 和 [Rasch03] 的算法。注意,[Rasch03] 中使用的公式包含大阶乘的交替和,因此不适合有限精度算术,仅适用于计算机代数系统 [Rasch03]。然而,该函数原则上可以使用计算 Gaunt 系数的任何算法,因此在计算 Gaunt 系数的算法适用于有限精度算术的情况下也是合适的。
sympy.physics.wigner.wigner_3j(j_1, j_2, j_3, m_1, m_2, m_3)
计算 Wigner 3j 符号 (\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3))。
参数:
j_1, j_2, j_3, m_1, m_2, m_3 :
整数或半整数。
返回:
有理数乘以有理数的平方根。
例子
>>> from sympy.physics.wigner import wigner_3j
>>> wigner_3j(2, 6, 4, 0, 0, 0)
sqrt(715)/143
>>> wigner_3j(2, 6, 4, 0, 0, 1)
0
如果参数不是整数或半整数值,则出错:
sage: wigner_3j(2.1, 6, 4, 0, 0, 0)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer
sage: wigner_3j(2, 6, 4, 1, 0, -1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer or half integer
注释
Wigner 3j 符号遵循以下对称性规则:
-
在列的任何置换下不变(除了在 (J:=j_1+j_2+j_3) 处的符号变化)
[\begin{split}\begin{aligned} \operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3) &=\operatorname{Wigner3j}(j_3,j_1,j_2,m_3,m_1,m_2) \ &=\operatorname{Wigner3j}(j_2,j_3,j_1,m_2,m_3,m_1) \ &=(-1)^J \operatorname{Wigner3j}(j_3,j_2,j_1,m_3,m_2,m_1) \ &=(-1)^J \operatorname{Wigner3j}(j_1,j_3,j_2,m_1,m_3,m_2) \ &=(-1)^J \operatorname{Wigner3j}(j_2,j_1,j_3,m_2,m_1,m_3) \end{aligned}\end{split}]
-
对于空间反射是不变的,即
[\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3) =(-1)^J \operatorname{Wigner3j}(j_1,j_2,j_3,-m_1,-m_2,-m_3)]
-
对称于基于 [Regge58] 的其他 72 种对称性
-
若 (j_1), (j_2), (j_3) 不满足三角关系,则为零
-
若 (m_1 + m_2 + m_3 \neq 0), 则为零
-
违反任何一个条件均为零
(m_1 \in {-|j_1|, \ldots, |j_1|}), (m_2 \in {-|j_2|, \ldots, |j_2|}), (m_3 \in {-|j_3|, \ldots, |j_3|})
算法
该函数使用 [Edmonds74] 的算法精确计算 3j 符号的值。注意,该公式包含大阶乘的交替和,因此不适合有限精度算术,仅适用于计算机代数系统 [Rasch03]。
作者
- Jens Rasch (2009-03-24): 初始版本
sympy.physics.wigner.wigner_6j(j_1, j_2, j_3, j_4, j_5, j_6, prec=None)
计算 Wigner 6j 符号 (\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6))。
参数:
j_1, …, j_6 :
整数或半整数。
prec :
精度,默认为
None。提供精度可以大大加快计算速度。
返回:
有理数乘以有理数的平方根
(如果prec=None),或者如果给出精度,则为实数。
例子
>>> from sympy.physics.wigner import wigner_6j
>>> wigner_6j(3,3,3,3,3,3)
-1/14
>>> wigner_6j(5,5,5,5,5,5)
1/52
参数必须是整数或半整数值,并满足三角关系,否则将出错:
sage: wigner_6j(2.5,2.5,2.5,2.5,2.5,2.5)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and fulfill the triangle relation
sage: wigner_6j(0.5,0.5,1.1,0.5,0.5,1.1)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and fulfill the triangle relation
注意事项
Wigner 6j 符号与 Racah 符号有关,但展示了更多的对称性,如下所述。
[\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6) =(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)]
Wigner 6j 符号遵循以下对称规则:
-
Wigner 6j 符号在列的任何排列下都是左不变的:
[\begin{split}\begin{aligned} \operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6) &=\operatorname{Wigner6j}(j_3,j_1,j_2,j_6,j_4,j_5) \ &=\operatorname{Wigner6j}(j_2,j_3,j_1,j_5,j_6,j_4) \ &=\operatorname{Wigner6j}(j_3,j_2,j_1,j_6,j_5,j_4) \ &=\operatorname{Wigner6j}(j_1,j_3,j_2,j_4,j_6,j_5) \ &=\operatorname{Wigner6j}(j_2,j_1,j_3,j_5,j_4,j_6) \end{aligned}\end{split}]
-
它们在每两列中交换上下参数的情况下是不变的,即
[\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6) =\operatorname{Wigner6j}(j_1,j_5,j_6,j_4,j_2,j_3) =\operatorname{Wigner6j}(j_4,j_2,j_6,j_1,j_5,j_3) =\operatorname{Wigner6j}(j_4,j_5,j_3,j_1,j_2,j_6)]
-
附加 6 种对称性[Regge59]总共产生 144 种对称性
-
仅在任意三个(j)的三元组满足三角关系时非零
算法
该函数使用[Edmonds74]的算法精确计算 6j 符号的值。请注意,公式包含大量阶乘的交错和,因此不适合有限精度算术,仅适用于计算代数系统[Rasch03]。
sympy.physics.wigner.wigner_9j(j_1, j_2, j_3, j_4, j_5, j_6, j_7, j_8, j_9, prec=None)
计算 Wigner 9j 符号 (\operatorname{Wigner9j}(j_1,j_2,j_3,j_4,j_5,j_6,j_7,j_8,j_9))。
参数:
j_1, …, j_9 :
整数或半整数。
prec : 精度,默认
None. 提供精度可以大大加快计算速度。
返回:
有理数乘以有理数的平方根
(如果prec=None),或者如果给出精度,则为实数。
例子
>>> from sympy.physics.wigner import wigner_9j
>>> wigner_9j(1,1,1, 1,1,1, 1,1,0, prec=64)
0.05555555555555555555555555555555555555555555555555555555555555555
>>> wigner_9j(1/2,1/2,0, 1/2,3/2,1, 0,1,1, prec=64)
0.1666666666666666666666666666666666666666666666666666666666666667
参数必须是整数或半整数值,并满足三角关系,否则将出错:
sage: wigner_9j(0.5,0.5,0.5, 0.5,0.5,0.5, 0.5,0.5,0.5,prec=64)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and fulfill the triangle relation
sage: wigner_9j(1,1,1, 0.5,1,1.5, 0.5,1,2.5,prec=64)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and fulfill the triangle relation
算法
该函数使用[Edmonds74]的算法精确计算 3j 符号的值。请注意,公式包含大量阶乘的交错和,因此不适合有限精度算术,仅适用于计算代数系统[Rasch03]。
sympy.physics.wigner.wigner_d(J, alpha, beta, gamma)
返回角动量 J 的 Wigner D 矩阵。
返回:
代表相应的欧拉角旋转矩阵(在基础上
(J_z)的特征向量)。
[\mathcal{D}_{\alpha \beta \gamma} = \exp\big( \frac{i\alpha}{\hbar} J_z\big) \exp\big( \frac{i\beta}{\hbar} J_y\big) \exp\big( \frac{i\gamma}{\hbar} J_z\big)]
这些分量是用通用形式计算的[Edmonds74],
方程 4.1.12。
解释
J:
表示被旋转角动量空间的总角动量的整数、半整数或 SymPy 符号。
alpha, beta, gamma - 表示欧拉旋转角的实数。
环绕所谓的垂直线、节点线和图形轴的旋转角。参见[Edmonds74]。
例子
最简单的例子:
>>> from sympy.physics.wigner import wigner_d
>>> from sympy import Integer, symbols, pprint
>>> half = 1/Integer(2)
>>> alpha, beta, gamma = symbols("alpha, beta, gamma", real=True)
>>> pprint(wigner_d(half, alpha, beta, gamma), use_unicode=True)
⎡ ⅈ⋅α ⅈ⋅γ ⅈ⋅α -ⅈ⋅γ ⎤
⎢ ─── ─── ─── ───── ⎥
⎢ 2 2 ⎛β⎞ 2 2 ⎛β⎞ ⎥
⎢ ℯ ⋅ℯ ⋅cos⎜─⎟ ℯ ⋅ℯ ⋅sin⎜─⎟ ⎥
⎢ ⎝2⎠ ⎝2⎠ ⎥
⎢ ⎥
⎢ -ⅈ⋅α ⅈ⋅γ -ⅈ⋅α -ⅈ⋅γ ⎥
⎢ ───── ─── ───── ───── ⎥
⎢ 2 2 ⎛β⎞ 2 2 ⎛β⎞⎥
⎢-ℯ ⋅ℯ ⋅sin⎜─⎟ ℯ ⋅ℯ ⋅cos⎜─⎟⎥
⎣ ⎝2⎠ ⎝2⎠⎦
sympy.physics.wigner.wigner_d_small(J, beta)
返回角动量 J 的小 Wigner d 矩阵。
返回:
表示对应欧拉角旋转的矩阵
的特征向量的整数(J_z))。
[\mathcal{d}_{\beta} = \exp\big( \frac{i\beta}{\hbar} J_y\big)]
这些分量是用通用形式计算的[Edmonds74],
方程 4.1.15。
解释
J 表示被旋转角动量空间的总角动量的整数、半整数或 SymPy 符号。
被旋转角动量空间的角动量。
beta 表示欧拉角的实数。
所谓的节点线。参见[Edmonds74]。
例子
>>> from sympy import Integer, symbols, pi, pprint
>>> from sympy.physics.wigner import wigner_d_small
>>> half = 1/Integer(2)
>>> beta = symbols("beta", real=True)
>>> pprint(wigner_d_small(half, beta), use_unicode=True)
⎡ ⎛β⎞ ⎛β⎞⎤
⎢cos⎜─⎟ sin⎜─⎟⎥
⎢ ⎝2⎠ ⎝2⎠⎥
⎢ ⎥
⎢ ⎛β⎞ ⎛β⎞⎥
⎢-sin⎜─⎟ cos⎜─⎟⎥
⎣ ⎝2⎠ ⎝2⎠⎦
>>> pprint(wigner_d_small(2*half, beta), use_unicode=True)
⎡ 2⎛β⎞ ⎛β⎞ ⎛β⎞ 2⎛β⎞ ⎤
⎢ cos ⎜─⎟ √2⋅sin⎜─⎟⋅cos⎜─⎟ sin ⎜─⎟ ⎥
⎢ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎥
⎢ ⎥
⎢ ⎛β⎞ ⎛β⎞ 2⎛β⎞ 2⎛β⎞ ⎛β⎞ ⎛β⎞⎥
⎢-√2⋅sin⎜─⎟⋅cos⎜─⎟ - sin ⎜─⎟ + cos ⎜─⎟ √2⋅sin⎜─⎟⋅cos⎜─⎟⎥
⎢ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎝2⎠⎥
⎢ ⎥
⎢ 2⎛β⎞ ⎛β⎞ ⎛β⎞ 2⎛β⎞ ⎥
⎢ sin ⎜─⎟ -√2⋅sin⎜─⎟⋅cos⎜─⎟ cos ⎜─⎟ ⎥
⎣ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎝2⎠ ⎦
从表 4 中[Edmonds74]
>>> pprint(wigner_d_small(half, beta).subs({beta:pi/2}), use_unicode=True)
⎡ √2 √2⎤
⎢ ── ──⎥
⎢ 2 2 ⎥
⎢ ⎥
⎢-√2 √2⎥
⎢──── ──⎥
⎣ 2 2 ⎦
>>> pprint(wigner_d_small(2*half, beta).subs({beta:pi/2}),
... use_unicode=True)
⎡ √2 ⎤
⎢1/2 ── 1/2⎥
⎢ 2 ⎥
⎢ ⎥
⎢-√2 √2 ⎥
⎢──── 0 ── ⎥
⎢ 2 2 ⎥
⎢ ⎥
⎢ -√2 ⎥
⎢1/2 ──── 1/2⎥
⎣ 2 ⎦
>>> pprint(wigner_d_small(3*half, beta).subs({beta:pi/2}),
... use_unicode=True)
⎡ √2 √6 √6 √2⎤
⎢ ── ── ── ──⎥
⎢ 4 4 4 4 ⎥
⎢ ⎥
⎢-√6 -√2 √2 √6⎥
⎢──── ──── ── ──⎥
⎢ 4 4 4 4 ⎥
⎢ ⎥
⎢ √6 -√2 -√2 √6⎥
⎢ ── ──── ──── ──⎥
⎢ 4 4 4 4 ⎥
⎢ ⎥
⎢-√2 √6 -√6 √2⎥
⎢──── ── ──── ──⎥
⎣ 4 4 4 4 ⎦
>>> pprint(wigner_d_small(4*half, beta).subs({beta:pi/2}),
... use_unicode=True)
⎡ √6 ⎤
⎢1/4 1/2 ── 1/2 1/4⎥
⎢ 4 ⎥
⎢ ⎥
⎢-1/2 -1/2 0 1/2 1/2⎥
⎢ ⎥
⎢ √6 √6 ⎥
⎢ ── 0 -1/2 0 ── ⎥
⎢ 4 4 ⎥
⎢ ⎥
⎢-1/2 1/2 0 -1/2 1/2⎥
⎢ ⎥
⎢ √6 ⎥
⎢1/4 -1/2 ── -1/2 1/4⎥
⎣ 4 ⎦
单位制度
此模块将单位制度集成到 SymPy 中,允许用户在计算时选择使用的系统,并提供显示和转换单位的实用工具。
单位(如米、磅、秒)和常数(如光年、玻尔兹曼常数)都被视为量。Quantity对象定义了单位和物理常数(尽管它的子类PhysicalConstant可能更适合物理常数)。
数量之间的关系由它们的维度和至少另一个相同维度的量的比例因子定义。这两种类型的关系通常在UnitSystem对象内定义,除了在每个单位制度中都有效的属性。例如,1 千米在所有单位制度中都等于 1000 米,其维度在所有维度系统中都是长度。另一方面,在 SI 单位制度中,光速等于 299792458 米每秒,而在自然单位中,光速等于 1(无单位)。在 SI 和自然单位中,光速的维度为速度,但在自然单位的维度系统中,速度是无维度的,因为长度和时间是等效的。类似地,在 SI 单位制度和 CGS 及高斯单位制度之间,在电磁量的维度和比例因子中存在差异,因为后两种单位制度不认为电流是一个基本维度。
与其他库中的实现相比,此实现的优势在于,它以不同的方式处理单位制度之间的关系,而不受 SI 单位制度对单位和物理常数关系的假设限制。
示例
单位模块中最重要的函数是convert_to,它允许将给定的量重新表示为某些目标量的幂的乘积。例如,要用米和秒表示光速:
>>> from sympy.physics.units import speed_of_light, meter, second
>>> from sympy.physics.units import convert_to
>>> convert_to(speed_of_light, [meter, second])
299792458*meter/second
如果无法用目标单位表示给定的量,将原样返回给定的量:
>>> convert_to(speed_of_light, [meter])
speed_of_light
数量之间的关系取决于单位制度。因此,convert_to接受一个可选的第三个参数,表示单位制度,默认为SI。根据所选的单位制度,转换可能返回不同的结果,例如,在cgs_gauss单位制度中,电流不是一个基本维度,而是可以表示为长度、时间和质量的组合:
>>> from sympy.physics.units.systems.si import SI
>>> from sympy.physics.units.systems.cgs import cgs_gauss
>>> from sympy.physics.units import ampere, gram, second
>>> convert_to(ampere, [meter, gram, second], SI)
ampere
>>> convert_to(ampere, [meter, gram, second], cgs_gauss)
149896229*sqrt(gram)*meter**(3/2)/(50*second**2)
相同维度的量不会自动简化,例如如果你将米除以千米,你会得到一个表示两个单位之间除法的对象。为了简化这类表达式,你可以调用.simplify()方法或导入quantity_simplify()函数,后者还可以接受一个单位系统作为可选参数。
>>> from sympy.physics.units.util import quantity_simplify
>>> from sympy.physics.units import kilometer
>>> meter/kilometer
meter/kilometer
>>> (meter/kilometer).simplify()
1/1000
>>> quantity_simplify(meter/kilometer)
1/1000
更多
有关未来发展的想法可以在Github wiki中找到。
-
单位制背后的哲学
-
更多例子
-
维度和维度系统
-
单位前缀
-
单位和单位制度
-
物理量
单位制背后的哲学
原文链接:
docs.sympy.org/latest/modules/physics/units/philosophy.html
维度
介绍
单位制背后的根源是维度系统,其结构主要决定了单位系统的结构。我们的定义可能看起来粗略,但对于我们的目的来说已经足够了。
维度被定义为可测量的并分配给特定现象的属性。在这个意义上,维度与纯数不同,因为它们具有一些额外的意义,因此不能将两个不同的维度相加。例如,时间或长度是维度,但对于我们有意义的任何其他事物,如角度、粒子数(摩尔...)或信息(比特...)也是如此。
从这个角度来看,唯一真正无量纲的量是纯数。无量纲的概念非常依赖于系统,正如在((c, \hbar, G))中所见的那样,所有单位在通常的常识中似乎都是无量纲的。这在通用单位系统的可计算性上是不可避免的(但最终我们可以告诉程序什么是无量纲的)。
通过取其乘积或其比率(在下文中定义)可以将维度组合在一起。例如,速度定义为长度除以时间,或者我们可以将长度看作速度乘以时间,取决于我们认为哪个更基本:一般来说,我们可以选择一组基础维度,从中我们可以描述所有其他维度。
群结构
在这个简短的介绍之后,旨在从直观的角度介绍维度之后,我们描述了数学结构。具有(n)个独立维度({d_i}_{i=1,\ldots,n})的维度系统由乘法群(G)描述:
-
存在一个纯数对应的单位元素(1);
-
两个元素(D_1, D_2 \in G)的乘积(D_3 = D_1 D_2)也在(G)中;
-
任何元素(D \in G)都有逆元(D^{-1} \in G)。
我们表示
[D^n = \underbrace{D \times \cdots \times D}_{\text{ times}},]
并且按定义(D⁰ = 1)。称为群生成元的({d_i}_{i=1,\ldots,n}),因为群中的任何元素(D \in G)都可以表示为生成元的幂的乘积:
[D = \prod_{i=1}^n d_i^{a_i}, \qquad a_i \in \mathbf{Z}.]
对于(a_i = 0, \forall i)给出了单位元,而对于(a_i = 1, a_j = 0, \forall j \neq i)我们恢复了生成元(d_i)。该群具有以下特性:
-
阿贝尔的,因为生成元交换,([d_i, d_j] = 0);
-
可数(无限但离散),因为元素按生成元的幂进行索引[1]。
可以通过取旧生成元的某些组合来改变维度基({d'i}{i=1,\ldots,n}):
[d'i = \prod{j=1}^n d_j^{P_{ij}}.]
线性空间表示
可以使用线性空间 (\mathbf{Z}^n) 作为群的表示,因为幂次系数 (a_i) 包含了所需的所有信息(我们不区分群的元素和其表示):
[\begin{split}(d_i)j = \delta{ij}, \qquad D = \begin{pmatrix} a_1 \ \vdots \ a_n \end{pmatrix}.\end{split}]
到 (d'i) 的基变换遵循线性空间的通常基变换规则,矩阵由新向量的系数 (P{ij}) 给出,这些系数简单地是旧基础下新向量的系数:
[d'i = P{ij} d_j.]
我们将在算法中使用这个最后的解决方案。
一个例子
为了说明所有这些形式主义,我们用一个具体例子来结束本节,即 MKSA 系统(m, kg, s),其维度为 (L: length, M: mass, T: time)。它们表示为(我们将始终按字母顺序排列向量)
[\begin{split}L = \begin{pmatrix} 1 \ 0 \ 0 \end{pmatrix}, \qquad M = \begin{pmatrix} 0 \ 1 \ 0 \end{pmatrix}, \qquad T = \begin{pmatrix} 0 \ 0 \ 1 \end{pmatrix}.\end{split}]
其他维度可以导出,例如速度 (V) 或作用量 (A)
[\begin{split}V = L T^{-1}, \qquad A = M L² T^{-2},\ V = \begin{pmatrix} 1 \ 0 \ -1 \end{pmatrix}, \qquad A = \begin{pmatrix} 2 \ 1 \ -2 \end{pmatrix}.\end{split}]
我们可以转换基础以转向自然系统 ((m, c, \hbar)),其维度为 (L: length, V: velocity, A: action) [2]。在此基础上,生成器为
[\begin{split}A = \begin{pmatrix} 1 \ 0 \ 0 \end{pmatrix}, \qquad L = \begin{pmatrix} 0 \ 1 \ 0 \end{pmatrix}, \qquad V = \begin{pmatrix} 0 \ 0 \ 1 \end{pmatrix},\end{split}]
而质量和时间分别由
[\begin{split}T = L V^{-1}, \qquad M = A V^{-2},\ T = \begin{pmatrix} 0 \ 1 \ -1 \end{pmatrix}, \qquad M = \begin{pmatrix} 1 \ 0 \ -2 \end{pmatrix}.\end{split}]
最终逆变换基矩阵 (P^{-1}) 通过将在旧基础下表达的向量粘合在一起而得到:
[\begin{split}P^{-1} = \begin{pmatrix} 2 & 1 & 1 \ 1 & 0 & 0 \ -2 & 0 & -1 \end{pmatrix}.\end{split}]
要找到基矩阵的变换,我们只需取其逆
[\begin{split}P = \begin{pmatrix} 0 & 1 & 0 \ 1 & 0 & 1 \ 0 & -2 & -1 \end{pmatrix}.\end{split}]
数量
一个量由其名称、维度和到相同维度的规范量的因子定义。规范量是单位模块的内部参考,不应影响最终用户。单位和物理常数都是数量。
单位
单位,如米、秒和千克,通常是人们选择的参考量,用于引用其他数量。
定义了几种不同维度的单位后,我们可以形成一个单位制度,这基本上是一个带有比例概念的维度系统。
常数
物理常数只是数量。它们表明我们以前并不理解两个维度实际上是相同的。例如,我们看到光速不是 1,因为我们没有意识到时间和空间是相同的(这是因为我们的感官;但在基本层面上它们是不同的)。例如,曾经有一个“热常数”,它允许在焦耳和卡路里之间进行转换,因为人们不知道热是能量。一旦他们理解了这一点,他们把这个常数固定为 1(这是一个非常简略的故事)。
现在我们固定国际单位制中基本常数的值,这表明它们是单位(我们用它们来定义其他常用单位)。
参考需求
不可能从头定义单位和单位系统:我们需要定义一些参考点,然后在它们上面建立其他内容。换句话说,我们需要一个刻度起源来定义我们单位的尺度(即一个因子为 1 的单位),并确保所有给定维度的单位都以一致的方式定义,这可以发生在我们希望在另一个系统中使用一个派生单位作为基本单位时:我们不应将其定义为具有比例 1,因为即使在系统内部是不一致的,我们也无法将其转换为第一个系统,因为从我们的角度来看,我们有两个不同单位的相同比例(这意味着它们在计算机中是相等的)。
我们将说,在系统外定义的维度和刻度是规范的,因为我们用它们进行所有计算。另一方面,相对于系统得到的维度和刻度称为物理的,因为它们最终带有一种意义。
让我们举一个具体(而且重要)的例子:质量单位的情况。我们希望将克定义为起点。我们希望将克定义为质量的规范起点,因此我们赋予它一个比例 1。然后我们可以定义一个系统(例如化学系统),将其作为基本单位。MKS 系统更喜欢使用千克;一个简单的选择是给它一个比例为 1 的标度,因为它是一个基本单位,但我们看到我们无法将其转换为化学系统,因为克和千克都被赋予了相同的因子。因此,我们需要把千克定义为 1000 克,然后才能在 MKS 系统中使用它作为基础。但是一旦我们问“千克在 MKS 中的因子是多少?”,我们得到的答案是 1,因为它是一个基本单位。
因此,我们将定义所有的计算而不涉及任何系统,在最后一步我们才将结果插入到系统中,以给出我们感兴趣的上下文。
文献
[Page52]
C. H. Page,SI 单位的类别,物理学美国杂志,20 卷,1 期(1952 年):1。
[Page78]
C. H. Page,物理学中的单位与维度,物理学美国杂志,46 卷,1 期(1978 年):78。
[deBoer79]
J. de Boer, 量和单位的群属性, Am. J. of Phys. 47, 9 (1979): 818.
[LevyLeblond77]
J.-M. Lévy-Leblond, 关于物理常数的概念性质, La Rivista Del Nuovo Cimento 7, no. 2 (1977): 187-214.
[NIST]
脚注
更多例子
原文链接:
docs.sympy.org/latest/modules/physics/units/examples.html
在接下来的几节中,我们将提供一些可以使用此模块完成的示例。
尺寸分析
我们将从牛顿第二定律开始
[m a = F]
其中 (m, a) 和 (F) 分别是质量、加速度和力。知道 (m) ((M)) 和 (a) ((L T^{-2})) 的维度后,我们将确定 (F) 的维度;显然,我们将发现它是一个力:(M L T^{-2})。
从那里我们将使用质量为 (m) 的粒子和质量为 (M) 的物体之间的引力表达式,距离为 (r)。
[F = \frac{G m M}{r²}]
以确定牛顿引力常数 (G) 的维度。结果应为 (L³ M^{-1} T^{-2})。
>>> from sympy import symbols >>> from sympy.physics.units.systems import SI >>> from sympy.physics.units import length, mass, acceleration, force >>> from sympy.physics.units import gravitational_constant as G >>> from sympy.physics.units.systems.si import dimsys_SI >>> F = mass*acceleration >>> F Dimension(acceleration*mass) >>> dimsys_SI.get_dimensional_dependencies(F) {Dimension(length): 1, Dimension(mass, M): 1, Dimension(time): -2} >>> dimsys_SI.get_dimensional_dependencies(force) {Dimension(length): 1, Dimension(mass): 1, Dimension(time): -2}尽管在国际单位制中它们相同,但尺寸不能直接比较:
>>> F == force False尺寸系统对象提供了测试尺寸等效性的方法:
>>> dimsys_SI.equivalent_dims(F, force) True>>> m1, m2, r = symbols("m1 m2 r") >>> grav_eq = G * m1 * m2 / r**2 >>> F2 = grav_eq.subs({m1: mass, m2: mass, r: length, G: G.dimension}) >>> F2 Dimension(mass*length*time**-2) >>> F2.get_dimensional_dependencies() {'length': 1, 'mass': 1, 'time': -2}
注意应先解方程,然后用尺寸进行替换。
具有数量的方程
使用开普勒第三定律
[\frac{T²}{a³} = \frac{4 \pi²}{GM}]
我们可以使用从维基百科获取的其他变量的已知值来找到金星的轨道周期。结果应为 224.701 天。
>>> from sympy import solve, symbols, pi, Eq >>> from sympy.physics.units import Quantity, length, mass >>> from sympy.physics.units import day, gravitational_constant as G >>> from sympy.physics.units import meter, kilogram >>> T = symbols("T") >>> a = Quantity("venus_a")在国际单位制中指定维度和比例:
>>> SI.set_quantity_dimension(a, length) >>> SI.set_quantity_scale_factor(a, 108208000e3*meter)添加太阳质量作为量:
>>> M = Quantity("solar_mass") >>> SI.set_quantity_dimension(M, mass) >>> SI.set_quantity_scale_factor(M, 1.9891e30*kilogram)现在是开普勒定律:
>>> eq = Eq(T**2 / a**3, 4*pi**2 / G / M) >>> eq Eq(T**2/venus_a**3, 4*pi**2/(gravitational_constant*solar_mass)) >>> q = solve(eq, T)[1] >>> q 2*pi*venus_a**(3/2)/(sqrt(gravitational_constant)*sqrt(solar_mass))
要转换为天数,使用 convert_to 函数(可能需要近似结果):
>>> from sympy.physics.units import convert_to
>>> convert_to(q, day)
71.5112118495813*pi*day
>>> convert_to(q, day).n()
224.659097795948*day
我们也可以使用来自天体物理系统的太阳质量和日子作为单位,但我们想展示如何创建一个所需的单位。
我们可以看到在这个例子中,中间维度可能不明确,比如 sqrt(G),但应检查最终结果 - 当所有维度组合在一起时 - 是否明确定义。
维度和维度系统
原文:
docs.sympy.org/latest/modules/physics/units/dimensions.html
物理维度的定义。
单位制将建立在这些维度之上。
文档中的大多数示例使用 MKS 系统,并且从计算机的角度来看:从人类的角度来看,在 MKS 中将长度添加到时间是不合法的,但在自然系统中是合法的;对于计算机在自然系统中不存在时间维度(而是速度维度代替)- 在基础中 - 因此将时间添加到长度的问题毫无意义。
class sympy.physics.units.dimensions.Dimension(name, symbol=None)
此类表示物理量的维度。
Dimension 构造函数以名称和可选符号作为参数。
例如,在经典力学中,我们知道时间与温度不同,并且维度使得这种差异明显(但它们不提供这些量的任何测量。
>>> from sympy.physics.units import Dimension
>>> length = Dimension('length')
>>> length
Dimension(length)
>>> time = Dimension('time')
>>> time
Dimension(time)
可以使用乘法、除法和指数运算(乘以数)来组合维度,以生成新的维度。仅当两个对象为相同维度时才定义加法和减法。
>>> velocity = length / time
>>> velocity
Dimension(length/time)
可以使用维度系统对象获取维度的维度依赖性,例如,可以使用 SI 单位约定中使用的维度系统:
>>> from sympy.physics.units.systems.si import dimsys_SI
>>> dimsys_SI.get_dimensional_dependencies(velocity)
{Dimension(length, L): 1, Dimension(time, T): -1}
>>> length + length
Dimension(length)
>>> l2 = length**2
>>> l2
Dimension(length**2)
>>> dimsys_SI.get_dimensional_dependencies(l2)
{Dimension(length, L): 2}
has_integer_powers(dim_sys)
检查维度对象是否仅具有整数幂。
所有维度幂应为整数,但在中间步骤中可能出现有理数幂。此方法可用于检查最终结果是否定义良好。
class sympy.physics.units.dimensions.DimensionSystem(base_dims, derived_dims=(), dimensional_dependencies={})
DimensionSystem 表示一组一致的维度。
构造函数接受三个参数:
-
基础维度;
-
派生维度:这些是以基础维度定义的(例如,速度是通过长度除以时间定义的);
-
维度依赖性:派生维度如何依赖于基础维度。
可选地,derived_dims 或 dimensional_dependencies 可能会被省略。
property can_transf_matrix
无用的方法,保持与先前版本的兼容性。
请勿使用。
返回从规范到基础维度基础的规范变换矩阵。
它是使用 inv_can_transf_matrix() 计算的矩阵的逆。
property dim
无用的方法,保持与先前版本的兼容性。
请勿使用。
给出系统的维度。
这是返回形成基础的维度数量。
dim_can_vector(dim)
无用的方法,保持与先前版本的兼容性。
请勿使用。
以规范基础维度表示的维度。
dim_vector(dim)
无用的方法,保持与先前版本的兼容性。
请勿使用。
以基础维度表示的向量。
property inv_can_transf_matrix
无用的方法,保持与先前版本的兼容性。
请勿使用。
计算从基础到规范维度基础的逆变换矩阵。
它对应于矩阵,其中列是规范基础维度的向量。
这个矩阵几乎不会被使用,因为维度总是相对于规范基定义的,因此不需要额外工作来在此基础上获取它们。尽管如此,如果此矩阵不是方阵(或不可逆),这意味着我们选择了一个不好的基。
property is_consistent
这个方法无用,仅为了与之前版本兼容而保留。
请勿使用。
检查系统是否定义良好。
is_dimensionless(dimension)
检查维度对象是否确实具有维度。
维度应该至少有一个具有非零幂的分量。
property list_can_dims
这个方法无用,仅为了与之前版本兼容而保留。
请勿使用。
列出所有规范维度名称。
print_dim_base(dim)
给出维度的字符串表达式,用基本符号表示。
单位前缀
原文:
docs.sympy.org/latest/modules/physics/units/prefixes.html
定义单位前缀类和一些常数的模块。
SI 和二进制前缀的常数字典被定义为 PREFIXES 和 BIN_PREFIXES。
class sympy.physics.units.prefixes.Prefix(name, abbrev, exponent, base=10, latex_repr=None)
这个类表示前缀,带有它们的名称、符号和因子。
前缀用于从给定单位创建导出单位。它们应始终封装到单位中。
该因子是从一个基数(默认为 10)构造到某个幂,并给出总倍数或分数。例如,千米 km 是从米(因子 1)和千(10 的 3 次方,即 1000)构造而成。基数可以更改以允许例如二进制前缀。
一个前缀乘以另一个对象总是返回另一个对象乘以这个因子的乘积,除非另一个对象:
-
是一个前缀,它们可以组合成一个新的前缀;
-
定义与前缀的乘法(这是单位类的情况)。
单位和单位制度
原文链接:
docs.sympy.org/latest/modules/physics/units/unitsystem.html
物理量的单位制;包括常数的定义。
class sympy.physics.units.unitsystem.UnitSystem(base_units, units=(), name='', descr='', dimension_system=None, derived_units: Dict[Dimension, Quantity] = {})
UnitSystem 表示一个连贯的单位集合。
单位系统基本上是一个具有比例概念的维度系统。许多方法都以相同的方式定义。
如果所有基本单位都有符号,那就更好了。
property dim
给出系统的维度。
这是返回形成基础的单位数量。
extend(base, units=(), name='', description='', dimension_system=None, derived_units: Dict[Dimension, Quantity] = {})
将当前系统扩展到一个新系统。
取当前系统的基本和标准单位,将它们与参数中给出的基本和标准单位合并。如果未提供,则名称和描述被覆盖为空字符串。
get_units_non_prefixed() → Set[Quantity]
返回该系统中没有前缀的单位。
property is_consistent
检查底层维度系统是否一致。
物理量
原文:
docs.sympy.org/latest/modules/physics/units/quantities.html
物理量。
class sympy.physics.units.quantities.Quantity(name, abbrev=None, latex_repr=None, pretty_unicode_repr=None, pretty_ascii_repr=None, mathml_presentation_repr=None, is_prefixed=False, **assumptions)
物理数量:可以是测量单位、常量或通用数量。
property abbrev
表示单位名称的符号。
如果定义了缩写词,则在前缀符号之前加上缩写词。
convert_to(other, unit_system='SI')
将数量转换为具有相同维度的另一个数量。
示例
>>> from sympy.physics.units import speed_of_light, meter, second
>>> speed_of_light
speed_of_light
>>> speed_of_light.convert_to(meter/second)
299792458*meter/second
>>> from sympy.physics.units import liter
>>> liter.convert_to(meter**3)
meter**3/1000
property free_symbols
返回数量的无返回符号。
property is_prefixed
数量是否带有前缀。例如,(kilogram) 带有前缀,但 (gram) 没有。
property scale_factor
相对于规范单位的整体数量。
set_global_relative_scale_factor(scale_factor, reference_quantity)
设置在所有单位系统中都有效的比例因子。
量之间的转换
几种简化涉及单位对象的表达式的方法。
sympy.physics.units.util.convert_to(expr, target_units, unit_system='SI')
将 expr 转换为其所有单位和数量表示为 target_units 的因子的相同表达式,只要维度兼容。
target_units 可以是单个单位/数量,也可以是单位/数量的集合。
示例
>>> from sympy.physics.units import speed_of_light, meter, gram, second, day
>>> from sympy.physics.units import mile, newton, kilogram, atomic_mass_constant
>>> from sympy.physics.units import kilometer, centimeter
>>> from sympy.physics.units import gravitational_constant, hbar
>>> from sympy.physics.units import convert_to
>>> convert_to(mile, kilometer)
25146*kilometer/15625
>>> convert_to(mile, kilometer).n()
1.609344*kilometer
>>> convert_to(speed_of_light, meter/second)
299792458*meter/second
>>> convert_to(day, second)
86400*second
>>> 3*newton
3*newton
>>> convert_to(3*newton, kilogram*meter/second**2)
3*kilogram*meter/second**2
>>> convert_to(atomic_mass_constant, gram)
1.660539060e-24*gram
转换为多个单位:
>>> convert_to(speed_of_light, [meter, second])
299792458*meter/second
>>> convert_to(3*newton, [centimeter, gram, second])
300000*centimeter*gram/second**2
转换为普朗克单位:
>>> convert_to(atomic_mass_constant, [gravitational_constant, speed_of_light, hbar]).n()
7.62963087839509e-20*hbar**0.5*speed_of_light**0.5/gravitational_constant**0.5
高能物理
伽玛矩阵
处理表示为张量对象的伽玛矩阵的模块。
示例
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
>>> from sympy.tensor.tensor import tensor_indices
>>> i = tensor_indices('i', LorentzIndex)
>>> G(i)
GammaMatrix(i)
请注意,四维空间中已经存在一个 GammaMatrixHead 实例:GammaMatrix,它只需声明为
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix
>>> from sympy.tensor.tensor import tensor_indices
>>> i = tensor_indices('i', LorentzIndex)
>>> GammaMatrix(i)
GammaMatrix(i)
访问度规张量
>>> LorentzIndex.metric
metric(LorentzIndex,LorentzIndex)
sympy.physics.hep.gamma_matrices.extract_type_tens(expression, component)
从 TensExpr 中提取所有具有 (component) 的张量。
返回两个张量表达式:
-
第一个包含所有
Tensor具有 (component)。 -
第二个包含所有其余。
sympy.physics.hep.gamma_matrices.gamma_trace(t)
一行伽玛矩阵的痕迹
示例
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, gamma_trace, LorentzIndex
>>> from sympy.tensor.tensor import tensor_indices, tensor_heads
>>> p, q = tensor_heads('p, q', [LorentzIndex])
>>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
>>> ps = p(i0)*G(-i0)
>>> qs = q(i0)*G(-i0)
>>> gamma_trace(G(i0)*G(i1))
4*metric(i0, i1)
>>> gamma_trace(ps*ps) - 4*p(i0)*p(-i0)
0
>>> gamma_trace(ps*qs + ps*ps) - 4*p(i0)*p(-i0) - 4*p(i0)*q(-i0)
0
sympy.physics.hep.gamma_matrices.kahane_simplify(expression)
此函数取消四维伽玛矩阵乘积中的收缩元素,导致一个等于给定表达式的表达式,没有收缩的伽玛矩阵。
参数:
expression 包含要简化的伽玛矩阵的张量表达式。
注意事项
如果给出旋量指标,则矩阵必须按照乘积中给定的顺序给出。
算法
该算法背后的思想是使用一些众所周知的身份,即用于包围偶数个 (\gamma) 矩阵的收缩
(\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N}} \gamma_\mu = 2 (\gamma_{a_{2N}} \gamma_{a_1} \cdots \gamma_{a_{2N-1}} + \gamma_{a_{2N-1}} \cdots \gamma_{a_1} \gamma_{a_{2N}} ))
对于奇数个 (\gamma) 矩阵
(\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N+1}} \gamma_\mu = -2 \gamma_{a_{2N+1}} \gamma_{a_{2N}} \cdots \gamma_{a_{1}})
而不是重复应用这些身份来取消所有收缩的指数,可以识别这种操作将导致的链接,因此问题简化为自由伽玛矩阵的简单重新排列。
示例
使用时,请记住原始表达式的系数必须单独处理
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
>>> from sympy.physics.hep.gamma_matrices import kahane_simplify
>>> from sympy.tensor.tensor import tensor_indices
>>> i0, i1, i2 = tensor_indices('i0:3', LorentzIndex)
>>> ta = G(i0)*G(-i0)
>>> kahane_simplify(ta)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
>>> tb = G(i0)*G(i1)*G(-i0)
>>> kahane_simplify(tb)
-2*GammaMatrix(i1)
>>> t = G(i0)*G(-i0)
>>> kahane_simplify(t)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
>>> t = G(i0)*G(-i0)
>>> kahane_simplify(t)
Matrix([
[4, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 4, 0],
[0, 0, 0, 4]])
如果没有收缩,将返回相同的表达式
>>> tc = G(i0)*G(i1)
>>> kahane_simplify(tc)
GammaMatrix(i0)*GammaMatrix(i1)
参考文献
[1] 降低伽玛矩阵收缩乘积的算法,Joseph Kahane,数学物理学杂志,第 9 卷,第 10 期,1968 年 10 月。
sympy.physics.hep.gamma_matrices.simplify_gpgp(ex, sort=True)
简化乘积 G(i)*p(-i)*G(j)*p(-j) -> p(i)*p(-i)
示例
>>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex, simplify_gpgp
>>> from sympy.tensor.tensor import tensor_indices, tensor_heads
>>> p, q = tensor_heads('p, q', [LorentzIndex])
>>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
>>> ps = p(i0)*G(-i0)
>>> qs = q(i0)*G(-i0)
>>> simplify_gpgp(ps*qs*qs)
GammaMatrix(-L_0)*p(L_0)*q(L_1)*q(-L_1)
物理向量模块
物理/向量参考资料
[WikiDyadics]
“Dyadics.” Wikipedia, the Free Encyclopedia. Web. 05 Aug. 2011. <en.wikipedia.org/wiki/Dyadics>.
[WikiDyadicProducts]
“Dyadic Product.” Wikipedia, the Free Encyclopedia. Web. 05 Aug. 2011. <en.wikipedia.org/wiki/Dyadic_product>.
[Likins1973]
Likins, Peter W. Elements of Engineering Mechanics. McGraw-Hill, Inc. 1973. Print.
向量指南
-
向量与参考系
-
向量:运动学
-
物理/向量模块的潜在问题/高级主题/未来功能
-
标量和向量场功能
-
物理向量 API
-
基本类
-
运动学(文档字符串)
-
打印(文档字符串)
-
基本函数(文档字符串)
-
基本场函数文档字符串
-
向量与参考框架
原文链接:
docs.sympy.org/latest/modules/physics/vector/vectors.html
在sympy.physics.vector中,向量和参考框架是动态系统的“基础”。本文档将从数学角度描述它们,并说明如何在此模块的代码中使用它们。
Vector
向量是具有大小(或长度)和方向的几何对象。在纸上,三维空间中的向量通常表示为:
image/svg+xml 页面上的向量 页面外的向量 页面内的向量
向量代数
向量代数是首要讨论的主题。
如果两个向量具有相同的大小和方向,则它们被称为相等。
Vector Operations
可以对向量进行多种代数操作:向量之间的加法,标量乘法和向量乘法。
向量加法是基于平行四边形法则。
image/svg+xml a b b a a+b
向量加法也是可交换的:
[\begin{split}\mathbf{a} + \mathbf{b} &= \mathbf{b} + \mathbf{a} \ (\mathbf{a} + \mathbf{b}) + \mathbf{c} &= \mathbf{a} + (\mathbf{b} + \mathbf{c})\end{split}]
标量乘法是向量与标量的乘积;结果是一个方向相同但大小按标量缩放的向量。注意,乘以-1 相当于围绕垂直于向量平面的任意轴旋转向量 180 度。
image/svg+xml a 2a -a
单位向量简单地说是其大小等于 1 的向量。对于任意向量(\mathbf{v}),我们可以定义一个单位向量:
[\mathbf{\hat{n}_v} = \frac{\mathbf{v}}{\Vert \mathbf{v} \Vert}]
注意,每个向量都可以写成标量和单位向量的乘积。
在sympy.physics.vector中实现了三种向量积:点积、叉积和外积。
点积运算将两个向量映射到一个标量。它的定义如下:
[\begin{split}\mathbf{a} \cdot \mathbf{b} = \Vert \mathbf{a} \Vert \Vert \mathbf{b} \Vert \cos(\theta)\\end{split}]
其中(\theta)为(\mathbf{a})和(\mathbf{b})之间的角度。
两个单位向量的点积代表了共同方向的大小;对于其他向量,它是共同方向的大小和两个向量大小的乘积。两个垂直向量的点积为零。下图显示了一些示例:
image/svg+xml a b a a a c 45° a ⋅ b = 0 a ⋅ a = 1 a ⋅ c = 1/sqrt(2)
点乘是交换的:
[\mathbf{a} \cdot \mathbf{b} = \mathbf{b} \cdot \mathbf{a}]
两个向量的叉乘向量乘法操作返回一个向量:
[\mathbf{a} \times \mathbf{b} = \mathbf{c}]
向量 (\mathbf{c}) 具有以下特性:它的方向与 (\mathbf{a}) 和 (\mathbf{b}) 垂直,其大小定义为 (\Vert \mathbf{c} \Vert = \Vert \mathbf{a} \Vert \Vert \mathbf{b} \Vert \sin(\theta))(其中 (\theta) 是 (\mathbf{a}) 和 (\mathbf{b}) 之间的角度),并且其方向由右手法则定义在 (\Vert \mathbf{a} \Vert \Vert \mathbf{b} \Vert) 之间。下图显示如下:
image/svg+xml a x a = 0 a a b d c c / sqrt(2) 45°
叉乘具有以下属性:
它不是交换的:
[\begin{split}\mathbf{a} \times \mathbf{b} &\neq \mathbf{b} \times \mathbf{a} \ \mathbf{a} \times \mathbf{b} &= - \mathbf{b} \times \mathbf{a}\end{split}]
并非可结合的:
[(\mathbf{a} \times \mathbf{b} ) \times \mathbf{c} \neq \mathbf{a} \times (\mathbf{b} \times \mathbf{c})]
两个平行向量的叉积为零。
两个向量的外积不在此讨论,而是在惯性部分(使用的地方)。其他有用的向量性质和关系包括:
[\begin{split}\alpha (\mathbf{a} + \mathbf{b}) &= \alpha \mathbf{a} + \alpha \mathbf{b}\ \mathbf{a} \cdot (\mathbf{b} + \mathbf{c}) &= \mathbf{a} \cdot \mathbf{b} + \mathbf{a} \cdot \mathbf{c}\ \mathbf{a} \times (\mathbf{b} + \mathbf{c}) &= \mathbf{a} \times \mathbf{b} + \mathbf{a} \times \mathbf{b}\ (\mathbf{a} \times \mathbf{b}) \cdot \mathbf{c} & \textrm{ 给出标量三重积。}\ \mathbf{a} \times (\mathbf{b} \cdot \mathbf{c}) & \textrm{ 不起作用,因为不能对向量和标量进行叉乘。}\ (\mathbf{a} \times \mathbf{b}) \cdot \mathbf{c} &= \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})\ (\mathbf{a} \times \mathbf{b}) \cdot \mathbf{c} &= (\mathbf{b} \times \mathbf{c}) \cdot \mathbf{a} = (\mathbf{c} \times \mathbf{a}) \cdot \mathbf{b}\ (\mathbf{a} \times \mathbf{b}) \times \mathbf{c} &= \mathbf{b}(\mathbf{a} \cdot \mathbf{c}) - \mathbf{a}(\mathbf{b} \cdot \mathbf{c})\ \mathbf{a} \times (\mathbf{b} \times \mathbf{c}) &= \mathbf{b}(\mathbf{a} \cdot \mathbf{c}) - \mathbf{c}(\mathbf{a} \cdot \mathbf{b})\\end{split}]
替代表示
如果我们有三个非共面的单位向量 (\mathbf{\hat{n}_x},\mathbf{\hat{n}_y},\mathbf{\hat{n}_z}),我们可以将任意向量 (\mathbf{a}) 表示为 (\mathbf{a} = a_x \mathbf{\hat{n}_x} + a_y \mathbf{\hat{n}_y} + a_z \mathbf{\hat{n}_z})。在这种情况下,(\mathbf{\hat{n}_x},\mathbf{\hat{n}_y},\mathbf{\hat{n}_z}) 被称为基底。(a_x, a_y, a_z) 被称为测量数。通常单位向量是相互垂直的,这时我们可以称它们为正交基,通常是右手的。
现在我们可以通过以下方式测试两个向量的相等性。对于向量:
[\begin{split}\mathbf{a} &= a_x \mathbf{\hat{n}_x} + a_y \mathbf{\hat{n}_y} + a_z \mathbf{\hat{n}_z}\ \mathbf{b} &= b_x \mathbf{\hat{n}_x} + b_y \mathbf{\hat{n}_y} + b_z \mathbf{\hat{n}_z}\\end{split}]
我们可以在以下情况下宣称相等:(a_x = b_x, a_y = b_y, a_z = b_z)。
同样的两个向量的向量加法表示为:
[\mathbf{a} + \mathbf{b} = (a_x + b_x)\mathbf{\hat{n}_x} + (a_y + b_y) \mathbf{\hat{n}_y} + (a_z + b_z) \mathbf{\hat{n}_z}]
现在定义乘法运算如下:
[\begin{split}\alpha \mathbf{b} &= \alpha b_x \mathbf{\hat{n}_x} + \alpha b_y \mathbf{\hat{n}_y} + \alpha b_z \mathbf{\hat{n}_z}\ \mathbf{a} \cdot \mathbf{b} &= a_x b_x + a_y b_y + a_z b_z\ \mathbf{a} \times \mathbf{b} &= \textrm{det }\begin{bmatrix} \mathbf{\hat{n}_x} & \mathbf{\hat{n}_y} & \mathbf{\hat{n}_z} \ a_x & a_y & a_z \ b_x & b_y & b_z \end{bmatrix}\ (\mathbf{a} \times \mathbf{b}) \cdot \mathbf{c} &= \textrm{det }\begin{bmatrix} a_x & a_y & a_z \ b_x & b_y & b_z \ c_x & c_y & c_z \end{bmatrix}\\end{split}]
要在给定的基础上写出一个向量,我们可以这样做:
[\begin{split}\mathbf{a} = (\mathbf{a}\cdot\mathbf{\hat{n}_x})\mathbf{\hat{n}_x} + (\mathbf{a}\cdot\mathbf{\hat{n}_y})\mathbf{\hat{n}_y} + (\mathbf{a}\cdot\mathbf{\hat{n}_z})\mathbf{\hat{n}_z}\\end{split}]
Examples
Some numeric examples of these operations follow:
[\begin{split}\mathbf{a} &= \mathbf{\hat{n}_x} + 5 \mathbf{\hat{n}_y}\ \mathbf{b} &= \mathbf{\hat{n}_y} + \alpha \mathbf{\hat{n}_z}\ \mathbf{a} + \mathbf{b} &= \mathbf{\hat{n}_x} + 6 \mathbf{\hat{n}_y} + \alpha \mathbf{\hat{n}_z}\ \mathbf{a} \cdot \mathbf{b} &= 5\ \mathbf{a} \cdot \mathbf{\hat{n}_y} &= 5\ \mathbf{a} \cdot \mathbf{\hat{n}_z} &= 0\ \mathbf{a} \times \mathbf{b} &= 5 \alpha \mathbf{\hat{n}_x} - \alpha \mathbf{\hat{n}_y} + \mathbf{\hat{n}_z}\ \mathbf{b} \times \mathbf{a} &= -5 \alpha \mathbf{\hat{n}_x} + \alpha \mathbf{\hat{n}_y} - \mathbf{\hat{n}_z}\\end{split}]
Vector Calculus
处理带有移动物体的向量微积分,我们必须引入参考框架的概念。一个经典的例子是火车沿着轨道运行,你和朋友都在里面。如果你和朋友都坐着,那么你们之间的相对速度为零。从火车外部的观察者来看,你们两个都会有速度。
我们现在将更严格地应用这一定义。一个参考框架是我们选择从中观察向量量的虚拟“平台”。如果我们有一个参考框架(\mathbf{N}),向量(\mathbf{a})在框架(\mathbf{N})中被认为是固定的,如果从(\mathbf{N})观察时它的任何特性都不会改变。我们通常会为每个参考框架分配一个固定的正交基向量集;(\mathbf{N})将具有(\mathbf{\hat{n}_x}, \mathbf{\hat{n}_y},\mathbf{\hat{n}_z})作为其基向量。
Derivatives of Vectors
一个在参考框架中不固定的向量,因此在从该框架观察时具有不同的特性。微积分是研究变化的学科,为了处理不同参考框架中固定和不固定向量的特殊性,我们需要在定义上更加明确。
image/svg+xml A B f e d c c xd xe f 修复在: xx A B
在上图中,我们有向量 (\mathbf{c,d,e,f})。如果我们对 (\mathbf{e}) 求关于 (\theta) 的导数:
[\frac{d \mathbf{e}}{d \theta}]
目前尚不清楚导数是什么。如果你是从框架 (\mathbf{A}) 观察的话,导数显然不为零。如果你是从框架 (\mathbf{B}) 观察的话,导数为零。因此,我们将引入框架作为导数符号的一部分:
[\begin{split}\frac{^{\mathbf{A}} d \mathbf{e}}{d \theta} &\neq 0 \textrm{,在参考框架 } \mathbf{A} \textrm{ 中,} \mathbf{e} \textrm{ 对 } \theta \textrm{ 的导数不为零。}\ \frac{^{\mathbf{B}} d \mathbf{e}}{d \theta} &= 0 \textrm{,在参考框架 } \mathbf{B} \textrm{ 中,} \mathbf{e} \textrm{ 对 } \theta \textrm{ 的导数为零。}\ \frac{^{\mathbf{A}} d \mathbf{c}}{d \theta} &= 0 \textrm{,在参考框架 } \mathbf{A} \textrm{ 中,} \mathbf{c} \textrm{ 对 } \theta \textrm{ 的导数为零。}\ \frac{^{\mathbf{B}} d \mathbf{c}}{d \theta} &\neq 0 \textrm{,在参考框架 } \mathbf{B} \textrm{ 中,} \mathbf{c} \textrm{ 对 } \theta \textrm{ 的导数不为零。}\\end{split}]
下面是特定参考框架中向量导数的一些额外性质:
[\begin{split}\frac{^{\mathbf{A}} d}{dt}(\mathbf{a} + \mathbf{b}) &= \frac{^{\mathbf{A}} d\mathbf{a}}{dt} + \frac{^{\mathbf{A}} d\mathbf{b}}{dt}\ \frac{^{\mathbf{A}} d}{dt}\gamma \mathbf{a} &= \frac{ d \gamma}{dt}\mathbf{a} + \gamma\frac{^{\mathbf{A}} d\mathbf{a}}{dt}\ \frac{^{\mathbf{A}} d}{dt}(\mathbf{a} \times \mathbf{b}) &= \frac{^{\mathbf{A}} d\mathbf{a}}{dt} \times \mathbf{b} + \mathbf{a} \times \frac{^{\mathbf{A}} d\mathbf{b}}{dt}\\end{split}]
关联基向量集合
现在我们需要定义两个不同参考框架之间的关系;或者如何将一个框架的基向量与另一个框架的基向量相关联。我们可以使用方向余弦矩阵(DCM)来实现这一点。方向余弦矩阵将一个框架的基向量与另一个框架的基向量相关联,如下所示:
[\begin{split}\begin{bmatrix} \mathbf{\hat{a}_x} \ \mathbf{\hat{a}_y} \ \mathbf{\hat{a}_z} \ \end{bmatrix} = \begin{bmatrix} ^{\mathbf{A}} \mathbf{C}^{\mathbf{B}} \end{bmatrix} \begin{bmatrix} \mathbf{\hat{b}_x} \ \mathbf{\hat{b}_y} \ \mathbf{\hat{b}_z} \ \end{bmatrix}\end{split}]
当两个框架(比如 (\mathbf{A}) 和 (\mathbf{B}))最初对齐时,其中一个框架的所有基向量围绕与一个基向量对齐的轴旋转,我们称这些框架通过简单旋转相关联。下图展示了这一点:
image/svg+xml A B θ θ azbz ax bx ay by
上述旋转是绕 Z 轴简单旋转,角度为 (\theta)。请注意,在旋转后,基向量 (\mathbf{\hat{a}_z}) 和 (\mathbf{\hat{b}_z}) 仍然对齐。
这个旋转可以由以下方向余弦矩阵来表征:
[\begin{split}^{\mathbf{A}}\mathbf{C}^{\mathbf{B}} = \begin{bmatrix} \cos(\theta) & - \sin(\theta) & 0\ \sin(\theta) & \cos(\theta) & 0\ 0 & 0 & 1\ \end{bmatrix}\end{split}]
简单的绕 X 和 Y 轴的旋转由以下方式定义:
[ \begin{align}\begin{aligned}\begin{split}\textrm{绕 x 轴旋转的 DCM: } \begin{bmatrix} 1 & 0 & 0\ 0 & \cos(\theta) & -\sin(\theta)\ 0 & \sin(\theta) & \cos(\theta) \end{bmatrix}\end{split}\\begin{split}\textrm{绕 y 轴旋转的 DCM: } \begin{bmatrix} \cos(\theta) & 0 & \sin(\theta)\ 0 & 1 & 0\ -\sin(\theta) & 0 & \cos(\theta)\ \end{bmatrix}\end{split}\end{aligned}\end{align} ]
正方向的旋转通过右手规则定义。
方向余弦矩阵还涉及到基向量组之间点积的定义。如果我们有两个带有关联基向量的参考坐标系,它们的方向余弦矩阵可以定义为:
[\begin{split}\begin{bmatrix} C_{xx} & C_{xy} & C_{xz}\ C_{yx} & C_{yy} & C_{yz}\ C_{zx} & C_{zy} & C_{zz}\ \end{bmatrix} = \begin{bmatrix} \mathbf{\hat{a}_x}\cdot\mathbf{\hat{b}_x} & \mathbf{\hat{a}_x}\cdot\mathbf{\hat{b}_y} & \mathbf{\hat{a}_x}\cdot\mathbf{\hat{b}_z}\ \mathbf{\hat{a}_y}\cdot\mathbf{\hat{b}_x} & \mathbf{\hat{a}_y}\cdot\mathbf{\hat{b}_y} & \mathbf{\hat{a}_y}\cdot\mathbf{\hat{b}_z}\ \mathbf{\hat{a}_z}\cdot\mathbf{\hat{b}_x} & \mathbf{\hat{a}_z}\cdot\mathbf{\hat{b}_y} & \mathbf{\hat{a}_z}\cdot\mathbf{\hat{b}_z}\ \end{bmatrix}\end{split}]
此外,方向余弦矩阵是正交的,即:
[\begin{split}^{\mathbf{A}}\mathbf{C}^{\mathbf{B}} = (^{\mathbf{B}}\mathbf{C}^{\mathbf{A}})^{-1}\ = (^{\mathbf{B}}\mathbf{C}^{\mathbf{A}})^T\\end{split}]
如果我们有参考坐标系 (\mathbf{A}) 和 (\mathbf{B}),在这个例子中经历了简单的绕 Z 轴旋转,旋转角度为 (\theta),我们将会有两组基向量。然后我们可以定义两个向量:(\mathbf{a} = \mathbf{\hat{a}_x} + \mathbf{\hat{a}_y} + \mathbf{\hat{a}_z}) 和 (\mathbf{b} = \mathbf{\hat{b}_x} + \mathbf{\hat{b}_y} + \mathbf{\hat{b}_z})。如果我们希望在 (\mathbf{A}) 坐标系中表达 (\mathbf{b}),我们执行以下操作:
[\begin{split}\mathbf{b} &= \mathbf{\hat{b}_x} + \mathbf{\hat{b}_y} + \mathbf{\hat{b}_z}\ \mathbf{b} &= \begin{bmatrix}\mathbf{\hat{a}_x}\cdot (\mathbf{\hat{b}_x} + \mathbf{\hat{b}_y} + \mathbf{\hat{b}_z})\end{bmatrix} \mathbf{\hat{a}_x} + \begin{bmatrix}\mathbf{\hat{a}_y}\cdot (\mathbf{\hat{b}_x} + \mathbf{\hat{b}_y} + \mathbf{\hat{b}_z})\end{bmatrix} \mathbf{\hat{a}_y} + \begin{bmatrix}\mathbf{\hat{a}_z}\cdot (\mathbf{\hat{b}_x} + \mathbf{\hat{b}_y} + \mathbf{\hat{b}_z})\end{bmatrix} \mathbf{\hat{a}_z}\ \mathbf{b} &= (\cos(\theta) - \sin(\theta))\mathbf{\hat{a}_x} + (\sin(\theta) + \cos(\theta))\mathbf{\hat{a}_y} + \mathbf{\hat{a}_z}\end{split}]
如果我们希望在B中表示(\mathbf{a}),我们可以这样做:
[\begin{split}\mathbf{a} &= \mathbf{\hat{a}_x} + \mathbf{\hat{a}_y} + \mathbf{\hat{a}_z}\ \mathbf{a} &= \begin{bmatrix}\mathbf{\hat{b}_x}\cdot (\mathbf{\hat{a}_x} + \mathbf{\hat{a}_y} + \mathbf{\hat{a}_z})\end{bmatrix} \mathbf{\hat{b}_x} + \begin{bmatrix}\mathbf{\hat{b}_y}\cdot (\mathbf{\hat{a}_x} + \mathbf{\hat{a}_y} + \mathbf{\hat{a}_z})\end{bmatrix} \mathbf{\hat{b}_y} + \begin{bmatrix}\mathbf{\hat{b}_z}\cdot (\mathbf{\hat{a}_x} + \mathbf{\hat{a}_y} + \mathbf{\hat{a}_z})\end{bmatrix} \mathbf{\hat{b}_z}\ \mathbf{a} &= (\cos(\theta) + \sin(\theta))\mathbf{\hat{b}_x} + (-\sin(\theta)+\cos(\theta))\mathbf{\hat{b}_y} + \mathbf{\hat{b}_z}\end{split}]
多个参考系下的导数
如果我们有参考系A和B,我们将会有两组基向量。然后我们可以定义两个向量:(\mathbf{a} = a_x\mathbf{\hat{a}_x} + a_y\mathbf{\hat{a}_y} + a_z\mathbf{\hat{a}_z}) 和 (\mathbf{b} = b_x\mathbf{\hat{b}_x} + b_y\mathbf{\hat{b}_y} + b_z\mathbf{\hat{b}_z})。如果我们想要在参考系A中对(\mathbf{b})进行导数运算,我们必须首先在A中表示它,然后对测量数字进行导数运算:
[\frac{^{\mathbf{A}} d\mathbf{b}}{dx} = \frac{d (\mathbf{b}\cdot \mathbf{\hat{a}_x} )}{dx} \mathbf{\hat{a}_x} + \frac{d (\mathbf{b}\cdot \mathbf{\hat{a}_y} )}{dx} \mathbf{\hat{a}_y} + \frac{d (\mathbf{b}\cdot \mathbf{\hat{a}_z} )}{dx} \mathbf{\hat{a}_z} +]
向量微积分的示例
一个矢量微积分的例子:
image/svg+xml A B l x ax ay bx by c θ
在本例中,我们有两个物体,每个物体都有一个附加的参考框架。我们将(\theta)和(x)视为时间的函数。我们希望知道向量(\mathbf{c})在(\mathbf{A})和(\mathbf{B})框架中的时间导数。
首先,我们需要定义(\mathbf{c});(\mathbf{c}=x\mathbf{\hat{b}_x}+l\mathbf{\hat{b}_y})。这在(\mathbf{B})框架中提供了一个定义。现在我们可以执行以下操作:
[\begin{split}\frac{^{\mathbf{B}} d \mathbf{c}}{dt} &= \frac{dx}{dt} \mathbf{\hat{b}_x} + \frac{dl}{dt} \mathbf{\hat{b}_y}\ &= \dot{x} \mathbf{\hat{b}_x}\end{split}]
要在(\mathbf{A})框架中进行导数运算,我们必须首先将这两个框架联系起来:
[\begin{split}^{\mathbf{A}} \mathbf{C} ^{\mathbf{B}} = \begin{bmatrix} \cos(\theta) & 0 & \sin(\theta)\ 0 & 1 & 0\ -\sin(\theta) & 0 & \cos(\theta)\ \end{bmatrix}\end{split}]
现在我们可以执行以下操作:
[\begin{split}\frac{^{\mathbf{A}} d \mathbf{c}}{dt} &= \frac{d (\mathbf{c} \cdot \mathbf{\hat{a}_x})}{dt} \mathbf{\hat{a}_x} + \frac{d (\mathbf{c} \cdot \mathbf{\hat{a}_y})}{dt} \mathbf{\hat{a}_y} + \frac{d (\mathbf{c} \cdot \mathbf{\hat{a}_z})}{dt} \mathbf{\hat{a}_z}\ &= \frac{d (\cos(\theta) x)}{dt} \mathbf{\hat{a}_x} + \frac{d (l)}{dt} \mathbf{\hat{a}_y} + \frac{d (-\sin(\theta) x)}{dt} \mathbf{\hat{a}_z}\ &= (-\dot{\theta}\sin(\theta)x + \cos(\theta)\dot{x}) \mathbf{\hat{a}_x} + (\dot{\theta}\cos(\theta)x + \sin(\theta)\dot{x}) \mathbf{\hat{a}_z}\end{split}]
注意,这是(\mathbf{A})框架中(\mathbf{c})的时间导数,并以(\mathbf{A})框架表达。然而,我们也可以在(\mathbf{B})框架中表达它,并且表达式仍然有效:
[\begin{split}\frac{^{\mathbf{A}} d \mathbf{c}}{dt} &= (-\dot{\theta}\sin(\theta)x + \cos(\theta)\dot{x}) \mathbf{\hat{a}_x} + (\dot{\theta}\cos(\theta)x + \sin(\theta)\dot{x}) \mathbf{\hat{a}_z}\ &= \dot{x}\mathbf{\hat{b}_x} - \theta x \mathbf{\hat{b}_z}\\end{split}]
注意两种表达形式在复杂度上的差异。它们是等价的,但其中一种要简单得多。这是一个非常重要的概念,因为在更复杂的形式中定义向量可能会大大减慢动力学方程的制定并增加其长度,有时甚至无法在屏幕上显示。
使用向量和参考框架
在所有相关数学关系被定义之后,我们才引入代码。这是由于向量的形成方式。在开始任何问题之前,必须定义一个参考框架(记得首先导入sympy.physics.vector):
>>> from sympy.physics.vector import *
>>> N = ReferenceFrame('N')
现在我们已经创建了一个参考框架(\mathbf{N})。要访问任何基向量,首先需要创建一个参考框架。现在我们已经创建了表示(\mathbf{N})的对象,我们可以访问它的基向量:
>>> N.x
N.x
>>> N.y
N.y
>>> N.z
N.z
物理向量中的向量代数
现在我们可以对这些向量进行基本的代数运算。
>>> N.x == N.x
True
>>> N.x == N.y
False
>>> N.x + N.y
N.x + N.y
>>> 2 * N.x + N.y
2*N.x + N.y
记住,不要将标量数量添加到向量中(N.x + 5);这将引发错误。在此时,我们将在我们的向量中使用 SymPy 的符号。在处理符号时,请记住参考 SymPy 的陷阱和问题。
>>> from sympy import Symbol, symbols
>>> x = Symbol('x')
>>> x * N.x
x*N.x
>>> x*(N.x + N.y)
x*N.x + x*N.y
在sympy.physics.vector中,已经实现了多个接口来进行矢量乘法,包括操作符级别、方法级别和函数级别。矢量点乘可以如下工作:
>>> N.x.dot(N.x)
1
>>> N.x.dot(N.y)
0
>>> dot(N.x, N.x)
1
>>> dot(N.x, N.y)
0
“官方”接口是函数接口;这是所有示例中将使用的内容。这是为了避免混淆,使属性和方法彼此紧邻,并且在操作符操作优先级的情况下。在sympy.physics.vector中用于矢量乘法的操作符没有正确的操作顺序;这可能导致错误。在使用操作符表示矢量乘法时,需要注意括号。
叉乘是另一种将讨论的矢量乘法。它提供了与点乘类似的接口,并伴随着相同的警告。
>>> N.x.cross(N.x)
0
>>> N.x.cross(N.z)
- N.y
>>> cross(N.x, N.y)
N.z
>>> cross(N.x, (N.y + N.z))
- N.y + N.z
向量还可以进行两种额外的操作:将向量归一化为长度 1,以及获取其大小。操作如下进行:
>>> (N.x + N.y).normalize()
sqrt(2)/2*N.x + sqrt(2)/2*N.y
>>> (N.x + N.y).magnitude()
sqrt(2)
向量通常以矩阵形式表达,特别是在数值计算中。由于矩阵形式不包含向量所定义的参考框架的任何信息,因此必须提供参考框架以从向量中提取测量数。有一个方便的函数可以实现这一点:
>>> (x * N.x + 2 * x * N.y + 3 * x * N.z).to_matrix(N)
Matrix([
[ x],
[2*x],
[3*x]])
物理向量中的矢量微积分
我们已经介绍了我们的第一个参考框架。如果我们愿意,我们可以在该框架中立即进行微分:
>>> (x * N.x + N.y).diff(x, N)
N.x
SymPy 有一个diff函数,但它目前无法与sympy.physics.vector中的向量一起使用,请使用Vector的diff方法。原因在于,在对Vector进行微分时,除了对什么进行微分外,还必须指定参考坐标系;SymPy 的diff函数不适合这种情况。
更有趣的情况出现在多个参考框架中。如果我们引入第二个参考框架(\mathbf{A}),我们现在有了两个框架。请注意,此时我们可以添加(\mathbf{N})和(\mathbf{A})的分量,但不能执行矢量乘法,因为尚未定义两个框架之间的关系。
>>> A = ReferenceFrame('A')
>>> A.x + N.x
N.x + A.x
如果我们想进行矢量乘法,首先必须定义一个方向。ReferenceFrame的orient方法提供了这个功能。
>>> A.orient(N, 'Axis', [x, N.y])
通过简单的绕 Y 轴旋转的方式,将(\mathbf{A})框架相对于(\mathbf{N})框架进行定位,旋转量为 x。这两个框架之间的方向余弦矩阵可以随时通过dcm方法查看:A.dcm(N)给出了方向余弦矩阵 (^{\mathbf{A}} \mathbf{C} ^{\mathbf{N}})。
其他更复杂的旋转类型包括体旋转、空间旋转、四元数和任意轴旋转。体和空间旋转等效于连续进行三个简单旋转,每个旋转围绕新框架中的基向量。以下是一个例子:
>>> N = ReferenceFrame('N')
>>> Bp = ReferenceFrame('Bp')
>>> Bpp = ReferenceFrame('Bpp')
>>> B = ReferenceFrame('B')
>>> q1,q2,q3 = symbols('q1 q2 q3')
>>> Bpp.orient(N,'Axis', [q1, N.x])
>>> Bp.orient(Bpp,'Axis', [q2, Bpp.y])
>>> B.orient(Bp,'Axis', [q3, Bp.z])
>>> N.dcm(B)
Matrix([
[ cos(q2)*cos(q3), -sin(q3)*cos(q2), sin(q2)],
[sin(q1)*sin(q2)*cos(q3) + sin(q3)*cos(q1), -sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), -sin(q1)*cos(q2)],
[sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3), sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2)]])
>>> B.orient(N,'Body',[q1,q2,q3],'XYZ')
>>> N.dcm(B)
Matrix([
[ cos(q2)*cos(q3), -sin(q3)*cos(q2), sin(q2)],
[sin(q1)*sin(q2)*cos(q3) + sin(q3)*cos(q1), -sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), -sin(q1)*cos(q2)],
[sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3), sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2)]])
空间方向与体方向类似,但是应用于从框架到体。体和空间旋转可以涉及两个或三个轴:‘XYZ’适用,‘YZX’、‘ZXZ’、‘YXY’等也适用。关键在于每个简单旋转是围绕不同轴进行的,与前一个不同;‘ZZX’不能完全定向一个三维空间的基向量集。
有时,创建一个新的参考框架并在一步中相对于现有框架进行定向会更方便。orientnew方法允许此功能,并且本质上包装了orient方法。在orient中可以做的所有事情,在orientnew中也可以做到。
>>> C = N.orientnew('C', 'Axis', [q1, N.x])
四元数(或欧拉参数)使用 4 个值来描述框架的方向。这些内容以及任意轴旋转可以在orient和orientnew方法的帮助下找到,或者在参考文献[Kane1983]中了解。
最后,在开始多框架计算操作之前,我们将介绍另一个sympy.physics.vector工具:dynamicsymbols。dynamicsymbols是在 SymPy 中创建时间未定义函数的快捷方式。这样一个‘dynamicsymbol’的导数如下所示。
>>> from sympy import diff
>>> q1, q2, q3 = dynamicsymbols('q1 q2 q3')
>>> diff(q1, Symbol('t'))
Derivative(q1(t), t)
上面的‘dynamicsymbol’打印结果不太清晰;我们还将在这里介绍其他一些工具。在非交互式会话中,我们可以使用vprint代替 print。
>>> q1
q1(t)
>>> q1d = diff(q1, Symbol('t'))
>>> vprint(q1)
q1
>>> vprint(q1d)
q1'
对于交互式会话,请使用init_vprinting。SymPy 的vprint、vpprint、latex和vlatex也存在类似的模拟。
>>> from sympy.physics.vector import init_vprinting
>>> init_vprinting(pretty_print=False)
>>> q1
q1
>>> q1d
q1'
在sympy.physics.vector中,任何时间变化的量都应使用‘dynamicsymbol’来表示,无论是坐标、变化的位置还是力量。‘dynamicsymbol’的主要用途是用于速度和坐标(在文档的运动学部分将有更多讨论)。
现在我们将用一个‘dynamicsymbol’定义我们新框架的方向,并且可以轻松地进行导数和时间导数。以下是一些例子。
>>> N = ReferenceFrame('N')
>>> B = N.orientnew('B', 'Axis', [q1, N.x])
>>> (B.y*q2 + B.z).diff(q2, N)
B.y
>>> (B.y*q2 + B.z).dt(N)
(-q1' + q2')*B.y + q2*q1'*B.z
注意,输出向量保持了它们提供时所在的框架不变。对于由多个框架的基向量组成分量的向量,这一点仍然成立:
>>> (B.y*q2 + B.z + q2*N.x).diff(q2, N)
N.x + B.y
向量的编码方式
接下来是如何在sympy.physics.vector中的代码中定义向量的简要说明。这是为了那些想要了解此部分如何工作的人提供的,不需要阅读它就可以使用此模块;除非您想了解如何实现此模块,否则不要阅读它。
每个Vector的主要信息存储在args属性中,该属性为每个相关帧中的每个基向量存储三个测量数。在创建ReferenceFrame之前,代码中不存在向量。此时,参考框架的x、y和z属性是不可变的Vector,其测量数分别为[1,0,0]、[0,1,0]和[0,0,1]。一旦可以访问这些向量,可以通过使用基向量进行代数运算来创建新向量。向量可以具有多个帧的分量。这就是为什么args是一个列表的原因;它的列表长度等于其组件中唯一ReferenceFrames的数量,即如果我们的新向量中有A和B帧基向量,则args的长度为 2;如果它有A、B和C帧基向量,则args的长度为 3。
args列表中的每个元素都是一个 2 元组;第一个元素是 SymPy Matrix(存储每组基向量的测量数),第二个元素是ReferenceFrame,用于将这些测量数与之关联。
ReferenceFrame存储了几个东西。首先,它存储您在创建时提供的名称(name属性)。它还存储方向余弦矩阵,使用orientnew方法在创建时定义,或在创建后调用orient方法。方向余弦矩阵由 SymPy 的Matrix表示,并且是一个字典的一部分,其中键是ReferenceFrame,值是Matrix;这些设置是双向的;当您将A定向到N时,您设置了A的方向字典,以包括N及其Matrix,但您也设置了N的方向字典,以包括A及其Matrix(该 DCM 是另一个的转置)。