3.3 测地线方程手推

3.3 测地线方程手推

摘要：测地线是弯曲空间中「走起来不自觉拐弯」的路线，数学身份是弧长（能量）泛函的驻定曲线，等价于切向量沿自身平行移动。本节从变分原理完整推导测地线方程，证明两种表述等价，在球面上核算大圆与纬线，最后用数值积分把测地线「走」出来，直观看到高纬纬线为何不是测地线。

老问题的新解法

「两点之间，哪条路最短？」平面上答案是线段；球面上呢——飞行航线图里那些向北拱起的弧线已经给出了经验答案：大圆劣弧。但航线图只是现象，本节要的是方程。换一个更本质的问法：如果一条曲线想成为「长度驻定」（邻近路线都不比它更短），它必须满足什么微分方程？

驻定性提示变分法。直接对弧长泛函 L = \int \sqrt{g_{ij}\dot x^i \dot x^j}\ dt 做变分可行但被开方号折磨；标准技巧是改用能量泛函 E = \tfrac{1}{2}\int g_{ij}\dot x^i \dot x^j\ dt（速度平方的积分，不开方）。由柯西不等式，L^2 \le 2E \cdot (t_1 - t_0)，等号当且仅当速度大小恒定；因此能量的驻定曲线配恒速参数化后恰是弧长的驻定曲线。参数恒速这一点今后默认遵守，称为仿射参数。

推导：从欧拉-拉格朗日方程到测地线方程

能量泛函的拉格朗日函数是 F(x, \dot x) = \tfrac{1}{2} g_{ij}(x)\dot x^i \dot x^j。套用欧拉-拉格朗日方程 \tfrac{d}{dt}\tfrac{\partial F}{\partial \dot x^k} - \tfrac{\partial F}{\partial x^k} = 0，逐项计算：

对 \dot x^k 求偏导：\tfrac{\partial F}{\partial \dot x^k} = g_{kj}\dot x^j，于是 \tfrac{d}{dt}\tfrac{\partial F}{\partial \dot x^k} = g_{kj}\ddot x^j + \partial_l g_{kj}\, \dot x^l \dot x^j；对 x^k 求偏导：\tfrac{\partial F}{\partial x^k} = \tfrac{1}{2}\partial_k g_{ij}\, \dot x^i \dot x^j。代入并移项，两边同乘 g^{mk} 解出二阶导：

这就是测地线方程。推导中「度量的导数组合恰好整理成克氏符号」一步建议亲手完成——变分法与第 2.2 节的联络公式在此会师，是本章最漂亮的两次握手之一：最短路线的方程里自然长出了搬运协议的系数，几何与动力学在此互认。

第二种表述与之等价：「曲线的切向量沿曲线自身平行移动」，即 \nabla_{\dot\gamma}\dot\gamma = 0。把协变导数展开：d\dot x^k/dt + \Gamma^k_{ij}\dot x^i\dot x^j = 0，与上式逐字相同。物理解读立现：测地线是弯曲空间中的惯性运动——不受力（无外加加速度）的粒子，其世界线自动就是测地线；所谓「引力使苹果下落」，在第六章将被改写为「地球质量弯曲时空，苹果沿测地线自由滑行」。方程右侧为零，意味着曲率并不作为「力项」出现在方程里——曲率藏在左边的 \Gamma 系数中，通过修改「什么是直线」起作用。

球面核算：大圆是测地线，纬线不是

以球面度量 ds^2 = d\theta^2 + \sin^2\theta\, d\varphi^2 与克氏符号 \Gamma^\theta_{\varphi\varphi} = -\sin\theta\cos\theta、\Gamma^\varphi_{\theta\varphi} = \cot\theta 代入方程的 \theta 分量：

大圆核验：取赤道 \theta \equiv \tfrac{\pi}{2}、\varphi = t。\dot\theta = \ddot\theta = 0，\sin\theta\cos\theta = 0，方程满足——赤道是测地线。任何大圆都是赤道的旋转像，而测地线方程在等距变换下不变，故一切大圆皆是测地线。再核 \varphi 分量 \ddot\varphi + 2\cot\theta\, \dot\theta\dot\varphi = 0：赤道上同样逐项为零。

纬线核验：取纬线 \theta \equiv \theta_0 \neq \tfrac{\pi}{2}、\varphi = t。\theta 分量给出 0 - \sin\theta_0\cos\theta_0 = 0，仅当 \theta_0 = \tfrac{\pi}{2} 成立——除赤道外的纬线都不是测地线。物理图像：飞机沿纬线平飞需要不断「向极点方向压杆」（这正是 \Gamma^\theta_{\varphi\varphi} 那一项），维持纬线要付出向内的加速度；松杆（加速度归零）后飞机自然滑入的路线就是大圆。远程航线偏爱高纬弧线，本质就是沿测地线飞才最短。

大圆与纬线的对照读图

数值实验：把测地线走出来

把方程交给数值积分器，从同一点出发、同速不同向，看球面上的测地线束如何铺开。下面的会话同时完成两件事：积分测地线方程；对比两条候选航线（大圆与纬线）的弧长。

import numpy as np
from scipy.integrate import solve_ivp

def geodesic_rhs(t, y):
    th, ph, vth, vph = y
    return [vth, vph,
            np.sin(th)*np.cos(th)*vph**2,          # theta 分量（含 Gamma 项）
            -2*np.cos(th)/np.sin(th)*vth*vph]      # varphi 分量（含 Gamma 项）

# 起点：北纬60度（余纬 pi/3）、东经0；初速指向正东（纬线切向，模长1）
y0 = [np.pi/3, 0.0, 0.0, 1.0/np.sin(np.pi/3)]
sol = solve_ivp(geodesic_rhs, [0, 2*np.pi], y0, rtol=1e-10, dense_output=True)

# 航线一：测地线飞到东经90度所需弧长（速度模恒为1，弧长=时间）
def phi_reach(t):
    return np.mod(sol.sol(t)[1], 2*np.pi)
ts = np.linspace(0, 2*np.pi, 200001)
phis = np.array([phi_reach(t) for t in ts])
t_arrive = ts[np.argmax(phis > np.pi/2)]
print(round(float(t_arrive), 4))     # 输出: 0.5456 —— 大圆弧长

# 航线二：沿北纬60度纬线走90度经度的弧长 = sin(theta0) * 弧度差
print(round(np.sin(np.pi/3)*np.pi/2, 4))   # 输出: 1.3603 —— 纬线弧长
​

两组读数与图上账单一致：同一对端点，纬线比测地线长约一倍半。改动初速方向（y0 的后两个分量）即可绘制任意出发方向的测地线；把 rhs 里的度量换成环面或伪球面，实验立即迁移——第四节将直接复用这段代码。

⚠️ 常见坑两处。坑一：忘记仿射参数，用任意参数套测地线方程。方程 \ddot x^m + \Gamma^m_{ij}\dot x^i\dot x^j = 0 只对恒速参数化成立；换参数会多出一项 f(t)\dot x^m，方程形式改变但曲线不变。坑二：把「测地线」与「最短线」画等号。测地线只是驻定曲线——大圆优弧也是测地线，却比劣弧长；测地线何时真正取到最短，要等第四章的 Hopf-Rinow 与割点理论回答。

本节要点回顾

• 推导主线：能量泛函 → 欧拉-拉格朗日 → 度量导数组合成克氏符号 → 测地线方程，全程无跳步；
• 双表述等价：变分驻定 ⇔ 切向量沿自身平行移动 \nabla_{\dot\gamma}\dot\gamma = 0；
• 球面判定：一切大圆是测地线；纬线仅赤道例外，\Gamma^\theta_{\varphi\varphi} 项是残留加速度；
• 数值实验：同一端点对，大圆弧长约为纬线的七分之三不到，航线北拱有了定量解释；
• 两个警示：方程依赖仿射参数；驻定不等于最短，最短性问题移交第四章。