跳轉到

01.3 平面與空間曲率計算#

本節目標#

  • 從定義出發推導平面曲線的經典曲率公式 \(\kappa = \dfrac{|y''|}{(1+y'^2)^{3/2}}\) 。
  • 用 SymPy 實作此公式（diff、norm/sqrt、simplify），並處理化簡技巧。
  • 掌握空間曲線曲率的一般公式 \(\kappa = \|\gamma'\times\gamma''\|/\|\gamma'\|^3\) 與 SymPy 的 cross 實作。
  • 實驗：拋物線、圓、螺旋線的曲率數值觀察，建立對曲率的「手感」。

幾何直覺：曲率到底在量什麼#

站在一條路上開車。方向盤轉得越兇、車速越快，你感受到的離心力越大。「曲率」就是把「方向盤轉得多兇」獨立出來的量——它不依賴車速，只依賴路本身。

用密切圓的語言：在曲線一點 \(P\) ，過 \(P\) 及鄰近兩點作一個圓，當鄰點趨近 \(P\) 時圓有極限位置（密切圓）。圓越小，曲線在 \(P\) 附近彎得越急。曲率 = 密切圓半徑的倒數。直線的密切圓半徑無窮大，曲率為 \(0\) 。

平面曲線曲率公式的推導#

設平面曲線寫成函數圖像 \(y = f(x)\) ，即參數化 \(\gamma(x) = (x, f(x))\) 。弧長參數下曲率是 \(\kappa = \|dT/ds\|\) ，我們把這個定義改寫到一般參數上。

步驟一：速度與速率

\[ \gamma'(x) = (1, f'), \qquad \|\gamma'\| = \sqrt{1 + f'^2} \]

步驟二：加速度

\[ \gamma''(x) = (0, f'') \]

步驟三：代入一般公式（下一小節推導） \(\kappa = \|\gamma'\times\gamma''\| / \|\gamma'\|^3\) 。平面外積的 \(z\) 分量就是 \(1\cdot f'' - f'\cdot 0 = f''\) ，因此

\[ \boxed{\ \kappa(x) = \frac{|f''(x)|}{\left(1 + f'(x)^2\right)^{3/2}}\ } \]

兩個立即的檢查：

  • 若 \(f\) 是線性函數， \(f''=0\) ，曲率為 0。直線果然不彎。
  • \(y'^2\) 項的存在表示：同一形狀「斜著放」時 \(y'\) 變大、分母變大，但分子也會變——最後幾何不變量扭回同樣的值（讀者可在實驗驗證）。

空間曲線的一般曲率公式#

回顧 01.2 節的結論：弧長參數下 \(\kappa = \|\gamma''(s)\|\) 。設 \(\gamma = \gamma(t)\) 是一般參數， \(s(t)\) 為弧長函數。由連鎖律：

\[ \frac{d\gamma}{ds} = \frac{\gamma'}{\| \gamma'\|} = \frac{v}{\|v\|}, \qquad \frac{d^2\gamma}{ds^2} = \frac{1}{\|v\|}\frac{d}{dt}\!\left(\frac{v}{\|v\|}\right) = \frac{\|v\|^2 a - (v\cdot a)\,v}{\|v\|^4} \]

其中 \(v=\gamma'\) 、 \(a=\gamma''\) 。取範數，利用 \(\|v\|^2\|a\|^2 - (v\cdot a)^2 = \|v\times a\|^2\) （Lagrange 恆等式），得：

\[ \kappa(t) = \frac{\|\gamma'(t)\times\gamma''(t)\|}{\|\gamma'(t)\|^3} \]

這個形式的妙處： 1. 只需要一階與二階導數。 2. 完全不涉及弧長的反解。 3. 分母 \(\|\gamma'\|^3\) 自動補償了「走太快」造成的砍彎假象。

SymPy 實作：平面曲線#

圖像式曲線 \(y=f(x)\)#

import sympy as sp

x = sp.symbols('x', real=True)

def curvature_graph(f, x):
    # 以 t=x 參數化：γ(t) = (t, f(t))
    d1 = sp.diff(f, x)
    d2 = sp.diff(f, x, 2)
    kappa = sp.Abs(d2) / (1 + d1**2)**sp.Rational(3, 2)
    return sp.simplify(kappa)

# 拋物線 y = x^2
f = x**2
print(curvature_graph(f, x))      # 2/(4*x**2 + 1)**(3/2)

注意在 \(x=0\) （頂點）處 \(\kappa = 2\) ，恰好等於「以曲率半徑 \(1/2\) 的圓貼合頂點」。

參數式平面曲線#

def curvature_parametric(gx, gy, t):
    d1 = sp.Matrix([sp.diff(gx, t), sp.diff(gy, t)])
    d2 = sp.Matrix([sp.diff(gx, t, 2), sp.diff(gy, t, 2)])
    cross_z = sp.diff(gx, t)*sp.diff(gy, t, 2) - sp.diff(gy, t)*sp.diff(gx, t, 2)
    speed = sp.sqrt(d1.dot(d1))
    return sp.simplify(sp.Abs(cross_z) / speed**3)

t = sp.symbols('t', real=True)
# 圓 r = (cos t, sin t)
print(curvature_parametric(sp.cos(t), sp.sin(t), t))    # 應為 1（單位圓）

\|norm\| 的替代寫法#

SymPy 的 sp.norm(M) 等價於 sp.sqrt(M.dot(M))，對 2D/3D 向量都適用。

M = sp.Matrix([sp.cos(t), sp.sin(t), t])
print(sp.norm(M.diff(t)))   # sqrt(2)

SymPy 實作：空間曲線與 cross product#

import sympy as sp

def curvature_3d(g, t):
    v = g.diff(t)         # γ'
    a = v.diff(t)         # γ''
    cross = v.cross(a)    # γ' × γ''
    return sp.simplify(sp.sqrt(cross.dot(cross)) / sp.norm(v)**3)

# 3D 位置向量記號
X, Y, Z = sp.symbols('X Y Z', real=True)
helix = sp.Matrix([sp.cos(t), sp.sin(t), t])
print(curvature_3d(helix, t))    # 1/2

Matrix 的 .cross 方法完整支援三維外積。對 2D 曲線用 .cross 前先把向量「升維」到 3D（補 \(z=0\) ），或直接用上面的行列式寫法。

進一步：simplify 的化簡力#

有時 curvature_3d 直接輸出帶一堆根號與平方項的龐大式子。常用策略：

def curvature_3d_pretty(g, t):
    v = g.diff(t); a = v.diff(t)
    c = v.cross(a)
    expr = sp.sqrt(sp.simplify(c.dot(c))) / sp.sqrt(sp.simplify(v.dot(v)))**3
    return sp.simplify(sp.radsimp(expr))

或者先 sp.trigsimp 消去 \(\sin^2+\cos^2\) ，這對圓類與螺旋類曲線特別有效。

動手實驗#

實驗 1：三條曲線的曲率速查表

import sympy as sp
t = sp.symbols('t', real=True)

def kappa(g):
    v = g.diff(t); a = v.diff(t)
    c = v.cross(a)
    return sp.simplify(sp.sqrt(sp.trigsimp(c.dot(c))) / sp.sqrt(sp.trigsimp(v.dot(v)))**3)

parabola_sp = sp.Matrix([t, t**2, 0])                 # 空間參數化的拋物線
circle_3d   = sp.Matrix([2*sp.cos(t), 2*sp.sin(t), 0])  # 半徑 2 的圓
helix       = sp.Matrix([sp.cos(t), sp.sin(t), t])     # 單位螺旋線

for name, g in [("拋物線", parabola_sp), ("圓", circle_3d), ("螺旋線", helix)]:
    print(name, ":", kappa(g))

觀察輸出：拋物線 \(\kappa(t)\) 隨 \(|t|\) 衰減（兩翼攤平）；半徑 2 的圓 \(\kappa = 1/2\) （半徑的倒數）；螺旋線 \(\kappa = 1/2\) （呼應公式 \(\kappa = a/(a^2+b^2)\) ，此處 \((a,b)=(1,1)\) ）。

實驗 2：拋物線頂點曲率的極限檢查

x = sp.symbols('x', real=True)
f = x**2
k = sp.Abs(sp.diff(f, x, 2)) / (1 + sp.diff(f, x)**2)**sp.Rational(3, 2)
print("頂點 κ =", k.subs(x, 0))          # 2
print("頂點密切圓半徑 =", 1/k.subs(x, 0))  # 1/2

再比較：在頂點畫半徑 \(1/2\) 的圓（圓心 \((0, 1/2)\) ），與拋物線 \(\gamma(x)=(x, x^2)\) 用 matplotlib 疊圖，目視貼合程度。

實驗 3：參數化不變性的驗證

t2 = sp.symbols('t2', real=True)
# 同一條拋物線，兩種走法
g_slow = sp.Matrix([t, t**2, 0])
g_fast = sp.Matrix([2*t2, 4*t2**2, 0])
print(kappa(g_slow))   # 2/(4*t**2+1)**(3/2)
print(sp.simplify(kappa(g_fast).subs(t2, t/2)))   # 化成 t 的式子後結果相同

實驗 4：（挑戰）用曲率公式反推一個生物問題

河口河道/河流的曲率大概滿足「波形 \(\kappa(s)\propto\sin\) 」的近似。給 \(\kappa(s)=\tfrac12\sin s\) 與 \(\tau=0\) ，使用 Frenet–Serret ODE（見 01.2 節）以 sp.dsolve 或數值求解重建曲線，畫出來看看像不像河道蜿蜒。

小結與練習#

本節重點

  1. 平面圖像曲線： \(\kappa = |y''| / (1+y'^2)^{3/2}\) ，分子量「二階彎」，分母補償「斜率放大」。
  2. 一般參數空間曲線： \(\kappa = \|\gamma'\times\gamma''\| / \|\gamma'\|^3\) ，由連鎖律與 Lagrange 恆等式導出。
  3. SymPy 實作三關鍵：diff、.cross、simplify/trigsimp。
  4. 曲率是參數化不變量：換走法不變（SymPy 可輕鬆驗證）。

練習

  1. 推導：對極坐標曲線 \(r(\theta)\) ，曲率公式為 \(\kappa = \dfrac{|r^2 + 2r'^2 - r\,r''|}{(r^2+r'^2)^{3/2}}\) ，並用 SymPy 以對數螺線 \(r=e^\theta\) 驗證。
  2. 計算三次曲線 \(y = x^3\) 在原點的曲率，並用 01.2 節的想法解釋為什麼「曲率為零卻不平」（這是拐點）。
  3. 對曲線 \(\gamma(t) = (t, t^2, t^3)\) 求 \(\kappa(t)\) 與 \(B'\) 的範數（即 \(|\tau|\) ），繪 \(\kappa, \tau\) 隨 \(t\) 的關係。
  4. 證明： \(\kappa\) 在參數變換 \(t \mapsto \varphi(t)\) （ \(\varphi'>0\) ）下不變。用 SymPy 對 \(\gamma(t)=(t, t^2, t^3)\) 與 \(\varphi(t)=t^3+t\) 驗證。
  5. 實驗：畫出 \(f(x) = \sin(x)\) 一個週期的曲率 \(\kappa(x)\) 曲線，找出曲率最大的位置（提示：sp.solve(sp.diff(k, x), x) 加數值驗證），與直覺對照。