03.4 sympy.vector 實戰#

本節目標#

前三節我們手搓了梯度、散度、旋度的符號計算；現在交給官方模組 sympy.vector，並體驗它最殺手級的功能：

CoordSys3D 座標系物件與 delop（ \(\nabla\) 算子）。
一行呼叫 delop.gradient / divergence / curl 完成所有微分幾何運算。
直角／圓柱／球座標自動切換：不用再背球座標的 curl 展開公式！
保守場判驗與位能還原（potential 函式）。
實驗：驗證重力場是保守場。

3.4.1 CoordSys3D 與 delop：把 ∇ 變成物件#

import sympy as sp
from sympy.vector import CoordSys3D, delop

C = CoordSys3D("C")                     # 直角座標系，基底為 C.i, C.j, C.k
x, y, z = C.x, C.y, C.z                 # 純量：座標函數

f = x**2 + y**2 + z**2
grad_f = delop.gradient(f)              # 取得未化簡的 del 表示
print(grad_f.doit())                    # 2*C.x*C.i + 2*C.y*C.j + 2*C.z*C.k

delop 就是 \(\nabla\) 的化身，它是一個算子：delop.gradient(f)、delop.dot(grad(F))、delop.cross(F) 分別對應 \(\nabla f\) 、 \(\nabla\cdot\mathbf{F}\) 、 \(\nabla\times\mathbf{F}\) ，回傳 sympy 物件後用 .doit() 求值。

向量場的表示用基底向量相加：

F = C.x*C.i - C.y*C.j + (x*y)*C.k       # 依賴座標的向量場
print(F)                                # C.x*C.i - C.y*C.j + C.x*C.y*C.k

3.4.2 gradient / divergence / curl 一行式#

sympy.vector 提供等價的函式呼叫：

from sympy.vector import gradient, divergence, curl, laplacian

print(gradient(f))                      # 2*C.x*C.i + 2*C.y*C.j + 2*C.z*C.k
print(divergence(C.x*C.i - C.y*C.j))    # 2 - 1? 來驗證：div(x, -y, 0) = 1 + (-1) = 0
Fv = C.x*C.i - C.y*C.j
print(sp.simplify(divergence(Fv)))      # 0（不可壓縮）

G = -C.y*C.i + C.x*C.j                  # 渦旋場
print(sp.simplify(curl(G)))             # 2*C.k  （上節手算一致！）
print(sp.simplify(laplacian(f)))        # 6

兩種寫法（delop.gradient(f) 與 gradient(f)）結果相同，差別是前者顯式展示「 \(\nabla\) 作用中」的過程，適合教學展示；後者一步到位，適合計算。

3.4.3 座標系自動切換：不用背公式！#

這是本模組的精華。手算球座標下的 curl 需要背長串公式：

\[ (\nabla\times\mathbf F)_\phi = \frac{1}{r}\!\left[\frac{\partial (r F_\theta)}{\partial r} - \frac{\partial F_r}{\partial \theta}\right] \ldots \]

在 sympy.vector 裡，只要宣告球座標系並在該系表達場，curl/gradient/divergence 全部自動處理度量與尺度因子：

from sympy.vector import CoordSys3D, curl, gradient, divergence

R, T, P = sp.symbols("r theta phi", positive=True)
S = CoordSys3D("S", transformation="spherical")   # 球座標
r, th, phi = S.r, S.theta, S.phi

# 場以球座標基底表示：grad(f) in spherical
f_sph = r**2 * sp.sin(th)
print(gradient(f_sph).doit())       # 自動給出 (2 r sinθ) e_r + (r cosθ) e_θ
                                     # 公式： ∂f/∂r e_r + (1/r ∂f/∂θ) e_θ ✓

# 球座標下的場的 curl：F = e_r × ? 取 F = (r cosφ) e_ti ／使用 radial 例
Fr = r*sp.cos(phi)
print(curl(Fr*S.i).doit())           # 自動處理尺度因子，與手算拱積公式一致

3.4.4 保守場與位能還原#

保守場的定義鏈： \(\nabla\times\mathbf{F} = 0\) （單連通區域）⇒ \(\mathbf{F} = \nabla \phi\) ⇒ 線積分與路徑無關，只看端點。sympy.vector 提供 potential 直接把 \(\phi\) 找回來：

from sympy.vector import potential, is_conservative, is_solenoidal

Fcons = 2*x*y*C.i + (x**2 + z)*C.j + y*C.k
print(is_conservative(Fcons))     # True（curl 為零）
print(is_solenoidal(Fcons))       # False（div = 2y ≠ 0）

phi = potential(Fcons)            # 還原位能
print(phi)                        # x**2*y + y*z
print(gradient(phi - phi).doit()) # 0 —— 檢驗 grad(phi) = F

若 is_conservative 回傳 False，代表場非保守，位能不存在。函式 is_solenoidal 則檢查無源。

3.4.5 動手實驗：驗證重力場是保守場#

理想（均勻）重力場 \(\mathbf{F} = -mg\,\mathbf{k}\) ，以及更真實的牛頓重力場 \(\mathbf{F} = -\dfrac{GMm}{r^3}\mathbf{r}\) ，都應是保守的：

m, g = sp.symbols("m g", positive=True)
Fg = -m*g*C.k
print(is_conservative(Fg))                 # True
print(potential(Fg))                       # m*g*z + C (位能 = mgh ✓)

# 牛頓重力（球座標優雅計算）
Gm = sp.symbols("G_m", positive=True)
N = CoordSys3D("N", transformation="spherical")
r_s = N.r
Fnewton = -Gm/(r_s**2) * N.i               # 徑向場 -G/r^2 e_r
print(is_conservative(Fnewton))            # True
print(sp.simplify(potential(Fnewton)))     # G_m/r + C  ⇒ 位能 -G/r ✓

物理意義：不同路徑把石頭抬高到同一點，重力做的功相同——這是所有「依高度定位能而不是依路徑」的物理直覺基礎。用第二類線積分交叉驗證（3.3 節方法）：

t = sp.symbols("t", real=True)
path = sp.Matrix([t, t**2, 1 - t])          # 任意蜿蜒路徑
Fv = sp.Matrix([0, 0, -m*g])
W = sp.integrate(Fv.dot(sp.diff(path, t)), (t, 0, 2))
print(W)                                    # 2*m*g = m*g*(z_start - z_end) ✓ 路徑無關

3.4.6 座標變換速查#

S = CoordSys3D("S", transformation="spherical")
print(S.position().doit())    # 位置向量在球座標下的笛卡爾形式
print(S.base_vectors())       # (e_r, e_theta, e_phi) 三個基底
print(S.jacobian_vectors())   # SHAPE 因子相關的資訊

小結與練習#

CoordSys3D + delop 把 \(\nabla\) 變成可操作物件；gradient/divergence/curl 一行完成。
對 transformation="spherical" / "cylindrical" 校準後，所有運算自動套用正確公式——這是最大賣點。
is_conservative + potential 完成「驗旋度為零 ⇒ 找位能」的整條鏈。

練習：

用 delop 與函式版兩種方法計算 \(f = xyz\) 的梯度，結果一致嗎？
對 \(\mathbf{F} = (yz, xz, xy)\) ：is_conservative 驗證、用 potential 求位能，再手算確認。
在圓柱座標系（transformation="cylindrical"）下計算 \(\nabla \cdot (r\,\mathbf{e}_r)\) 與 \(\nabla\times(r\,\mathbf{e}_r)\) 。
用 potential 找出一個保守場的位能，畫出其三維曲面（matplotlib plot_surface），說明位能谷對應向下方向。
挑戰：構造 \(\mathbf{F} = \nabla(\frac{1}{r})\) ，is_conservative 應回 True，但通量不為零（奇點）——用 divergence 與手工通量積分對照，體會「局部 vs 全域」的差異。