02.4 測地線：曲面上的「直線」#

本節目標#

給出「直線」的內在定義：測地線 = 局部最短（或至少「不轉彎」）的曲線。
從第一基本形式 \(E, F, G\) 寫出測地線方程（Christoffel 符號版本）。
用 SymPy 驗證三個經典事實：平面上測地線是直線、球面上是大圓、柱面上是螺旋線。
實驗：在球面上畫出台北到舊金山的大圓航線，解釋為什麼飛機要飛「北邊」。

幾何直覺：什麼是曲面上的直線？#

在平面上，直線有兩個等價刻劃：

最短：連接兩點的所有曲線中，直線長度最短。
不轉彎：直線的加速度向量為零；在曲面上，「加速度」只能談切平面裡的分量，所以我們要求加速度的切向分量為零——法向分量被曲面強迫，我們管不著，但小蟲自己「感覺不到轉彎」就夠了。

刻劃 2 比較好算，就叫它測地線：曲線 \(\gamma(s) = \mathbf{r}(u(s), v(s))\) （ \(s\) 為弧長）是測地線，若

\[ \gamma'' \perp T_p S \quad \text{（二階導完全落在法向方向）}. \]

刻劃 1 在「足夠短的區間」上與它等價（全域最短可能不成立：球面上兩點間的大圓有優弧與劣弧兩段，優弧是測地線卻不是最短——所以說局部最短）。

測地線方程#

把條件 \(\gamma''\) 無切向分量寫成座標式。令 \(q^1 = u,\ q^2 = v\) ，度量矩陣 \(g = \begin{pmatrix} E & F \\ F & G\end{pmatrix}\) ，定義 Christoffel 符號

\[ \Gamma^i_{jk} = \frac{1}{2} g^{il} \left( \frac{\partial g_{kl}}{\partial q^j} + \frac{\partial g_{jl}}{\partial q^k} - \frac{\partial g_{jk}}{\partial q^l} \right), \]

則測地線方程是兩條二階 ODE：

\[ \frac{d^2 q^i}{ds^2} + \Gamma^i_{jk}\, \frac{dq^j}{ds}\, \frac{dq^k}{ds} = 0, \qquad i = 1, 2. \]

注意方程只用到 \(E, F, G\) ——測地線和 02.3 的 \(K\) 一樣是內在概念：住在曲面上的二維生物也能定義並找出牠的「直線」。

SymPy 實作：由第一基本形式產生測地線方程#

import sympy as sp

s = sp.symbols('s', real=True)
u1, u2 = sp.symbols('u1 u2', real=True)          # 座標 q1, q2
q = [sp.Function('q1')(s), sp.Function('q2')(s)] # 待解的曲線

def christoffel(E, F, G):
    """回傳 Gamma[i][j][k]（i,j,k ∈ {0,1}），由 E,F,G 組成。
    約定：E 是 du1^2 的係數，G 是 du2^2 的係數。"""
    g    = sp.Matrix([[E, F], [F, G]])
    ginv = g.inv()
    Gam = [[[0]*2 for _ in range(2)] for _ in range(2)]
    for i in range(2):
        for j in range(2):
            for k in range(2):
                term = 0
                for l in range(2):
                    term += ginv[i, l]*(sp.diff(g[k, l], (u1, u2)[j])
                                      + sp.diff(g[j, l], (u1, u2)[k])
                                      - sp.diff(g[j, k], (u1, u2)[l]))
                Gam[i][j][k] = sp.simplify(term/2)
    return Gam

def geodesic_eqs(E, F, G):
    """回傳兩條測地線方程（皆應為 0）。"""
    Gam = christoffel(E, F, G)
    eqs = []
    for i in range(2):
        eq = sp.diff(q[i], s, 2)
        for j in range(2):
            for k in range(2):
                eq += Gam[i][j][k]*sp.diff(q[j], s)*sp.diff(q[k], s)
        eqs.append(sp.simplify(eq))
    return eqs

驗證一：平面 → 直線#

平面 \(\mathbf{r} = (u, v, 0)\) 的度量是 \(E = G = 1,\ F = 0\) ，所有 \(\Gamma = 0\) ，方程退化為 \(q_1'' = q_2'' = 0\) ，通解顯然是直線：

a1, a2, b1, b2 = sp.symbols('a1 a2 b1 b2')
eqs = geodesic_eqs(1, 0, 1)
print(eqs)  # [Derivative(q1(s), (s, 2)), Derivative(q2(s), (s, 2))]
sol = [a1*s + b1, a2*s + b2]        # 候選解：等速直線
check = [sp.simplify(e.subs({q[i]: sol[i] for i in range(2)})) for e in eqs]
print(check)                        # [0, 0]  → 直線是測地線 ✓

驗證二：球面 → 大圓#

球面（半徑 \(R\) ， \(u\) 為經度、 \(t\) 為餘緯）的度量是 \(ds^2 = R^2 \sin^2 t\, du^2 + R^2\, dt^2\) 。以 \(u_1 = u\) 、 \(u_2 = t\) 代入：

R = sp.symbols('R', positive=True)
Es = R**2*sp.sin(u2)**2   # du^2 的係數
Gs = R**2                 # dt^2 的係數
eqs = geodesic_eqs(Es, 0, Gs)
# 候選解一：赤道 t = pi/2（常數），u = s/R（等速）
rep1 = {q[0]: s/R, q[1]: sp.pi/2 + 0*s, u1: 0, u2: sp.pi/2}  # 座標值也要代入
print([sp.simplify(e.subs(rep1)) for e in eqs])   # [0, 0]
# 候選解二：經線 u = 常數，t = s/R
rep2 = {q[0]: 0*s, q[1]: s/R}
print([sp.simplify(e.subs(rep2)) for e in eqs])   # [0, 0]

赤道與經線都是大圓，確實是測地線。（一般的大圓可由球面對稱旋轉得到：測地方程是內在的，旋轉不改變它。）

驗證三：柱面 → 螺旋線#

柱面 \(\mathbf{r} = (R\cos u, R\sin u, v)\) 的度量 \(E = R^2,\ G = 1,\ F = 0\) ， Christoffel 也全為 \(0\) （係數是常數），方程又是 \(q'' = 0\) ，解為 \(u = a s,\ v = b s\) ： \(b = 0\) 時是水平圓、 \(a = 0\) 時是豎直母線、兩者皆非零就是螺旋線。

print(geodesic_eqs(R**2, 0, 1))      # 又是純二階導 → 螺旋線 u = a s, v = b s

有趣對照：柱面外在看起來彎，但內在（ \(K = 0\) ，見 02.3）是平的，所以牠的直線攤平後就是平面直線——螺旋線攤平正好是一條斜直線。

畫圖：三種測地線#

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(15, 4))
ss = np.linspace(0, 4*np.pi, 200)

ax = fig.add_subplot(131, projection='3d')          # 平面：直線
ax.plot(ss, 0.6*ss, 0*ss, 'r', lw=2)
ax.set_title('plane: straight line')

ax = fig.add_subplot(132, projection='3d')          # 球面：大圓
th = np.linspace(0, np.pi, 60)
TH, PH = np.meshgrid(th, np.linspace(0, 2*np.pi, 60))
ax.plot_surface(np.sin(TH)*np.cos(PH), np.sin(TH)*np.sin(PH), np.cos(TH),
                alpha=0.15, color='gray')
ang = np.linspace(0, 2*np.pi, 200)
gc = np.stack([np.cos(ang), 0.6*np.sin(ang), 0.8*np.sin(ang)])
gc = gc/np.linalg.norm(gc, axis=0)   # 大圓 = 球面與過球心平面的交線
ax.plot(*gc, 'r', lw=2)
ax.set_title('sphere: great circle')

ax = fig.add_subplot(133, projection='3d')          # 柱面：螺旋線
ax.plot(np.cos(ss), np.sin(ss), 0.5*ss, 'r', lw=2)
ax.set_title('cylinder: helix')
plt.tight_layout(); plt.show()

動手實驗：台北 → 舊金山的大圓#

台北（TPE，北緯 \(25.03^\circ\) ，東經 \(121.57^\circ\) ）與舊金山（SFO，北緯 \(37.77^\circ\) ，西經 \(122.42^\circ\) ）。飛機走的「最短線」是過這兩點與地心的平面切出的大圓——它會明顯偏向阿留申群島北方，而不是沿緯線走。

import numpy as np
import matplotlib.pyplot as plt

R = 1.0
def ll_to_xyz(lat_deg, lon_deg):
    la, lo = np.radians(lat_deg), np.radians(lon_deg)
    return np.array([R*np.cos(la)*np.cos(lo),
                     R*np.cos(la)*np.sin(lo),
                     R*np.sin(la)])

TPE = ll_to_xyz( 25.03,  121.57)
SFO = ll_to_xyz( 37.77, -122.42)

# 大圓參數化：n = 兩點張成平面（過地心）的法向量
n = np.cross(TPE, SFO); n = n/np.linalg.norm(n)
# 從 TPE 沿大圓轉到 SFO：夾角 omega，轉軸 n
omega = np.arccos(np.clip(np.dot(TPE, SFO)/(R*R), -1, 1))
ang = np.linspace(0, omega, 100)
# Rodrigues 公式：TPE 繞 n 旋轉 ang
path = np.array([TPE*np.cos(a)
                 + np.cross(n, TPE)*np.sin(a)
                 + n*np.dot(n, TPE)*(1 - np.cos(a)) for a in ang])

# 對照組：沿「等緯度」走的路線（把兩點經度連起來、緯度線性內插）
lats = np.linspace(25.03, 37.77, 100)
lons = np.linspace(121.57, -122.42, 100)
rhumb = np.array([ll_to_xyz(la, lo) for la, lo in zip(lats, lons)])

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(projection='3d')
th = np.linspace(-np.pi/2, np.pi/2, 40)
ph = np.linspace(-np.pi, np.pi, 80)
TH, PH = np.meshgrid(th, ph)
ax.plot_surface(R*np.cos(TH)*np.cos(PH), R*np.cos(TH)*np.sin(PH), R*np.sin(TH),
                alpha=0.15, color='gray')
ax.plot(*path.T,  'r', lw=3, label='great circle (geodesic)')
ax.plot(*rhumb.T, 'b--', lw=2, label='along latitude (longer!)')
ax.scatter(*TPE, color='k', s=40); ax.text(*TPE*1.1, ' TPE')
ax.scatter(*SFO, color='g', s=40); ax.text(*SFO*1.1, ' SFO')
ax.legend(); ax.set_box_aspect((1, 1, 1)); plt.show()

# 量長度（單位：地球半徑倍數 → 乘 6371 km）
L_gc    = omega*R
L_rhumb = np.sum(np.linalg.norm(np.diff(rhumb, axis=0), axis=1))
print(L_gc, L_rhumb)   # 大圓較短；乘上 6371 km 可得實際里程

紅色大圓明顯「凸向北方」——這正是 02.3 的剪刀效應在導航上的化身：沿緯線走看似「直」，但在 \(K > 0\) 的球面幾何裡它不是測地線。

小結與練習#

測地線 = 加速度無切向分量 = 局部最短；方程 \(q'' + \Gamma^i_{jk} q'^j q'^k = 0\) 完全由第一基本形式決定（內在！）。
平面直線、球面大圓、柱面螺旋線是同一條方程的三種面貌；度量常數為零時測地線就是等速直線的座標像。
全域最短要注意「優弧陷阱」：測地線未必最短，最短的測地線只是局部最短。

練習：

把平面測地線方程 \(q'' = 0\) 的通解 \(q_i = a_i s + b_i\) 代回，確認 SymPy 給 \([0, 0]\) ；再證明弧長參數下 \(a_1^2 + a_2^2 = 1\) （ \(s\) 是弧長）。
球面上除了赤道與經線，驗證一般大圓也是解：取任意過原點平面的單位法向量 \(\mathbf{n} = (n_1, n_2, n_3)\) ，寫出 \(\gamma(s)\) 的參數式代入 geodesic_eqs（提示：可能需要 sp.trigsimp 加上 \(\sin^2+\cos^2 = 1\) ）。
用 sp.dsolve 嘗試解平面與柱面的測地線方程，觀察輸出。
（挑戰）計算圓錐面 \(\mathbf{r} = (v\cos u, v\sin u, av)\) 的測地線方程；它有個守恆量 \(v^2 u'\) （對應角動量）。用 sp.checkodesol 或數值驗證。
把台北→東京（北緯 \(35.68^\circ\) 、東經 \(139.69^\circ\) ）加入上圖，比較大圓與沿緯線的長度差，估算省下的公里數。