07.4 平行移動與測地線#

本節目標#

理解平行移動的定義 \(\nabla_{\dot\gamma}\dot\gamma = 0\) ：向量沿曲線「保持自己」的精確說法。
建立「測地線 = 自己與自己平行 = 直線的推廣」的幾何直覺。
親手做球面上的平行移動實驗：向量沿赤道繞一圈回到起點，轉了一個角度——這是曲率的預兆（Gauss–Bonnet 的伏筆）。
用 SymPy 求解測地線方程 \(u'' + \Gamma u'u' = 0\) ，驗證球面的大圓是測地線。
用數值 ODE solver（scipy.integrate.solve_ivp）在一般曲面上求測地線並畫圖。

平行移動：向量沿線「不變化」#

在 \(\mathbb{R}^n\) 中，說一個向量沿曲線 \(\gamma(t)\) 保持不變，就是每個分量都是常數： \(\dot v^i = 0\) 。但在曲面或彎曲座標上，這個說法失效了——回憶 07.2：對分量求導 \(\dot v^i\) 會漏掉基底本身的變化。正確的「變化率」是協變導數沿曲線的版本 \(\nabla_{\dot\gamma} v = \dot v^i + \Gamma^i_{jk}\, \dot\gamma^j v^k\) （第 \(i\) 分量； \(\Gamma^i_{jk}\) 在曲線上的點取值）。於是平行移動的定義：

\[ \nabla_{\dot\gamma} v = 0 \quad\Longleftrightarrow\quad \dot v^i + \Gamma^i_{jk}\, \dot\gamma^j v^k = 0 \ \text{（對每個 } i\text{）} \]

這是一條一階線性 ODE：給定曲線 \(\gamma\) 與初始向量 \(v(0)\) ，沿線的 \(v(t)\) 被唯一決定。由度量相容性 \(\nabla g = 0\) 得三個推論：保持長度 \(\frac{d}{dt}|v|^2 = \nabla_{\dot\gamma}(g(v,v)) = 0\) ；兩個向量一起平行移動時保持內積（夾角不變）；以及最重要的——平行移動依賴於路徑，這正是彎曲空間的靈魂所在。

幾何直覺與曲率預兆：球面繞圈實驗#

在平面上，沿任何閉曲線平行移動一圈，向量回到原方向。但在球面上：在赤道上取一個指向北極的向量，沿赤道平行移動一圈——赤道上 \(\Gamma^\phi_{\theta\phi} = \cot\theta\) 在 \(\theta = \pi/2\) 處為零，向量保持分量不變即可，回到起點時仍然指向北。

這看起來沒事？換一條路徑：先沿一條經線走到北極，再沿另一條經線走下來，兩經線夾角 \(\alpha\) 。因為經線本身是測地線（切向量沿自身平行移動），三角形繞一圈後向量**轉了角度 \(\alpha\) **。一般地，球面上沿閉曲線平行移動一圈的轉角等於圍住面積乘上高斯曲率 \(K = 1/r^2\) ，即 \(\Delta\alpha = \iint K\, dA\) （Gauss–Bonnet 的特例）。

這是曲率的預兆：平行移動的路徑依賴性，正是「空間彎曲」的可操作定義。第 8 章的 Riemann 張量，就是把「繞無窮小迴圈轉了角度」的效應精確編碼的機器。

測地線：自己與自己平行#

曲線 \(\gamma\) 若自己的切向量沿自己平行移動， \(\nabla_{\dot\gamma}\dot\gamma = 0\) ，就稱為測地線。直覺：這是「不轉彎、不加速」的曲線——自由落體、橡皮筋繃直、飛機大圓航線，都是測地線。寫成分量方程（ \(\gamma\) 的第 \(i\) 個座標記為 \(u^i(t)\) ）：

\[ \ddot u^i + \Gamma^i_{jk}\, \dot u^j \dot u^k = 0 \]

這是二階非線性 ODE 系統。注意結構： \(\ddot u^i\) 被「Christoffel 產生的假力」抵消，才能保持直行。兩個重要性質：等速性 \(\frac{d}{dt} g(\dot\gamma, \dot\gamma) = 0\) （初始速度的長度決定參數化）；變分觀點——測地線是長度（能量）泛函的臨界曲線，兩點之間取極值者。

SymPy 求解測地線方程：球面大圓驗證#

球面 \((\theta, \phi)\) 座標，度量 \(g = \mathrm{diag}(1, \sin^2\theta)\) ，非零 Christoffel 符號（07.2 算過）為 \(\Gamma^\theta_{\phi\phi} = -\sin\theta\cos\theta\) 、 \(\Gamma^\phi_{\theta\phi} = \cot\theta\) ，代入測地線方程：

\[ \ddot\theta - \sin\theta\cos\theta\, \dot\phi^2 = 0, \qquad \ddot\phi + 2\cot\theta\, \dot\theta\dot\phi = 0 \]

用 SymPy 的 dsolve 驗證：赤道（ \(\theta \equiv \pi/2\) ， \(\phi = t\) ）與經線（ \(\phi \equiv\) 常數）都滿足方程；一般大圓也滿足。

import sympy as sp

t = sp.symbols('t')
th, phi = sp.Function('theta')(t), sp.Function('phi')(t)
# 測地線方程（單位球面）
eq1 = sp.diff(th, t, 2) - sp.sin(th)*sp.cos(th)*sp.diff(phi, t)**2
eq2 = sp.diff(phi, t, 2) + 2*sp.cos(th)/sp.sin(th)*sp.diff(th, t)*sp.diff(phi, t)
# 驗證 1：赤道 θ = π/2, φ = t；驗證 2：經線 φ = c, θ = t
print("赤道:", sp.simplify(eq1.subs({th: sp.pi/2, phi: t})),
      sp.simplify(eq2.subs({th: sp.pi/2, phi: t})))
print("經線:", sp.simplify(eq1.subs({th: t, phi: sp.Symbol('c')})),
      sp.simplify(eq2.subs({th: t, phi: sp.Symbol('c')})))

輸出全部為 0——赤道與經線（都是大圓）是測地線。緯線不是：把 \(\theta(t) = \theta_0 \neq \pi/2\) 、 \(\phi = t\) 代入，eq1 = -sin θ₀ cos θ₀ ≠ 0——緯線上走需要「轉彎力」，這正是縮尺地球儀上飛機航線向極點彎的原因。

數值法：用 ODE solver 在曲面上求測地線#

一般曲面上 dsolve 無能為力，改用 scipy.integrate.solve_ivp：把二階系統化成一階系統（狀態：位置 + 速度），Christoffel 符號先用 SymPy 算好、再 lambdify 成數值函數。以單位球為例，解測地線並畫在 3D 球面上：

import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# --- SymPy 部分：符號算出 Christoffel，再轉數值函數 ---
s_th, s_ph = sp.symbols('theta phi', real=True)
g = sp.Matrix([[1, 0], [0, sp.sin(s_th)**2]])
ginv, coords = g.inv(), [s_th, s_ph]
Gamma_sym = [[[(sp.Rational(1,2)*sum(ginv[i,m]*(sp.diff(g[m,k], coords[j])
               + sp.diff(g[m,j], coords[k]) - sp.diff(g[j,k], coords[m]))
               for m in range(2))).simplify()
               for k in range(2)] for j in range(2)] for i in range(2)]
Gamma_num = sp.lambdify((s_th, s_ph), sp.Array(Gamma_sym), 'numpy')

# --- 數值部分：一階 ODE 系統 y = [θ, φ, θ̇, φ̇] ---
def geodesic_rhs(t, y):
    th_, ph_, dth, dph = y
    G = Gamma_num(th_, ph_)
    return [dth, dph,
            -G[0][0][0]*dth*dth - 2*G[0][0][1]*dth*dph - G[0][1][1]*dph*dph,
            -G[1][0][0]*dth*dth - 2*G[1][0][1]*dth*dph - G[1][1][1]*dph*dph]

sol = solve_ivp(geodesic_rhs, [0, 2*np.pi],
                [np.pi/3, 0.0, 0.0, 1.0], rtol=1e-10, atol=1e-12)

# --- 3D 畫圖：球面 + 測地線 ---
u, v = np.meshgrid(np.linspace(0, np.pi, 40), np.linspace(0, 2*np.pi, 80))
xs, ys, zs = np.sin(u)*np.cos(v), np.sin(u)*np.sin(v), np.cos(u)
gx, gy, gz = (np.sin(sol.y[0])*np.cos(sol.y[1]),
              np.sin(sol.y[0])*np.sin(sol.y[1]), np.cos(sol.y[0]))
fig = plt.figure(figsize=(7, 6))
ax = fig.add_subplot(projection='3d')
ax.plot_surface(xs, ys, zs, alpha=0.15, color='steelblue')
ax.plot(gx, gy, gz, 'r-', lw=2)
ax.set_title('Geodesic on the unit sphere (numerical)')
plt.show()

觀察：紅色曲線是一條大圓——數值解從微分方程自動「長出」大圓，程式裡從頭到尾沒有大圓的任何公式。（順帶檢查： \(g_{ij}\dot u^i\dot u^j\) 沿解為常數——等速性。）

動手實驗#

實驗 1：球面上一條一般測地線的數值解（主實驗）

把上面的初始條件換成斜的方向，例如 y0 = [np.pi/3, np.pi/4, 0.7, 0.9]，積分時間拉長到 \(4\pi\) （繞球兩圈的時間尺度）：

y0 = [np.pi/3, np.pi/4, 0.7, 0.9]
sol = solve_ivp(geodesic_rhs, [0, 4*np.pi], y0, rtol=1e-11, atol=1e-13)
xyz = np.stack([np.sin(sol.y[0])*np.cos(sol.y[1]),
                np.sin(sol.y[0])*np.sin(sol.y[1]), np.cos(sol.y[0])])
d = np.linalg.norm(xyz - xyz[:, :1], axis=0)
print("起點:", sol.y[:2, 0], " 回歸距離:", d[d > 1e-3].min())  # 應很小：閉合大圓

確認：不管初始方向為何，球面上的測地線永遠是閉合大圓；再用 \(g_{ij}\dot u^i\dot u^j\) 驗證等速性（沿解為常數）。

實驗 2：平行移動繞赤道一圈，向量轉了角度嗎？

沿緯線 \(\theta = \theta_0\) （非大圓）平行移動一個向量一圈，數值解 \(\dot v^i = -\Gamma^i_{jk}\dot\gamma^j v^k\) ，比較起終向量的夾角：

th0 = np.pi/3   # 緯度 60° 的緯線
def rhs(t, y):  # y = [φ, vθ, vφ]；γ̇ = (0, 1)
    G = Gamma_num(th0, y[0])
    return [1.0, -G[0][1][1]*y[2], -G[1][1][0]*y[1]]
sol = solve_ivp(rhs, [0, 2*np.pi], [0.0, 0.0, 1.0], rtol=1e-12, atol=1e-14)
v_end, v_start = sol.y[1:, -1], sol.y[1:, 0]
s2 = np.sin(th0)**2   # 度量沿緯線為常數：內積 = vθ wθ + sin²θ0 vφ wφ
num = v_end[0]*v_start[0] + s2*v_end[1]*v_start[1]
den = np.sqrt(v_start[0]**2 + s2*v_start[1]**2) * np.sqrt(v_end[0]**2 + s2*v_end[1]**2)
cosang = num/den      # 注意：長度要用度量算，不是 np.linalg.norm！
print("轉角（弧度）:", np.arccos(np.clip(cosang, -1, 1)))
print("理論值 2π(1-cos θ0) =", 2*np.pi*(1-np.cos(th0)))  # 和樂 = 面積 × K

兩者應相符（ \(\theta_0 = \pi/3\) 時恰好為 \(\pi\) ）；換 \(\theta_0 = \pi/6\) 再驗一次公式的一般性。這就是holonomy（和樂）：轉角 = 圍住的球冠面積 × \(K\) ——第 8 章 Riemann 張量要精確描述的效應。

實驗 3：測地線偏離的加速（曲率的另一種味道）

在球面上取兩條初始平行的測地線（起點相差 \(\epsilon\) 、初速相同），數值追蹤距離 \(s(t)\) ：先收斂（對蹠點附近相交）再發散——正曲率的特徵。畫 \(s(t)\) ，與 \(s(t) = \epsilon\cos t\) 對照。

實驗 4（挑戰）：換一個曲面

把 g 換成嵌入甜甜圈的度量 \(g = \mathrm{diag}(b^2,\ (a+b\cos u)^2)\) ，重跑同一 pipeline：測地線在外側（ \(K>0\) ）向內彎、內側（ \(K<0\) ）散開。只換度量、不改架構——座標不變性的紅利。

小結與練習#

本節重點

平行移動 \(\nabla_{\dot\gamma} v = 0\) 是「向量沿線不變化」的精確版本，由 Christoffel 符號與曲線共同決定。
平行移動保持長度與內積，但依賴路徑：球面上繞一圈向量轉角 \(\Delta\alpha = \iint K\,dA\) ——曲率的預兆。
測地線 \(\nabla_{\dot\gamma}\dot\gamma = 0\) 即 \(\ddot u^i + \Gamma^i_{jk}\dot u^j\dot u^k = 0\) ；等速；是長度泛函的臨界線。
球面的測地線恰為大圓（赤道、經線可符號驗證；緯線不是）。
符號（SymPy 算 \(\Gamma\) ）＋數值（solve_ivp 積分）的混合 pipeline 可處理任意曲面，只換度量即可。

練習

用 SymPy 把一般大圓 \(\theta(t), \phi(t)\) 的參數式（由 \(n\cdot x = 0\) 參數化）代入測地線方程，驗證恆為零。
推導緯線 \(\theta = \theta_0\) 的 \(\ddot\phi = 0\) 、 \(\ddot\theta = \sin\theta_0\cos\theta_0\,\dot\phi^2\) ，說明緯線非測地線；求緯線的測地曲率 \(|\cot\theta_0|\) 並用 SymPy 驗證。
把實驗 1 的初速長度改為 2，確認路徑形狀不變、只是「時間變快兩倍」——驗證測地線的仿射參數化自由度。
用 scipy.optimize.minimize 對離散折線的總長度做數值最佳化（兩端固定在球面上兩點），確認收斂結果貼合大圓——變分觀點的數值驗證。
（挑戰）旋轉曲面 \(X(u,v) = (f(u)\cos v, f(u)\sin v, g(u))\) 的測地線有 Clairaut 定理 \(f(u)\sin\alpha =\) 常數。用數值解驗證：在球面上此常數即大圓到某軸的距離。