【理解黎曼几何】6. 曲率的计数与计算(Python)
By 苏剑林 | 2016-10-19 | 51935位读者 |曲率的独立分量 #
黎曼曲率张量是一个非常重要的张量,当且仅当它全部分量为0时,空间才是平直的。它也出现在爱因斯坦的场方程中。总而言之,只要涉及到黎曼几何,黎曼曲率张量就必然是核心内容。
已经看到,黎曼曲率张量有4个指标,这也意味着它有$n^4$个分量,$n$是空间的维数。那么在2、3、4维空间中,它就有16、81、256个分量了,可见,要计算它,是一件相当痛苦的事情。幸好,这个张量有很多的对称性质,使得独立分量的数目大大减少,我们来分析这一点。
首先我们来导出黎曼曲率张量的一些对称性质,这部分内容是跟经典教科书是一致的。定义
$$R_{\mu\alpha\beta\gamma}=g_{\mu\nu}R^{\nu}_{\alpha\beta\gamma} \tag{50} $$
定义这个量的原因,要谈及逆变张量和协变张量的区别,我们这里主要关心几何观,因此略过对张量的详细分析。这个量被称为完全协变的黎曼曲率张量,有时候也直接叫做黎曼曲率张量,只要不至于混淆,一般不做区分。通过略微冗长的代数运算(在一般的微分几何、黎曼几何或者广义相对论教材中都有),可以得到
$$\begin{aligned}&R_{\mu\alpha\beta\gamma}=-R_{\mu\alpha\gamma\beta}\\
&R_{\mu\alpha\beta\gamma}=-R_{\alpha\mu\beta\gamma}\\
&R_{\mu\alpha\beta\gamma}=R_{\beta\gamma\mu\alpha}\\
&R_{\mu\alpha\beta\gamma}+R_{\mu\beta\gamma\alpha}+R_{\mu\gamma\alpha\beta}=0
\end{aligned} \tag{51} $$
前两个式子很快就告诉我们,如果$\beta=\gamma$或$\mu=\alpha$,那么$R_{\mu\alpha\gamma\beta}=0$,这是反对称的结果。下面的计数表明,独立分量的数目还可以以进一步减少。我们可以将每个等式看成一个约束,有一个约束,独立分量就减少1,我们要考虑有多少个独立的约束,来得出有多少个独立分量。
首先,$R_{\mu\alpha\beta\gamma}$可以看成一个$n\times n$矩阵(前两个指标),由这个矩阵的每个元素都是一个$n\times n$矩阵(后两个指标),由$R_{\mu\alpha\beta\gamma}=-R_{\mu\alpha\gamma\beta}$可知,作为矩阵每个元素的矩阵,是一个反对称矩阵,因此只有$n(n-1)/2$个独立分量,而由$R_{\mu\alpha\beta\gamma}=-R_{\alpha\mu\beta\gamma}$知道,作为矩阵的矩阵,它也是反对称的,因此也只有$n(n-1)/2$个独立分量,也就是说,根据前两个式子就可以知道总的独立分量数不超过$[n(n-1)/2]^2$。
接着,第三、第四个约束可以为我们进一步减少分量数目。$R_{\mu\alpha\beta\gamma}=R_{\beta\gamma\mu\alpha}$表明,如果将剩下的$[n(n-1)/2]^2$个分量排成$n(n-1)/2\times n(n-1)/2$的矩阵,那么它是对称的,这时候它的独立分量数目就只有
$$\frac{1}{2}\frac{n(n-1)}{2}\left[\frac{n(n-1)}{2}+1\right] \tag{52} $$
最后一个约束$R_{\mu\alpha\beta\gamma}+R_{\mu\beta\gamma\alpha}+R_{\mu\gamma\alpha\beta}=0$是针对四个不相同指标的,因为一旦有一对相同指标,就能够转化为前三个约束的线性组合。因此,独立分量的数目还需要减少$\binom{n}{4}$:
$$\frac{1}{2}\frac{n(n-1)}{2}\left[\frac{n(n-1)}{2}+1\right]-\binom{n}{4}=\frac{n^2(n^2-1)}{12} \tag{53} $$
于是,在2、3、4、5维空间中,黎曼曲率张量的独立分量数目分别是1、6、20、50。特别地,在二维空间中(即研究三维空间中的二维曲面之时),曲率张量只有一个独立分量。
曲率张量的计算 #
下述代码有误,属于SymPy的bug,已经提交给官方,请等待修正。
很遗憾,在这里我们只能略去这个曲率张量的人工计算了,对于这个曲率张量,并没有明显简单的计算方法。利用变分方法可以简化$\Gamma^{\mu}_{\alpha\beta}$的计算,进而间接简化了曲率张量的计算。而一个自上而下的简化计算方法需要用到外微分的知识,我们目前在这个系列中没有办法做到。
但在计算机方面,利用Python中的符号计算库SymPy可以帮助我们快速地计算联络系数和曲率张量。SymPy是Python的一个符号计算库,具有小巧、开源、强大等特点。当然,如果跟商业软件Mathematica相比那肯定是没法比的,但是在某些特定的方面,其实尤有过之,也就是“不胜在全、只胜在专”了。SymPy内置了专门用来处理微分几何和张量的模块,这使得它可以非常方便地计算$\Gamma_{\alpha\beta}^{\mu}$和$R_{\alpha\beta\gamma}^{\mu}$。
一个简单的代码如下:
from sympy.diffgeom import Manifold, Patch, CoordSystem
from sympy.diffgeom import TensorProduct as TP
from sympy.diffgeom import metric_to_Riemann_components as Riemann
from sympy.diffgeom import metric_to_Christoffel_2nd as Christoffel
from sympy import Symbol, Function, latex, sin
n = 2
M = Manifold('M', n)
P = Patch('P', M)
coord = CoordSystem('coord', P, ['x%s'%i for i in range(n)])
x = coord.coord_functions()
dx = coord.base_oneforms()
g = [[1, 0], [0, sin(x[0])**2]]
metric = sum([g[i][j]*TP(dx[i], dx[j]) for i in range(n) for j in range(n)])
C = Christoffel(metric)
R = Riemann(metric)
这是以二维球面坐标为例的,其中n定义了维数,g定义了度规矩阵,如果读者要修改,只需要修改这两个参数就好。里边的x变量是坐标的集合,dx变量是微分元的集合。为了一致,这里的指标是从0开始的,也就是遍历$0\sim n-1$。得到结果是
[[[0, 0], [0, -\sin(x0)\cdot \cos(x0)]], [[0, \cos(x0)/\sin(x0)], [\cos(x0)/\sin(x0), 0]]]
[[[[0, 0], [0, 0]], [[0, \sin(x0)\cdot \cdot 2], [-\sin(x0)\cdot \cdot 2, 0]]], [[[0, -1], [1, 0]], [[0, 0], [0, 0]]]]
可以看到,程序会把所有的分量都列出来,如果我们只关心非0分量,那么可以用以下代码输出:
for i1 in range(n):
for i2 in range(n):
for i3 in range(n):
if C[i1, i2, i3]:
print 'Gamma^%s_%s%s:'%(i1, i2, i3), C[i1, i2, i3]
for i1 in range(n):
for i2 in range(n):
for i3 in range(n):
for i4 in range(n):
if R[i1, i2, i3, i4]:
print 'Riemann^%s_%s%s%s:'%(i1, i2, i3, i4), R[i1, i2, i3, i4]
得到
Gamma^0_11: -\sin(x0)\cdot \cos(x0)
Gamma^1_01: \cos(x0)/\sin(x0)
Gamma^1_10: \cos(x0)/\sin(x0)
Riemann^0_101: \sin(x0)\cdot \cdot 2
Riemann^0_110: -\sin(x0)\cdot \cdot 2
Riemann^1_001: -1
Riemann^1_010: 1
注意,这里的曲率张量算的是$R_{\alpha\beta\gamma}^{\mu}$,如果读者关心$R_{\mu\alpha\beta\gamma}$,但SymPy并没有现成的函数,我们可以自己用张量运算库来算,即
from sympy.tensor.array import tensorproduct, tensorcontraction
tensorcontraction(tensorproduct(g, R), (1, 2))
其中tensorproduct是用来做张量积的,tensorcontraction是用来缩并的。得到
[[[[0, 0], [0, 0]], [[0, \sin(x0)\cdot \cdot 2], [-\sin(x0)\cdot \cdot 2, 0]]], [[[0, -\sin(x0)\cdot \cdot 2], [\sin(x0)\cdot \cdot 2, 0]], [[0, 0], [0, 0]]]]
最后如果要输出为LaTeX代码,只需要用latex命令转就行了,比如
for i1 in range(n):
for i2 in range(n):
for i3 in range(n):
if C[i1, i2, i3]:
print 'Gamma^%s_%s%s:'%(i1, i2, i3), \
latex(C[i1, i2, i3]).replace('\\mathrm{x_', '{x_').replace('\\boldsymbol{{x_', '{{x_')
最后的replace是做一些简单的微调,使得结果更符合我们平时的写法,结果是
Gamma^0_11: - \sin{\left ({{x_{0}}} \right )} \cos{\left ({{x_{0}}} \right )}
Gamma^1_01: \frac{\cos{\left ({{x_{0}}} \right )}}{\sin{\left ({{x_{0}}} \right )}}
Gamma^1_10: \frac{\cos{\left ({{x_{0}}} \right )}}{\sin{\left ({{x_{0}}} \right )}}
最后我们可能需要对含有待定函数的度量求曲率,比如$ds^2=f(x,y)(dx^2+dy^2)$,这时候大概的代码是:
f = Function('f')
g = [[f(x[0], x[1]), 0], [0, f(x[0], x[1])]]
metric = sum([g[i][j]*TP(dx[i], dx[j]) for i in range(n) for j in range(n)])
C = Christoffel(metric)
R = Riemann(metric)
这时候在$\Gamma_{\alpha\beta}^{\mu}$中,会出现这样的项:
Subs(Derivative(f(_xi_1, x1), _xi_1), (_xi_1,), (x0,))/(2\cdot f(x0, x1))
转换为LaTeX是(跟前面一样做了一下简单的处理,替换掉boldsymbol和mathrm标签)
\frac{\left. \frac{d}{d \xi_{1}} f{\left (\xi_{1},{{x_{1}}} \right )} \right|_{\substack{ \xi_{1}={{x_{0}}} }}}{2 f{\left ({{x_{0}}},{{x_{1}}} \right )}}
输出效果是
$$\frac{\left. \frac{d}{d \xi_{1}} f{\left (\xi_{1},{{x_{1}}} \right )} \right|_{\substack{ \xi_{1}={{x_{0}}} }}}{2 f{\left ({{x_{0}}},{{x_{1}}} \right )}}$$
实际上就是
$$\frac{ \frac{\partial }{\partial x_0} f\left (x_0, x_1 \right ) }{2 f{\left ({{x_{0}}},{{x_{1}}} \right )}}$$
可见,SymPy的输出还有待改善呀,很多地方还是需要人工优化一下的。更糟糕的是,$R_{\alpha\beta\gamma}^{\mu}$中出现了这样的项:
-Subs(Derivative(f(_#_1(_Dummy_38), x1), _#_1(_Dummy_38), _#_1(_Dummy_38)), (_#_1(_Dummy_38),), (x0,))/f(x0, x1) - Subs(Derivative(f(x0, _#_0(_Dummy_38)), _#_0(_Dummy_38), _#_0(_Dummy_38)), (_#_0(_Dummy_38),), (x1,))/f(x0, x1) + Subs(Derivative(f(_xi_1, x1), _xi_1), (_xi_1,), (x0,))\cdot \cdot 2/(2\cdot f(x0, x1)\cdot \cdot 2) + Subs(Derivative(f(x0, _xi_2), _xi_2), (_xi_2,), (x1,))\cdot \cdot 2/(2\cdot f(x0, x1)\cdot \cdot 2)
转换为LaTeX后是
- \frac{1}{f{\left ({{x_{0}}},{{x_{1}}} \right )}} \left. \frac{d^{2}}{d \operatorname{_#_{1}}{\left (Dummy_{38} \right )}^{2}} f{\left (\operatorname{_#_{1}}{\left (Dummy_{38} \right )},{{x_{1}}} \right )} \right|_{\substack{ \operatorname{_#_{1}}{\left (Dummy_{38} \right )}={{x_{0}}} }} - \frac{1}{f{\left ({{x_{0}}},{{x_{1}}} \right )}} \left. \frac{d^{2}}{d \operatorname{_#_{0}}{\left (Dummy_{38} \right )}^{2}} f{\left ({{x_{0}}},\operatorname{_#_{0}}{\left (Dummy_{38} \right )} \right )} \right|_{\substack{ \operatorname{_#_{0}}{\left (Dummy_{38} \right )}={{x_{1}}} }} + \frac{\left. \frac{d}{d \xi_{1}} f{\left (\xi_{1},{{x_{1}}} \right )} \right|_{\substack{ \xi_{1}={{x_{0}}} }}^{2}}{2 f^{2}{\left ({{x_{0}}},{{x_{1}}} \right )}} + \frac{\left. \frac{d}{d \xi_{2}} f{\left ({{x_{0}}},\xi_{2} \right )} \right|_{\substack{ \xi_{2}={{x_{1}}} }}^{2}}{2 f^{2}{\left ({{x_{0}}},{{x_{1}}} \right )}}
事实上这还会报错!这时候我们还需要把诸如\operatorname{_#_{0}}{\left (Dummy_{38} \right )}这些内容替换为普通的符号,比如
print latex(R[0,1,0,1]).replace('\\mathrm{x_', '{x_').replace('\\boldsymbol{{x_', '{{x_')\
.replace('\\operatorname{_#_{0}}{\\left (Dummy_{38} \\right )}', '\\xi')\
.replace('\\operatorname{_#_{1}}{\\left (Dummy_{38} \\right )}', '\\xi')
这才得到一个没报错、有点靠谱的输出:
$$\begin{aligned}&- \frac{1}{f{\left ({{x_{0}}},{{x_{1}}} \right )}} \left. \frac{d^{2}}{d \xi^{2}} f{\left (\xi,{{x_{1}}} \right )} \right|_{\substack{ \xi={{x_{0}}} }} - \frac{1}{f{\left ({{x_{0}}},{{x_{1}}} \right )}} \left. \frac{d^{2}}{d \xi^{2}} f{\left ({{x_{0}}},\xi \right )} \right|_{\substack{ \xi={{x_{1}}} }} \\
&+ \frac{\left. \frac{d}{d \xi_{1}} f{\left (\xi_{1},{{x_{1}}} \right )} \right|_{\substack{ \xi_{1}={{x_{0}}} }}^{2}}{2 f^{2}{\left ({{x_{0}}},{{x_{1}}} \right )}} + \frac{\left. \frac{d}{d \xi_{2}} f{\left ({{x_{0}}},\xi_{2} \right )} \right|_{\substack{ \xi_{2}={{x_{1}}} }}^{2}}{2 f^{2}{\left ({{x_{0}}},{{x_{1}}} \right )}}\end{aligned}$$
最后整理为
$$-\frac{\frac{\partial^2 f(x_0, x_1)}{\partial x_0^2}}{f(x_0,x_1)} -\frac{\frac{\partial^2 f(x_0, x_1)}{\partial x_1^2}}{f(x_0,x_1)} + \frac{\left(\frac{\partial f(x_0, x_1)}{\partial x_0}\right)^2}{2f^2 (x_0,x_1)} + \frac{\left(\frac{\partial f(x_0, x_1)}{\partial x_1}\right)^2}{2f^2 (x_0,x_1)}$$
或者
$$-\frac{\nabla^2 f}{f}+\frac{1}{2}\left|\frac{\nabla f}{f}\right|^2$$
转载到请包括本文地址:https://kexue.fm/archives/4026
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Oct. 19, 2016). 《【理解黎曼几何】6. 曲率的计数与计算(Python) 》[Blog post]. Retrieved from https://kexue.fm/archives/4026
@online{kexuefm-4026,
title={【理解黎曼几何】6. 曲率的计数与计算(Python)},
author={苏剑林},
year={2016},
month={Oct},
url={\url{https://kexue.fm/archives/4026}},
}
October 23rd, 2016
(53)中,等号左边的第一个式子还得乘一个1/2 ^-^
October 23rd, 2016
或者说,(52)式差一个1/2
你是对的,已经修正。感谢阅读并指出^_^