Understanding Riemannian Geometry: 6. Counting and Computing Curvature (Python)

【理解黎曼几何】6. 曲率的计数与计算(Python) 苏剑林 · 2016-10-19 · Mathematics geometry curvature Python riemannian geometry permutation and combination

Independent Components of Curvature

Independent components of curvature

The Riemann curvature tensor is an extremely important tensor: a space is flat if and only if all of its components vanish. It also appears in Einstein's field equations. In short, whenever Riemannian geometry is involved, the Riemann curvature tensor is inevitably at the heart of things.

We've already seen that the Riemann curvature tensor has 4 indices, which means it has $n^4$ components, where $n$ is the dimension of the space. So in 2, 3, and 4 dimensions it has 16, 81, and 256 components respectively — clearly, computing it by hand would be quite painful. Fortunately, this tensor has a great deal of symmetry, which drastically reduces the number of independent components. Let's work this out.

First let's derive some symmetry properties of the Riemann curvature tensor; this part is consistent with what you'll find in classical textbooks. Define

$$R_{\mu\alpha\beta\gamma}=g_{\mu\nu}R^{\nu}_{\alpha\beta\gamma} \tag{50} $$

The reason for defining this quantity involves the distinction between contravariant and covariant tensors; since we're mainly concerned with the geometric picture here, we'll skip the detailed tensor analysis. This quantity is called the fully covariant Riemann curvature tensor, and is sometimes simply called the Riemann curvature tensor — as long as no confusion arises, we generally won't distinguish between the two. Through a somewhat tedious algebraic derivation (found in any standard textbook on differential geometry, Riemannian geometry, or general relativity), one obtains

$$\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} $$more

The first two equations immediately tell us that if $\beta=\gamma$ or $\mu=\alpha$, then $R_{\mu\alpha\gamma\beta}=0$ — this is the antisymmetry result. The counting below shows that the number of independent components can be reduced even further. We can treat each equation as a constraint: each constraint reduces the number of independent components by (at least) one, and we want to figure out how many independent constraints there are in order to determine the number of independent components.

First, $R_{\mu\alpha\beta\gamma}$ can be viewed as a $n\times n$ matrix (indexed by the first two indices), where every entry of this matrix is itself a $n\times n$ matrix (indexed by the last two indices). From $R_{\mu\alpha\beta\gamma}=-R_{\mu\alpha\gamma\beta}$ we know that, as a matrix, each entry-matrix is antisymmetric, so it has only $n(n-1)/2$ independent components; and from $R_{\mu\alpha\beta\gamma}=-R_{\alpha\mu\beta\gamma}$ we know that the outer matrix (of matrices) is also antisymmetric, so it too has only $n(n-1)/2$ independent components. That is to say, based on the first two equations alone, the total number of independent components is at most $[n(n-1)/2]^2$.

Next, the third and fourth constraints let us reduce the component count further. $R_{\mu\alpha\beta\gamma}=R_{\beta\gamma\mu\alpha}$ shows that if we arrange the remaining $[n(n-1)/2]^2$ components into a $n(n-1)/2\times n(n-1)/2$ matrix, this matrix is symmetric, so the number of independent components becomes just

$$\frac{1}{2}\frac{n(n-1)}{2}\left[\frac{n(n-1)}{2}+1\right] \tag{52} $$

The last constraint, $R_{\mu\alpha\beta\gamma}+R_{\mu\beta\gamma\alpha}+R_{\mu\gamma\alpha\beta}=0$, applies to four mutually distinct indices, since as soon as a pair of indices coincide, it reduces to a linear combination of the first three constraints. So the number of independent components must be reduced by a further $\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} $$

Thus, in 2, 3, 4, and 5 dimensions, the number of independent components of the Riemann curvature tensor is 1, 6, 20, and 50 respectively. In particular, in two dimensions (i.e., when studying a 2D surface embedded in 3D space), the curvature tensor has only a single independent component.

Computing the curvature tensor

The code below has an issue due to a SymPy bug, which has already been reported upstream — please wait for a fix.

Unfortunately, we'll have to skip the manual computation of the curvature tensor here, since there's no obviously simple way to compute it by hand. Variational methods can simplify the computation of $\Gamma^{\mu}_{\alpha\beta}$, which indirectly simplifies the computation of the curvature tensor. A top-down simplified computation method requires knowledge of exterior differentiation, which we're not yet in a position to cover in this series.

On the computational side, however, we can use Python's symbolic computation library SymPy to quickly compute the connection coefficients and the curvature tensor. SymPy is a symbolic computation library for Python, known for being lightweight, open-source, and powerful. Of course, it's no match for commercial software like Mathematica overall, but in certain specific respects it actually holds its own — it doesn't try to win on breadth, only on depth in particular areas. SymPy has a built-in module specifically for differential geometry and tensors, which makes it very convenient to compute $\Gamma_{\alpha\beta}^{\mu}$ and $R_{\alpha\beta\gamma}^{\mu}$.

Here's a simple piece of code:

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)

This uses spherical coordinates on a 2-sphere as an example, where n defines the dimension and g defines the metric matrix; if you want to modify it, you only need to change these two parameters. The variable x is the set of coordinates, and dx is the set of differential elements. For consistency, indices here start from 0, i.e., they range over $0\sim n-1$. The result is

[[[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]]]]

As you can see, the program will list all components; if we only care about the nonzero ones, we can use the following code to print them out:

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]

which gives

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

Note that the curvature tensor computed here is $R_{\alpha\beta\gamma}^{\mu}$; if the reader is instead interested in $R_{\mu\alpha\beta\gamma}$, SymPy doesn't have a ready-made function for that, but we can compute it ourselves using the tensor operations library, i.e.,

from sympy.tensor.array import tensorproduct, tensorcontraction
tensorcontraction(tensorproduct(g, R), (1, 2))

where tensorproduct performs the tensor product and tensorcontraction performs the contraction. This gives

[[[[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]]]]

Finally, if you want to output LaTeX code, you just need to call the latex command to convert it, for example

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_')

The final replace calls do a bit of simple touch-up, to make the result match our usual notation more closely; the result is

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 )}}

Finally, we may sometimes need to compute the curvature of a metric containing an undetermined function, such as $ds^2=f(x,y)(dx^2+dy^2)$; in that case, the code looks roughly like this:

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)

At this point, $\Gamma_{\alpha\beta}^{\mu}$ will contain a term like this:

Subs(Derivative(f(_xi_1, x1), _xi_1), (_xi_1,), (x0,))/(2\cdot f(x0, x1))

Converted to LaTeX, this is (again with the same minor cleanup, replacing the boldsymbol and mathrm tags):

\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 )}}

which renders as

$$\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 )}}$$

and is, in effect, just

$$\frac{ \frac{\partial }{\partial x_0} f\left (x_0, x_1 \right ) }{2 f{\left ({{x_{0}}},{{x_{1}}} \right )}}$$

As you can see, SymPy's output still leaves something to be desired — a fair amount of manual cleanup is still needed. Worse still, $R_{\alpha\beta\gamma}^{\mu}$ contains a term like this:

-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)

Converted to LaTeX, this becomes

- \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 )}}

In fact, this even throws an error! At this point we also need to replace things like \operatorname{_#_{0}}{\left (Dummy_{38} \right )} with ordinary symbols, for example

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')

which finally gives us an output that doesn't error out and is at least somewhat trustworthy:

$$\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}$$

Finally, tidying this up, we get

$$-\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)}$$

or equivalently

$$-\frac{\nabla^2 f}{f}+\frac{1}{2}\left|\frac{\nabla f}{f}\right|^2$$