9.2 拉普拉斯方程與勢論（1843 年）#
學習動機#
9.1 節的兩個方程——熱方程與波方程——都有「方向」：熱朝著低溫流，波朝著某個方向傳。它們是演化方程：給了初值，就一路演化下去。
但物理世界裡還有一類問題沒有時間：
- 天體的引力場在哪裡？給定一堆質量的分佈，求每一點的引力加速度。
- 靜電場在哪裡？給定一堆電荷，求每一點的電場。
- 流體以不可壓縮方式流過什麼形狀的管道？求速度場。
這三個問題的共同點是：沒有初值，只有邊界（或分佈），而且它們的方程都是橢圓型的。第二個問題的答案，是拉普拉斯方程：
它是全書最重要的一條方程，因為它是唯一一個「無源無匯」的場」的方程。任何沒有散度的向量場，都可以用一個標量位能 \(u\) 表示成 \(\mathbf{F} = -\nabla u\) ，而 \(\nabla \cdot \mathbf{F} = 0\) 就變成了 \(\nabla^2 u = 0\) 。
歷史動機：1843 年，皇家學會（Royal Society）重印了 George Green 1828 年那篇論文的英譯本——〈論把數學分析應用於吸引與電的理論〉。這篇被冷落了整整十五年的短文，裡面包含了 Green 恆等式、Green 函數、單極定理、勢函數的基本思想，以及「用邊界值決定內部場」的整套方法。勢論（potential theory）從此成為一門獨立的學科，並在二十年内長成數學、物理、天文的共同語言。
本章的 SymPy 目標是：
- 從 Poisson 方程出發，逐字驗證均質球、球殼、點源的勢函數；
- 用 Green 恆等式符號地證明 Dirichlet 問題解的唯一性（能量法）；
- 用最大原理說明「為什麼橢圓型方程只認邊界」；
- 用符號運算檢查單位圓盤的 Green 函數與 Poisson 核的每一個性質；
- 用 SOR 迭代實際解一個 Poisson 邊界值問題，並觀察二階收斂。
9.2.1 Poisson 的方程（1843）：把萬有引力變成偏微分方程#
從「力」到「位能」#
牛頓的萬有引力定律給的是力： \(\mathbf{F}(\mathbf{r}) = -Gm\,\dfrac{\mathbf{r}}{|\mathbf{r}|^3}\) 。這是一個點對點的關係：要知道 \(\mathbf{r}\) 處的力，得知道所有源的位置，然後逐一加總。對兩個物體很方便，對一千個就很笨拙。
位能（potential）的想法來自一個觀察：位能可以「疊加」，而且由位能取梯度就得到力。定義比勢 \(\Phi(\mathbf{r})\) 滿足 \(\mathbf{g}(\mathbf{r}) = -\nabla \Phi(\mathbf{r})\) ，其中 \(\mathbf{g}\) 是重力加速度。對點質量 \(m\) 而言， \(\Phi(\mathbf{r}) = -\dfrac{Gm}{|\mathbf{r}|}\) （這裡 \(\Phi\) 是「單位質量」的勢能，所以 \(\nabla\Phi\) 直接就是加速度。）
現在關鍵的一步：把「由源疊加」換成「解一個方程」。
Poisson 方程的推導#
設質量分佈為 \(\rho(\mathbf{r})\) ，比勢為 \(\Phi\) 。我們要找的不是 \(\Phi\) 的值，而是它的散度：
由牛頓定律，點源的散度是「負的 \(4\pi G\) 乘上 delta 函數」，用分部積分可以推廣到連續分佈：
這就是 Poisson 方程。（慣例差異：有些書定義 \(\Phi\) 為「每單位質量的勢」並寫 \(\nabla^2\Phi = 4\pi G\rho\) ；有些定義勢能 \(U = m\Phi\) 寫 \(\nabla^2 U = -4\pi Gm\rho\) ；符號相反，物理內容相同。）
Poisson 方程的威力在於：它把「對每個位置求和」變成「求一個函數」。而且它是局部的： \(\Phi\) 在某點的值，只由那附近的 \(\rho\) 直接決定（更精確地說，由一個局部加權平均決定）。
用 SymPy 驗證：均質球的勢#
我們完整算一遍。物理結論是教科書上的標準結果：
- 球外（ \(r > R\) ）： \(\Phi = -GM/r\) ，等同於把全部質量集中在球心；
- 球內（ \(r < R\) ）： \(\Phi = -\dfrac{G}{2}\left(3R^2 - \dfrac{r^2}{3}\right)\rho\) 。
import sympy as sp
r, R, G, rho = sp.symbols('r R G rho', positive=True)
eps0 = sp.Symbol('eps_0', positive=True)
M = sp.Rational(4, 3)*sp.pi*rho*R**3
# ---------- (1) Poisson 方程：均質球的比勢 ----------
# 慣例 Phi = 單位質量勢能，g = -dPhi/dr，Phi(∞) = 0
Phi_out = -G*M/r
Phi_in = -G*(2*sp.pi*rho)*(R**2 - r**2/3)
print('球外 Phi =', sp.simplify(Phi_out))
print('球內 Phi =', sp.simplify(Phi_in))
lap = lambda f: sp.simplify(sp.diff(f, r, 2) + 2/r*sp.diff(f, r))
print('球內 ∇²Phi =', lap(Phi_in), ' Poisson: ∇²Phi = 4πG rho ?',
sp.simplify(lap(Phi_in) - 4*sp.pi*G*rho) == 0)
print('球外 ∇²Phi =', lap(Phi_out))
# ---------- (2) 邊界上的連續性 ----------
print('Phi 在 r=R 連續？ 差 =', sp.simplify((Phi_in - Phi_out).subs(r, R)))
print("Phi' 在 r=R 連續？ 差 =",
sp.simplify((sp.diff(Phi_in, r) - sp.diff(Phi_out, r)).subs(r, R)))
# ---------- (3) 重力加速度 ----------
g_out = sp.simplify(-sp.diff(Phi_out, r))
g_in = sp.simplify(-sp.diff(Phi_in, r))
print('球外 g =', g_out, ' 等於 -GM/r² ?', sp.simplify(g_out + G*M/r**2) == 0)
print('球內 g =', g_in, ' 等於 -(4πG/3) rho r ?',
sp.simplify(g_in + sp.Rational(4, 3)*sp.pi*G*rho*r) == 0)
# ---------- (4) 靜電學的對應：點電荷的電位與電場 ----------
q = sp.Symbol('q', positive=True)
phi = q/(4*sp.pi*eps0*r)
E = sp.simplify(-sp.diff(phi, r))
print()
print('點電荷的電位 phi =', phi)
print('電場 E = -dphi/dr =', E, ' 等於 q/(4π ε0 r²) ?',
sp.simplify(E - q/(4*sp.pi*eps0*r**2)) == 0)
print('球面通量 ∮E·n dS =', sp.simplify(4*sp.pi*r**2*E), ' （Gauss 定律：等於 q/eps0）')
# ---------- (5) 單極定理：單位點源的通量恆為 1 ----------
P = -1/(4*sp.pi*r)
print()
print('單位點源 Phi =', P, ' r>0 處 ∇²Phi =', lap(P))
print('球面通量 ∮∂Phi/∂n dS =', sp.simplify(4*sp.pi*r**2*sp.diff(P, r)),
' （恆為 1，與 r 無關）')
預期輸出：
球外 Phi = -4*pi*G*R**3*rho/(3*r) 球內 Phi = 2*pi*G*rho*(-3*R**2 + r**2)/3 球內 ∇²Phi = 4*pi*G*rho Poisson: ∇²Phi = 4πG rho ? True 球外 ∇²Phi = 0 Phi 在 r=R 連續？ 差 = 0 Phi' 在 r=R 連續？ 差 = 0 球外 g = -4*pi*G*R**3*rho/(3*r**2) 等於 -GM/r² ? True 球內 g = -4*pi*G*r*rho/3 等於 -(4πG/3) rho r ? True 點電荷的電位 phi = q/(4*pi*eps_0*r) 電場 E = -dphi/dr = q/(4*pi*eps_0*r**2) 等於 q/(4π ε0 r²) ? True 球面通量 ∮E·n dS = q/eps_0 （Gauss 定律：等於 q/eps0） 單位點源 Phi = -1/(4*pi*r) r>0 處 ∇²Phi = 0 球面通量 ∮∂Phi/∂n dS = 1 （恆為 1，與 r 無關）
逐行解讀：
- 第一、二行：符號輸出的形式 \(\dfrac{-4\pi G R^3 \rho}{3r}\) 正是 \(-GM/r\) （因為 \(M = \frac{4}{3}\pi\rho R^3\) ）。我們不是「宣稱」球外勢等於點源的勢，而是把它寫成一個符號表達式，再用
sp.simplify化簡出可核對的形式。 - 第三行「球內 \(\nabla^2\Phi = 4\pi G\rho\) 」：這是 Poisson 方程的直接驗證。注意 \(\nabla^2\) 在球座標下對球對稱函數是 \(\frac{\mathrm{d}^2}{\mathrm{d}r^2} + \frac{2}{r}\frac{\mathrm{d}}{\mathrm{d}r}\) —— 這正是第 7 章座標變換的成果。「Poisson 方程在球內自動成立」這件事，用 SymPy 三行就驗完了。
- 第四、五行是本章的關鍵物理： \(\Phi\) 與 \(\Phi'\) 在 \(r = R\) 都連續。這意味著：球面上沒有「勢的跳躍」，也就沒有「表面上的力」。對照後面的電場版本——電場在帶電面上確實會跳躍（這就是高斯定律的物理內容）。勢連續、場可跳躍，這個區分在 9.2.3 節會派上用場。
- 第六、七行：球內重力加速度 \(\mathbf{g} = -\frac{4\pi G\rho}{3}r\,\hat{\mathbf{r}}\) 是線性的，這就是牛頓「把地球當成一個巨大的球」能算出 \(\mathbf{g} \approx 9.8\ \mathrm{m/s^2}\) 的原因。
- 第八到十行（靜電學）：把 \(G\rho\) 換成 \(q/\varepsilon_0\) ，一切照樣成立。這不是巧合：靜電學就是牛頓引力學換掉一個常數。這也解釋了為什麼 19 世紀的人相信「電和引力是同一種東西的不同表現」——直到實驗否定了這個想法。
- 最後兩行是單極定理：單位點源的球面通量恆為 1，與球半徑無關。這是 Gauss 定律的原型，也是 9.2.2 節的主角。
9.2.2 Green 公式與單極定理：1828 年的那篇短文#
Green 1828：被忽略十五年的短文#
1828 年 12 月，George Green（1793–1841）在《皇家學會哲學會刊》發表了一篇只有 31 頁的小論文〈論把數學分析應用於吸引與電的理論〉。論文的目的是證明地球周圍的引力位 behaves like a point mass（牛頓的「地球是球殼」猜想），方法則是把位能表示成面積分——也就是今天的 Green 函數表示。論文在生前幾乎無人問津；1843 年皇家學會重印了它的英譯本，此後 William Thomson（Kelvin 勳爵）把它與 Faraday 的實驗工作聯繫起來，勢論才真正「活」了過來。
Green 恆等式#
Green 的核心工具是分部積分的一個巧妙組合。對兩個足夠光滑的函數 \(u, v\) 在區域 \(D\) 上：
這就是 Green 第二恆等式。它的三個用法，都將在後面出現：導出 Green 函數（把 \(v\) 換成「奇點解」）；證明 Dirichlet 問題解的唯一性（把兩個解相減）；導出 Green 公式（解 Poisson 方程的積分公式）。
我們用 SymPy 在一個矩形上把它從頭算一遍：
import sympy as sp
x, y = sp.symbols('x y', real=True)
lap = lambda f: sp.diff(f, x, 2) + sp.diff(f, y, 2)
# ---------- (1) Green 第二恆等式：矩形 D = [0,1]^2 ----------
u = x*y
v = x**2 + y**2
lhs = sp.integrate((u*lap(v) - v*lap(u)), (y, 0, 1), (x, 0, 1))
print('左邊 ∫∫(uΔv - vΔu) dA =', sp.simplify(lhs))
def edge_int(u, v, coord, value, nx, ny, direction):
if coord == 'x':
f = u*sp.diff(v, x) - v*sp.diff(u, x)
return sp.integrate((nx*f).subs(x, value), (direction, 0, 1))
f = u*sp.diff(v, y) - v*sp.diff(u, y)
return sp.integrate((ny*f).subs(y, value), (direction, 0, 1))
rhs = (edge_int(u, v, 'x', 1, 1, 0, y) + edge_int(u, v, 'x', 0, -1, 0, y)
+ edge_int(u, v, 'y', 1, 0, 1, x) + edge_int(u, v, 'y', 0, 0, -1, x))
print('右邊 ∮(u∂v/∂n - v∂u/∂n) ds =', sp.simplify(rhs))
print('左右相等？', sp.simplify(lhs - rhs) == 0)
# ---------- (2) 唯一性定理（能量法）----------
print()
print('Green 第一恆等式：∫∫|∇w|² dA = ∮ w ∂w/∂n ds - ∫∫ w Δw dA')
print('若 w|邊界 = 0 且 Δw = 0，則 ∫∫|∇w|² dA = 0，故 w ≡ 0')
w0 = x*(1 - x)*y*(1 - y)
print('反例 w0 = x(1-x)y(1-y)：邊界為 0，但 ∇²w0 =', sp.simplify(lap(w0)))
# ---------- (3) 最大原理與五點格式 ----------
avg = ((x + 1)**2 - y**2 + (x - 1)**2 - y**2 + x**2 - (y + 1)**2 + x**2 - (y - 1)**2)/4
print()
print('五點平均值算子作用於 u = x²-y² 的誤差 =', sp.simplify(avg - (x**2 - y**2)))
print('（差分算子只取二階項，所以對二次多項式恰為零）')
預期輸出：
左邊 ∫∫(uΔv - vΔu) dA = 1 右邊 ∮(u∂v/∂n - v∂u/∂n) ds = 1 左右相等？ True Green 第一恆等式：∫∫|∇w|² dA = ∮ w ∂w/∂n ds - ∫∫ w Δw dA 若 w|邊界 = 0 且 Δw = 0，則 ∫∫|∇w|² dA = 0，故 w ≡ 0 反例 w0 = x(1-x)y(1-y)：邊界為 0，但 ∇²w0 = 2*x*(x - 1) + 2*y*(y - 1) 五點平均值算子作用於 u = x²-y² 的誤差 = 0 （差分算子只取二階項，所以對二次多項式恰為零）
注意 edge_int 這個小函數：它把「對邊界積分」拆成四條邊，每條邊手動指定外法線 \((\pm 1, 0)\) 或 \((0, \pm 1)\) —— 這是 8.2 節格林定理的應用。符號運算在這裡的價值是：外法線的方向符號不會出錯，而這正是手算最容易搞錯的地方。
第一恆等式與唯一性：令 \(v = u\) ，Green 第二恆等式變成
移項即得 Green 第一恆等式（上半部形狀的那個）。這就是能量法：
定理（Dirichlet 問題解的唯一性） 若 \(u_1, u_2\) 都是 \(\nabla^2 u = 0\) 在 \(D\) 中的解，且在 \(\partial D\) 上取值相同，則 \(u_1 \equiv u_2\) 。
證明（用 SymPy 的輸出就能讀懂）：令 \(w = u_1 - u_2\) ，則 \(\nabla^2 w = 0\) 且 \(w |_{\partial D} = 0\) 。用第一恆等式：
所以 \(\nabla w \equiv 0\) ，即 \(w\) 是常數；又 \(w\) 在邊界上為零，所以 \(w \equiv 0\) 。∎
「反例 \(w_0\) 」那一行的意義：它提醒我們，唯一性需要兩個條件同時成立——邊界值相同且滿足方程。 \(w_0 = x(1-x)y(1-y)\) 在邊界上確實為零，但它不滿足 \(\nabla^2 w_0 = 0\) （輸出顯示 \(\nabla^2 w_0 = 2x(x-1) + 2y(y-1)\) ，在內部不為零），所以它不是反例。
單極定理：Gauss 定律的原型#
單極定理（Green 1828）：若 \(u\) 在區域 \(D\) 內除了一點 \(P\) 之外處處調和，則
特別地，單位點源 \(u = -\dfrac{1}{2\pi}\log r\) 的通量恆為 \(-1\) 。
我們在 9.2.1 的程式碼裡已經看到了這個結果的最簡版本： \(-\frac{1}{4\pi r}\) 的球面通量恆為 \(1\) 。這個「通量與 r 無關」的性質，是「點源在外部看起來像什麼」的完整定量答案——它對應物理上的兩件事：
- 萬有引力：球外重力 \(-GM/r^2\) （平方反比）；
- 庫侖定律：球外電場 \(q/(4\pi\varepsilon_0 r^2)\) （平方反比）。
平方反比定律不是猜的，是從 Green 恆等式推出來的。 這是一個深刻的歷史事實：牛頓在《原理》（1687）裡用平方反比是經驗歸納；Green 在 1828 年證明了只有平方反比與拉普拉斯方程相容。
9.2.3 勢論如何連結電磁學：Maxwell 的橋樑#
從「力」到「場」的轉向#
19 世紀的物理學有一場悄悄但深刻的轉變：不再問「這個力和那個力相等嗎」，而是問「場在各處是多少」。法拉第（Faraday）在 1830 年代的《皇家學院講義》中提出了「力線」的概念；Maxwell 在 1865 年的〈電磁場的動力學理論〉中把它數學化。
Maxwell 的核心洞見是：電場與磁場可以由兩個「位能」產生——標量電位 \(\phi\) 與向量位 \(\mathbf{A}\) ：
這裡的 \(\phi\) 就是靜電學的位能，而它滿足的方程正是我們剛學的 Poisson 方程（在無自由電荷處是 Laplace 方程）： \(\nabla^2 \phi = -\dfrac{\rho}{\varepsilon_0}\) 。Poisson 方程 + 格林定理 + 邊界條件 = 電磁場的完整理論。
兩個勢論的對照表#
| 引力（Newton 1687） | 靜電（Coulomb 1785） | |
|---|---|---|
| 源 | 質量密度 \(\rho\) | 電荷密度 \(\rho_e\) |
| 方程 | \(\nabla^2 \Phi = 4\pi G\rho\) | \(\nabla^2 \phi = -\rho_e/\varepsilon_0\) |
| 場 | \(\mathbf{g} = -\nabla \Phi\) | \(\mathbf{E} = -\nabla \phi\) |
| 單位點源 | \(-G/r\) | \(q/(4\pi\varepsilon_0 r)\) |
| 球面通量 | \(-4\pi GM\) | \(q/\varepsilon_0\) |
| 邊界條件 | 連續（無表面密度） | 可跳躍（有面電荷） |
為什麼「勢連續、場可跳躍」？#
這一點值得單獨說清楚，因為它是 9.2.1 節程式碼輸出的直接推論。在物體的表面上，如果有面密度 \(\sigma\) （例如一張帶電的薄殼），則：
- 電位 \(\phi\) 在跨越表面時連續（因為 \(\phi\) 是 \(\frac{1}{4\pi\varepsilon_0}\int \frac{\sigma}{|\mathbf{r}-\mathbf{r}'|}\,\mathrm{d}A'\) ，積分沒有奇異性）；
- 電場 \(\mathbf{E}\) 在跨越表面時跳躍 \(\dfrac{\sigma}{\varepsilon_0}\) （因為 \(\mathbf{E}\) 的核 \(\frac{1}{4\pi\varepsilon_0}\frac{\mathbf{r}-\mathbf{r}'}{|\mathbf{r}-\mathbf{r}'|^3}\) 在源點有奇異性）。
這正是高斯定律： \(\oint \mathbf{E}\cdot\mathbf{n}\,\mathrm{d}S = q_{\text{enc}}/\varepsilon_0\) —— 跨越一個面電荷時法向通量的跳躍。
物理直覺：電位像溫度，電場像熱流——溫度在界面連續，熱流則會跳變。對照引力的情形：地球的質量分佈只有體密度、沒有面密度，所以 \(\Phi\) 與 \(\mathbf{g}\) 都連續，沒有跳躍。
9.2.4 邊界值問題與唯一性：橢圓型方程的兩大定理#
Dirichlet 問題#
Dirichlet 問題：給定區域 \(D\) 的邊界 \(\partial D\) 上的連續函數 \(f\) ，求 \(u\) 使得
這是橢圓型方程的基本問題。與 ODE 的 Dirichlet 問題（ \(u'' = 0\) 在 \([a,b]\) 上， \(u(a) = \alpha, u(b) = \beta\) ）平行，但「邊界條件」從兩個數字變成一個函數。
三個基本問題：
| 問題 | 給定 | 適用方程類型 | 解的性質 |
|---|---|---|---|
| Dirichlet | \(u\) 在邊界上的值 | 橢圓型 | 通常唯一存在 |
| Neumann | \(\partial u/\partial n\) 在邊界上的值 | 橢圓型 | 只差一個常數（否則無解） |
| Cauchy | \(u\) 與 \(\partial u/\partial n\) 都在邊界上 | 雙曲型 | 局部唯一（Illposed） |
為什麼橢圓型只有 Dirichlet 與 Neumann？ 因為橢圓型方程沒有實特徵線（9.1.3 節：特徵斜率是虛數），初值資料無法傳播。物理上，Laplace 方程描述的是「平衡態」，平衡態只由邊界條件決定。
定理一：唯一性（能量法）#
已在 9.2.2 節證明。讓我們把它寫成完整形式：
定理 若 \(u_1, u_2 \in C^2(D) \cap C(\overline{D})\) 都滿足 \(\nabla^2 u = 0\) 且在 \(\partial D\) 上兩者相等，則 \(u_1 \equiv u_2\) 。
能量法的本質：把偏導數的平方積積起來，得到一個非負量，然後證明它必須為零。這是 19 世紀「能量方法」的典型風格，在後來的存在性證明中扮演核心角色。
定理二：最大原理#
定理（強最大原理） 若 \(u\) 在 \(D\) 內調和且在 \(\overline{D}\) 上連續，則 \(u\) 的最大值與最小值都在 \(\partial D\) 上取得。若 \(u\) 在內部某點達到最大值，則 \(u\) 是常數。
這個定理是極大原理，它是 Laplace 方程最重要的定性性質。用一個「剝削」（barrier） argument 來證明：設 \(\max_{\overline D} u = M\) 在內部某點 \(P\) 達成，取以 \(P\) 為中心的球 \(B_r \subset D\) ；若 \(u\) 不是常數，則在 \(\partial B_r\) 上 \(u < M\) ；構造 \(v = u - \varepsilon + c(\text{以 } P \text{ 為心的球面函數})\) ，讓 \(v\) 在 \(\partial B_r\) 上小於等於 \(u\) 但 \(v(P) > u(P)\) ；則 \(v\) 也在內部達到最大值，與 \(v\) 調和矛盾。
最大原理的物理意義：流體的穩態溫度不會超過邊界溫度的最大值。這就是為什麼「熱方程的平衡態」與「Laplace 方程的 Dirichlet 解」是同一件事——而最大原理保證了這個平衡態是唯一的。
實際求解一個 Poisson 邊界值問題#
我們用逐次超松弛法（SOR）真的解一個 Poisson 問題，並觀察二階收斂與最大原理。
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
# ---------- 單位正方形的 Poisson 問題：Δu = -2π² sin(πx) sin(πy)，u|邊界 = 0 ----------
x, y = sp.symbols('x y', real=True)
u_exact = sp.sin(sp.pi*x)*sp.sin(sp.pi*y)
print('解析解 u = sin(πx) sin(πy)：')
print(' Δu =', sp.simplify(sp.diff(u_exact, x, 2) + sp.diff(u_exact, y, 2)))
print(' 邊界 u(0,y) =', sp.simplify(u_exact.subs(x, 0)),
' u(1,y) =', sp.simplify(u_exact.subs(x, 1)))
def sor(N, src, omega=1.9, iters=2000):
h = 1.0/(N - 1)
xg = np.linspace(0, 1, N)
X, Y = np.meshgrid(xg, xg, indexing='ij')
f = src(X, Y)
u = np.zeros((N, N))
for it in range(iters):
for i in range(1, N - 1):
for j in range(1, N - 1):
s = (u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1] - f[i, j]*h*h)/4.0
u[i, j] = (1 - omega)*u[i, j] + omega*s
return u, h, X, Y
src = lambda X, Y: -2*np.pi**2*np.sin(np.pi*X)*np.sin(np.pi*Y)
print()
print('網格收斂（SOR, omega = 1.9）：')
for N in (21, 41, 81):
u, h, X, Y = sor(N, src)
Uex = np.array([[float(u_exact.subs({x: X[i, j], y: Y[i, j]})) for j in range(N)] for i in range(N)])
print(' N = %3d h = %.4f max error = %.4e' % (N, h, np.max(np.abs(u - Uex))))
u, h, X, Y = sor(81, src)
Uex = np.array([[float(u_exact.subs({x: X[i, j], y: Y[i, j]})) for j in range(81)] for i in range(81)])
interior = u[1:-1, 1:-1]
print()
print('最大原理檢查：')
print(' 內部 最大值 = %.6f，最小值 = %.6f' % (interior.max(), interior.min()))
print(' 解析 最大值 = %.6f' % Uex.max())
print(' 邊界值全為 0，內部解嚴格為正 ?', interior.min() > 0)
print(' 內部解不超過解析最大值（誤差內）?', interior.max() <= Uex.max() + 1e-3)
fig, ax = plt.subplots(figsize=(5.4, 4.6))
im = ax.pcolormesh(X, Y, u, cmap='viridis', shading='auto')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Poisson problem on the unit square')
fig.colorbar(im, ax=ax)
fig.tight_layout()
plt.show()
預期輸出：
解析解 u = sin(πx) sin(πy)： Δu = -2*pi**2*sin(pi*x)*sin(pi*y) 邊界 u(0,y) = 0 u(1,y) = 0 網格收斂（SOR, omega = 1.9）： N = 21 h = 0.0500 max error = 2.0587e-03 N = 41 h = 0.0250 max error = 5.1420e-04 N = 81 h = 0.0125 max error = 1.2852e-04 最大原理檢查： 內部 最大值 = 1.000129，最小值 = 0.001542 解析 最大值 = 1.000000 邊界值全為 0，內部解嚴格為正 ? True 內部解不超過解析最大值（誤差內）? True（圖形由
matplotlib顯示於螢幕；若在無 GUI 環境執行，請改用plt.savefig('out.png')。）
三個觀察：
- 收斂階是 \(O(h^2)\) 。誤差從 \(2.06\times10^{-3}\) 降到 \(5.14\times10^{-4}\) 再到 \(1.29\times10^{-4}\) ，每次網格加密（ \(h\) 減半）誤差大約降為原來的 \(1/4\) —— 這正是二階收斂的特徵。對照第 6 章「一階收斂」的判別，是同一套語言。
- 最大原理在數值上被驗證：邊界值全為零，但內部解嚴格為正（最小值 \(0.0015 > 0\) ），且不超過解析最大值。符號方法給出的是「為什麼」，數值方法給出的是「確實如此」。
- **超松弛參數 \(\omega = 1.9\) ** ：Gauss–Seidel 是 \(\omega = 1\) ；超松弛讓迭代「走得有點過頭」，通常能加速收斂。 \(0 < \omega < 2\) 是穩定範圍。
順帶一提：Neumann 問題的相容條件#
Neumann 問題給定 \(\partial u/\partial n = g\) 。此時對 Green 第一恆等式取 \(v = 1\) ：
所以必須有 \(\oint_{\partial D} g\,\mathrm{d}s = 0\) （相容條件）。而且若有一個解 \(u\) ，則 \(u + C\) 也是解（ \(C\) 任意）——Neumann 問題的解只差一個常數。這是因為 \(\nabla^2 C = 0\) 而 \(\partial C/\partial n = 0\) 。物理上這就是「把整個系統的電位抬高 \(\Delta V\) ，場完全不變」。
9.2.5 Green 函數與 Perron 定理：把邊界條件變成一個積分#
Green 函數：把「解方程」變成「做積分」#
Green 函數的想法：既然 Poisson 方程的解在區域內部只「看見」邊界值，那麼我們應該能寫出
其中 \(G\) 滿足：在 \(\mathbf{r} \neq \mathbf{s}\) 處 \(\nabla^2_{\mathbf{r}} G = 0\) ；在 \(\mathbf{r} = \mathbf{s}\) 處有單位奇點 \(G \sim -\dfrac{1}{2\pi}\log|\mathbf{r}-\mathbf{s}|\) ；在邊界上 \(G = 0\) 。條件 2 的 \(-1/(2\pi)\) 是正規化的選擇，使得第 1 條的「源強度」為一。
Green 函數把「求一個函數」變成了「對邊界資料做一次加權積分」——這是 19 世紀數學最重大的「積分化」成果之一。
單位圓盤的 Green 函數#
我們不只引用公式，而是用 SymPy 從零驗證它的每一個性質。
單位圓盤 \(D = \{z : |z| < 1\}\) 的 Dirichlet Green 函數是
這個公式來自鏡像法（把 \(w\) 關於單位圓反演成 \(1/\bar{w}\) ，然後用點源的 Green 函數 \(G = -\frac{1}{2\pi}\log|z-w|\) 相減）。
import sympy as sp
x, y, a, b = sp.symbols('x y a b', real=True)
# ---------- 單位圓盤的 Dirichlet Green 函數 ----------
# G(z, w) = (1/2π) log|(1 - z w̄)/(z - w)|, w = a+ib, z = x+iy
A = (a*y - b*x)**2 + (a*x + b*y - 1)**2 # |1 - z w̄|^2
B = (x - a)**2 + (y - b)**2 # |z - w|^2
lap = lambda f: sp.simplify(sp.diff(f, x, 2) + sp.diff(f, y, 2))
print('|1 - z w̄|^2 =', A)
print('|z - w|^2 =', B)
print('A - B =', sp.factor(sp.expand(A - B)))
print(' -> 在 |z| = 1 上 A = B，故 G = 0（邊界條件）')
G = sp.log(A)/(4*sp.pi) - sp.log(B)/(4*sp.pi) # 等價於 (1/2π) log(A/B)
print('∇²G（z ≠ w 處）= 0 ?', lap(G) == 0)
# 對稱性 G(z,w) = G(w,z)
A2 = (x*b - y*a)**2 + (x*a + y*b - 1)**2
B2 = (a - x)**2 + (b - y)**2
print('G(z,w) - G(w,z) =', sp.simplify(sp.log(A)/(4*sp.pi) - sp.log(B)/(4*sp.pi)
- sp.log(A2)/(4*sp.pi) + sp.log(B2)/(4*sp.pi)))
# ---------- Poisson 核（Green 函數的邊界極限） ----------
r, th = sp.symbols('r theta', real=True)
P = (1 - r**2)/(1 - 2*r*sp.cos(th) + r**2)
print()
print('Poisson 核 P_r(θ) =', P)
print('P_r(0) =', sp.simplify(P.subs(th, 0)), ' （等於 (1+r)/(1-r)）')
for rv in (sp.Rational(1, 2), sp.Rational(1, 3)):
Pv = P.subs(r, rv)
print(' r = %s : ∫P_r dθ = %s' % (rv, sp.simplify(sp.integrate(Pv, (th, -sp.pi, sp.pi)))))
預期輸出：
|1 - z w̄|^2 = (a*y - b*x)**2 + (a*x + b*y - 1)**2 |z - w|^2 = (-a + x)**2 + (-b + y)**2 A - B = (a**2 + b**2 - 1)*(x**2 + y**2 - 1) -> 在 |z| = 1 上 A = B，故 G = 0（邊界條件） ∇²G（z ≠ w 處）= 0 ? True G(z,w) - G(w,z) = 0 Poisson 核 P_r(θ) = (1 - r**2)/(r**2 - 2*r*cos(theta) + 1) P_r(0) = (-r - 1)/(r - 1) （等於 (1+r)/(1-r)） r = 1/2 : ∫P_r dθ = 2*pi r = 1/3 : ∫P_r dθ = 2*pi
逐行解讀（這是本章最漂亮的一次符號驗證）：
A - B = (a² + b² - 1)(x² + y² - 1)：這個因式分解一箭雙雕地證明了兩件事：- 在邊界 \(|z| = 1\) 上（ \(x^2+y^2-1 = 0\) ）有 \(A = B\) ，所以 \(G = \frac{1}{4\pi}\log\frac{A}{B} = 0\) —— 邊界條件；
- 在 \(|w| = 1\) 上（ \(a^2+b^2-1 = 0\) ）同樣有 \(A = B\) —— 這說明公式對 \(w\) 也在「圓盤內部」是對稱的（配合下一行的對稱性檢查）。
∇²G = 0 ? True：這是符號地確認 Green 函數在 \(z \neq w\) 處調和。注意lap(G)裡的sp.simplify是關鍵——未化簡的 \(\frac{\partial^2}{\partial x^2}\log\big((a y - b x)^2 + (a x + b y - 1)^2\big)\) 是極其可憎的式子，sp.simplify把它化成了零。這就是符號代數的價值：不是我們算得比人快，而是機器不會在化簡時厭倦。G(z,w) - G(w,z) = 0：Green 函數的自對稱性。它在物理上意味著「源與觀測點的角色可以互換」，在數學上則是 Green 恆等式的推論。- Poisson 核： \(P_r(\theta) = \dfrac{1-r^2}{1 - 2r\cos\theta + r^2}\) 正是 Green 函數對法線求導的邊界極限。它有兩個關鍵性質：
- \(P_r > 0\) （正值性）；
- \(\int_{-\pi}^{\pi} P_r(\theta)\,\mathrm{d}\theta = 2\pi\) （守恆性，輸出以 \(r = 1/2, 1/3\) 兩個有理數值確認）。
守恆性的物理意義：把邊界上每點的「權重」加總，總共是一份完整的「份量」；而 \(P_r(0) = \frac{1+r}{1-r}\) 隨 \(r\to 1\) 爆炸——越靠近邊界，正對著你的那一小段邊界權重越大。這就是「邊界值對靠近邊界的內部點影響最大」這件事的定量表達。
Perron 定理：不用 Green 函數也能證明存在#
Perron（1923）給出了一個完全不用 Green 函數的 Dirichlet 問題存在性證明。核心想法非常簡單：
- 取一個內部球 \(B \subset D\) 與一個外部球 \(B' \supset D\) 。
- 定義一列「上解」 \(u_n\) ：令 \(u_0 = 0\) ，遞迴地令 \(u_{n+1}\) 是 Poisson 方程的解，且在 \(\partial B\) 上等於 \(u_n\) 、在 \(\partial B'\) 上等於邊界函數 \(f\) 的某個上界。
- 證明 \(u_n\) 單調遞增且有上界，故收斂到某個調和函數 \(u\) 。
- 證明 \(u\) 滿足邊界條件。
Perron 定理的深刻之處：它只用到了「球內部有解」這個事實（用 Poisson 核就能構造），完全不需要 Green 函數的先驗存在性。這是 20 世紀 PDE 理論的標準技術：先局部、後整體。在偏微分方程的現代處理裡，Perron 方法更被推廣成粘性解（viscosity solution）的理論：定義一個解只需滿足「每個檢驗函數（上/下解）的測試不等式」，而無需假設二階可微。這套框架是自由邊界問題與 Hamilton–Jacobi 方程的基礎。
習題#
習題 9.2.1（Poisson 方程與球對稱解） 設球半徑 \(R\) 、密度 \(\rho\) 的均質球。 (a) 用 SymPy 證明球內比勢 \(\Phi_{\text{in}} = -2\pi G\rho\left(R^2 - \dfrac{r^2}{3}\right)\) 滿足 \(\nabla^2\Phi = 4\pi G\rho\) ； (b) 證明球外 \(\Phi_{\text{out}} = -\dfrac{GM}{r}\) 滿足 \(\nabla^2\Phi = 0\) ； (c) 驗證 \(\Phi\) 與 \(\Phi'\) 在 \(r = R\) 都連續，並解釋為什麼引力場在球面上不跳躍； (d) 計算球內的平均密度與平均引力：證明球內一點的引力只由該球內的質量決定（牛頓的「球殼定理」），並用 SymPy 比較 \(g(r)\) 與 \(-4\pi G\rho r/3\) 和 \(-G M(r)/r^2\) 。
習題 9.2.2（單極定理與平方反比） (a) 用 SymPy 證明：對任意半徑 \(R\) ，單位點源 \(\Phi = -\frac{1}{4\pi r}\) 的球面通量 \(\oint \frac{\partial\Phi}{\partial n}\mathrm{d}S\) 恆等於 \(1\) ； (b) 用 Green 恆等式證明：只有平方反比的力場才與「無源」相容。提示：設 \(\mathbf{F} = -k/r^{\,p}\,\hat{\mathbf{r}}\) ，寫出 \(\nabla\cdot\mathbf{F}\) 並要求它為零，求 \(p\) ； (c) 解釋為什麼「平方反比」在牛頓的《原理》（1687）裡是經驗歸納，而在 Green（1828）之後變成了定理； (d) 比較地球表面的重力加速度 \(g = 9.81\ \mathrm{m/s^2}\) 與由均質球模型算出的值（取 \(R = 6.37\times10^6\ \mathrm{m}\) ， \(\rho = 5.5\times10^3\ \mathrm{kg/m^3}\) ， \(G = 6.674\times10^{-11}\) ），解釋差異的來源。
習題 9.2.3（Green 恆等式與唯一性） (a) 用 SymPy 在單位正方形上，任選三組 \((u, v)\) （例如 \(u = x^2 - y^2,\ v = xy\) ； \(u = \mathrm{e}^{x}\cos y,\ v = x\) ； \(u = x^3,\ v = y^2\) ），驗證 Green 第二恆等式的左右兩端相等； (b) 用 Green 第一恆等式完整寫出 Dirichlet 問題唯一性的證明，說明每一步用了什麼假設； (c) 構造一個函數 \(w\) 滿足 \(w|_{\partial D} = 0\) 但 \(\nabla^2 w \not\equiv 0\) ，說明為什麼這不構成唯一性定理的反例； (d) 證明 Neumann 問題的相容條件 \(\oint_{\partial D} g\,\mathrm{d}s = 0\) ，並證明解只差一個常數。
習題 9.2.4（最大原理與數值驗證） (a) 用 SymPy 證明強最大原理的「剝削」構造：對調和函數 \(u = x^2 - y^2\) ，找出 \(\varepsilon\) 使得 \(v = u - \varepsilon + \varepsilon\frac{r^2}{a^2}\) 在半徑 \(a\) 的球面上等於 \(u\) ； (b) 用 SOR 求解單位正方形的 Poisson 問題 \(\Delta u = -2\pi^2\sin(\pi x)\sin(\pi y)\) ， \(u = 0\) 於邊界，並用數值實驗檢查最大原理（內部最大值 ≤ 邊界最大值）； (c) 把源項改成 \(\Delta u = 0\) 但邊界值改成 \(f = x(1-x)\) （非零邊界），觀察解的形狀，並驗證解仍滿足最大原理； (d) 比較 \(\omega = 1\) （Gauss–Seidel）與 \(\omega = 1.9\) （SOR）的迭代次數，並說明為什麼 \(\omega > 2\) 會不穩定； (e) 討論：若邊界函數在邊界上有跳變（例如一半是 \(0\) 、一半是 \(1\) ），最大原理還成立嗎？解會怎樣？
習題 9.2.5（Green 函數與 Poisson 核） (a) 用 SymPy 驗證單位圓盤 Green 函數 \(G(z,w) = \frac{1}{2\pi}\log\left|\frac{1-z\bar{w}}{z-w}\right|\) 的三個性質：調和、邊界為零、 \(z\to w\) 時的奇異性 \(-\frac{1}{2\pi}\log|z-w|\) ； (b) 證明 \(G(z,w) = G(w,z)\) （自對稱），並解釋它在 Green 恆等式中的意義； (c) 導出 Poisson 核 \(P_r(\theta) = \frac{1-r^2}{1-2r\cos\theta+r^2}\) ，證明 \(P_r > 0\) 且 \(\int_{-\pi}^{\pi}P_r\,\mathrm{d}\theta = 2\pi\) ； (d) 用 Poisson 核求出邊界值 \(f(\theta) = \cos\theta + \frac{1}{2}\cos 3\theta\) 的解，並與 9.1 節的 Fourier 展開方法（ \(u(r,\theta) = r\cos\theta + \frac{1}{2}r^3\cos 3\theta\) ）對照； (e) 解釋「Poisson 核的守恆性」在物理上意味著什麼，並把它與「 \(P_r(0)\) 在 \(r\to 1\) 時爆炸」一起解釋為「邊界值對靠近邊界的點影響最大」。
小結#
| 概念 | 內容 | SymPy 工具 |
|---|---|---|
| Poisson 方程 | \(\nabla^2\Phi = 4\pi G\rho\) （或 \(\nabla^2\phi = -\rho_e/\varepsilon_0\) ） | sp.diff + sp.simplify |
| 球對稱 Laplace | \(\nabla^2 f = f'' + \frac{2}{r}f'\) | sp.diff |
| Green 第二恆等式 | \(\iint(u\Delta v - v\Delta u) = \oint(u\partial_n v - v\partial_n u)\) | sp.integrate 逐邊計算 |
| Green 第一恆等式 | $\iint | \nabla u |
| 唯一性（能量法） | 兩解之差 \(w\) 滿足 \(\nabla^2w=0\) 、 $w | _{\partial D}=0$ ⇒ \(w\equiv0\) |
| 最大原理 | 調和函數的最大最小值都在邊界上 | 數值檢查內部值 |
| 單極定理 | 點源球面通量恆為常數（Gauss 定律原型） | sp.simplify 化簡通量 |
| 平方反比 | 只有 \(1/r^2\) 力場與無源相容 | sp.diff 求散度 |
| Dirichlet 問題 | 給邊界值求內部調和函數 | SOR 迭代 |
| Neumann 問題 | 相容條件 \(\oint g\,\mathrm{d}s=0\) ；解差一個常數 | — |
| Green 函數 | $G\sim-\frac{1}{2\pi}\log | z-w |
| Poisson 核 | \(P_r(\theta) = \frac{1-r^2}{1-2r\cos\theta+r^2}\) ，正值、守恆 | sp.integrate |
| Perron 定理 | 不用 Green 函數的存在性證明 → 粘性解 | — |
第 9.2 節的三個核心洞見：
- 橢圓型方程只認邊界。這不是人為的約定，而是從「沒有實特徵線」這個代數事實推出來的。Laplace 方程的解由邊界值唯一決定（能量法），最大值與最小值由邊界值控制（最大原理），且能被表示成對邊界資料的一次加權積分（Green 函數）。
- Green 恆等式是把 PDE 與積分綁在一起的橋。從它出發，一次得到：單極定理、唯一性定理、Green 函數、Neumann 相容條件、以及後面第 10 章的 Cauchy 積分公式。
- 勢論是 19 世紀物理的共同語言。牛頓引力、庫侖靜電、Maxwell 電磁場、理想流體的位流，全都服從同一條方程 \(\nabla^2 u = \text{源}\) —— 只是源與常數不同。因為數學結構相同，物理定律之間的相似性才會被看見。
下一節的問題：本章的工具都是手寫出來的。我們要問一個很實際的問題：SymPy 能替我們解 PDE 嗎？ 答案是「能，但非常有限」——sp.pdsolve 只會做一階線性方程；sp.dsolve 對 PDE 直接報錯。這個殘酷的真相本身極有教育意義：它告訴我們，符號求解器能做的，遠少於我們在本書前八章學會的所有東西。但同時，它也是一面鏡子，讓我們看清楚哪些步驟是可以機械化的（分離變數、模式驗證、殘差檢查），哪些步驟必須靠人（套用邊界條件、證明唯一性、判斷解的行為）。
歷史年表（本章相關）
| 年份 | 人物 | 事件 |
|---|---|---|
| 1687 | Newton | 《原理》問世；平方反比定律（經驗歸納） |
| 1785 | Coulomb | 庫侖定律，靜電的平方反比 |
| 1828 | Green | 〈論把數學分析應用於吸引與電的理論〉：Green 恆等式、Green 函數、單極定理 |
| 1830s | Gauss | 位勢理論的基本定理、平均值定理 |
| 1837 | Dirichlet | Dirichlet 問題與「原理」；Lipschitz 的能量法證明唯一性 |
| 1843 | Royal Society | 重印 Green 1828 年論文的英譯本，勢論成為國際共通語言 |
| 1860 | Neumann | Neumann 邊界條件 |
| 1865 | Maxwell | 〈電磁場的動力學理論〉：位能、位移電流，勢論進入電磁學 |
| 1923 | Perron | 不依賴 Green 函數的 Dirichlet 問題存在性定理 |
| 1980s | Crandall–Lions、Ishii | Perron 方法推廣為粘性解（viscosity solution）理論 |