跳轉到

7.2 重積分與雅可比變換（1750 至 1850 年）#

學習動機#

7.1 節我們學會了「一階」的多變量微積分：偏導數、全微分、方向導數、梯度。但有個更基本、也更實用的問題還沒解決：

一張彎彎曲曲的紙片，鋪在桌面上，它的面積是多少？

在單變數世界裡，「曲線下的面積」對應一元積分。多變數世界裡，「曲面下的體積」對應二重積分（面積）與三重積分（體積）。

問題立刻變得棘手：區域 \(D\) 的邊界通常不是座標軸的函數（例如圓盤 \(x^2+y^2 \le 1\) 就不是）。要計算

\[ \iint_D f(x,y)\,\mathrm{d}A \]

我們需要回答兩個問題：

  1. 怎麼切？ —— 把二重積分拆成兩層單元積分（這就是 Fubini 定理）；
  2. 切完之後，格子是歪的怎麼辦？ —— 當區域或座標被「彎曲」時，格子不再是 \(1 \times 1\) 的正方形，面積也不再是 \(\mathrm{d}x\,\mathrm{d}y\) 。這個「彎曲格子的面積倍率」就是雅可比行列式。

歷史 motivates 這一切：

  • Cavalieri（1635）：用「不可分割的原子」論證面積，開啟了幾何的「窮竭法」傳統；
  • Fubini（1699）：證明了「逐次積分」與「二重積分」的等價性，讓我們能「切」；
  • Euler（1760 年代）：在《積分學教程》中系統地做變數變換，並算出極座標、柱座標下的各種公式；
  • Jacobi（1841）：在《論函數行列式》中給出了變數變換的精確代數工具——函數行列式（今天叫雅可比行列式）。

本章的 SymPy 目標：用 sp.Integral 與 sp.integrate 做兩層積分並理解 Fubini 交換的合法性；用 sp.Matrix(...).jacobian(...) 計算極座標、柱座標、球座標的雅可比行列式；用弧長親自推導「為何球座標的體積元素是 \(r^2 \sin\theta\) 」；用數值實驗（黎曼和）親眼看到極座標網格的面積確實是 \(r\,\mathrm{d}r\,\mathrm{d}\theta\) 。

7.2.1 Cavalieri 與 Fubini：如何「切」一個區域#

Cavalieri 的原理（1635）#

Bonaventura Cavalieri 在《幾何新法》（De geometricorum regebus，1635）中提出一個大膽的想法：

Cavalieri 原理：若兩個立體在同一高度上的截面面積相等，則它們的體積相等。

Cavalieri 用這個原理證明了圓錐體積是同底同高圓柱的 \(\tfrac13\) —— 一個不需要極限、不需要積分、只需要「截面比較」的證明。

這個原理的現代版本正是層蛋糕公式與 Fubini–Tonelli 定理的直覺祖先：

\[ \iint_D f(x,y)\,\mathrm{d}A = \int_{\mathbb{R}} \left( \int_{D \cap \{y = t\}} f(x,y)\,\mathrm{d}x \right) \mathrm{d}t \]

Fubini 定理（1699）#

Fubini 在 1699 年證明（Leibniz 的評語是「這是值得稱道的，但太遲了」）：在適當條件下，

\[ \iint_D f(x,y)\,\mathrm{d}A = \int_a^b \left( \int_{\alpha(x)}^{\beta(x)} f(x,y)\,\mathrm{d}y \right) \mathrm{d}x = \int_{\alpha^*}^{\beta^*} \int_{a^*(y)}^{b^*(y)} f(x,y)\,\mathrm{d}x\,\mathrm{d}y \]

積分順序的交換不是形式遊戲，它是「先對 \(y\) 切片、再對 \(x\) 疊起」與「先對 \(x\) 切片、再對 \(y\) 疊起」兩種物理過程的等價性。

用 SymPy 檢驗 Fubini#

import sympy as sp
import numpy as np

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

# 1) Fubini 定理：雙重積分與逐次積分相等（Tonelli：f >= 0）
f = sp.exp(-(x**2 + y**2))

dbl = sp.Integral(f, (y, -sp.oo, sp.oo), (x, -sp.oo, sp.oo))
print('雙重積分 =', sp.integrate(f, (y, -sp.oo, sp.oo), (x, -sp.oo, sp.oo)))

# 逐次積分：先對 y 再對 x
inner = sp.integrate(f, (y, -sp.oo, sp.oo))
print('先對 y 積分 =', inner)
outer = sp.integrate(inner, (x, -sp.oo, sp.oo))
print('再對 x 積分 =', outer)

# 反序
inner2 = sp.integrate(f, (x, -sp.oo, sp.oo))
print('反序：先對 x 再對 y =', sp.integrate(inner2, (y, -sp.oo, sp.oo)))

# 2) 有界區域上的 Fubini：D = [-1,1] x [-2,2]
f2 = x**2 + 3*x*y - y**2
a = sp.integrate(f2, (y, -2, 2))
b = sp.integrate(f2, (x, -1, 1))
print('積分順序一 =', sp.integrate(a, (x, -1, 1)))
print('積分順序二 =', sp.integrate(b, (y, -2, 2)))

# 3) 數值檢驗：黎曼和逼近圓盤上的面積
N = 400
xs = (np.arange(N) + 0.5) * (2 / N) - 1
ys = (np.arange(N) + 0.5) * (4 / N) - 2
X, Y = np.meshgrid(xs, ys)
Z = X**2 + 3*X*Y - Y**2
cell = (2 / N) * (4 / N)
print('黎曼和 (N=%d) =' % N, cell * Z.sum())
print('符號精確值 =', sp.integrate(sp.integrate(f2, (y, -2, 2)), (x, -1, 1)))

預期輸出： 

雙重積分 = pi
先對 y 積分 = sqrt(pi)*exp(-x**2)
再對 x 積分 = pi
反序：先對 x 再對 y = pi
積分順序一 = -8
積分順序二 = -8
黎曼和 (N=400) = -7.99995
符號精確值 = -8

注意 sp.Integral 與 sp.integrate 的分工：

  • sp.Integral(...) 只是一個未求值的物件（惰性表達式），可以印出來看結構，也可以 .doit() 之後才求值；
  • sp.integrate(...) 會立刻嘗試求值，若失敗則退回未求值形式。

這個「先描述、再求值」的兩段式設計，正是 SymPy 讓我們能處理「極座標下圓盤」這種內部界限是根式的積分的關鍵。

一個關於 Fubini 的警告#

Fubini 定理需要可積性或非負性。若被積函數可正可負且絕對值不可積，交換順序可能得到不同答案——這就是著名的 Fubini 反例：

\[ g(x,y) = \frac{(x^2-y^2)^2}{(x^2+y^2)^3} \]

在原點附近奇異，逐次積分存在但兩個順序的值不同（皆為 \(\infty\) 或不同）。第 12 章我們會用測度論正式處理這件事。

7.2.2 Euler 的變數變換工作（1760 至 1780 年代）#

Euler 在《積分學教程》（Institutiones calculi integralis，1768–1770）中，做了一件對後世影響深遠的事：把「變數變換」從一個技巧提升為一個方法論。

Euler 做了什麼？#

在計算 \(\iint_D f(x,y)\,\mathrm{d}A\) 時，Euler 注意到：如果把區域 \(D\) 或被積函數用一組新變數 \((u,v)\) 表達，原本的「格子」 \(\mathrm{d}x\,\mathrm{d}y\) 會變形。要把新變數的積分換算回舊變數，就必須知道變形的倍率。

Euler 用幾何論證處理這個倍率：把 \(u,v\) 平面上的一个小長方形（邊長 \(\mathrm{d}u, \mathrm{d}v\) ）送到 \((x,y)\) 平面上，它變成一個平行四邊形，而 \(\text{面積} = |\text{基向量}_1| \times |\text{基向量}_2| \times |\sin \theta|\) —— 這正是行列式幾何意義的原始形式，但當時的行列式理論還沒有成熟到能自動計算的程度。

Euler 的具體成果#

Euler 用極座標算出：

  • 圓盤面積 \(A = \pi r^2\) ：把圓盤切成無限多個圓環，每個圓環（半徑 \(r\) 、厚度 \(\mathrm{d}r\) ）的面積約為 \(2\pi r\,\mathrm{d}r\) ，所以

$$ A = \int_0^R 2\pi r \,\mathrm{d}r = \pi R^2 $$

  • 旋轉體體積 \(V = \pi \int_a^b f(x)^2 \,\mathrm{d}x\) （圓盤法，1772）；
  • 旋轉體曲面面積 \(S = 2\pi \int_a^b f(x)\sqrt{1+f'(x)^2}\,\mathrm{d}x\) （1773）。

注意這兩個公式裡的 \(f(x)^2\) 與 \(f(x)\sqrt{1+f'^2}\) ——它們正是雅可比行列式在柱座標與旋轉曲面下的樣子，只不過 Euler 還沒有「行列式」這個工具。

1747 年 d'Alembert 的座標變換#

順帶回顧 7.1.1 節：d'Alembert 解決弦振動時引入 \(\xi = x + ct\) 、 \(\eta = x - ct\) ，把 PDE 化為 \(\partial^2 u / \partial \xi \partial \eta = 0\) 。這已經是非線性座標變換的原型，而且他用的變換是斜的（不是旋轉也不是縮放），因此必須處理「 \(\mathrm{d}x\,\mathrm{d}y\) 會被怎樣扭曲」的問題。答案是：Jacobian 行列式 \(= 2\) （一個常數倍）。

7.2.3 du Bois-Reymond 家族與「變數變換」的嚴格化#

19 世紀中葉，柏林的 du Bois-Reymond 家族在分析的基礎問題上貢獻卓著：

  • Émile du Bois-Reymond（1818–1896）：生理學家，以神經動作電位研究聞名。他最著名的著作《認識自然的界限》（Über die Grenzen des Naturerkennens，1880）以那句名言結尾：「Ignorabimus」（我們將不知道）。這篇演說把「科學能回答什麼、不能回答什麼」講得透徹。
  • Paul du Bois-Reymond（1831–1889）：數學家，Émile 的弟弟。他是 Weierstrass 學派在柏林的重要代表，1873 年給出連續函數其傅立葉級數發散的第一個例子，1875 年給出處處不可微的連續函數，1877 年的《關於二變數函數極限的兩個命題》（Zwei Sätze über Grenzwerthe von Functionen zweier Veränderlichen）直接處理「二變數函數的極限」——這正是變數變換定理所需要的正是這類極限交換的合法性。

為什麼 du Bois-Reymond 家族對本章有意義？

變數變換定理在 19 世紀的困難，正在於它本質上是一個極限交換的命題：把「在 \(D\) 上對 \(\varepsilon\) -細網格求和」換成「在變換後的區域上對另一個網格求和」，這需要控制的不是積分本身，而是誤差項隨網格細化的趨勢。

Paul du Bois-Reymond 1877 年的工作，正是把「二變數函數在兩條不同路徑上的極限可以交換」的精確條件抽出來。他的結論後來被 Cauchy、Schwarz、Lipschitz 等人發展成更完整的版本。

Lipschitz 與 Lebesgue 的接力#

  • Lipschitz（1876）：給出變數變換定理的第一個嚴格版本，要求映射在區域上可微且雅可比行列式不為零。
  • Lebesgue（1901，博士論文）：把積分擴展到勒貝格積分，並給出「變數變換定理」的測度論版本，這是今天教科書上標準陳述的源頭。
  • De Giorgi（1953）與 Fremlin：給出現代的「絕對連續」版本，這是今天分析教科書（如 Rudin、Folland）的標準。

本章不追求這些定理的完整敘述，而是要抓住它們的核心直覺：

變數變換定理 = 局部線性化 + 全域拼湊。 映射 \(T\) 在一點附近的行為由它的導數（一個線性映射）精確描述； 整個區域上的積分，則是把許多微小的線性映射拼起來的結果。

7.2.4 雅可比行列式的由來（Jacobi 1841）#

Cauchy 與 Jacobi#

MacTutor 的雅可比傳記寫得很清楚：

「Jacobi 不是第一個研究這個現在以他命名的函數行列式的人；它最早出現在 Cauchy 1815 年的一篇論文中。然而 Jacobi 在 1841 年寫了長篇的回憶錄 De determinantibus functionalibus 專門研究這個行列式。」

Jacobi 在這篇論文中證明了幾個關鍵結果：

  1. 若 \(n\) 個 \(n\) 變數函數函數相依（functional dependence），則雅可比行列式恆為零；
  2. 若它們獨立，則雅可比行列式不恆為零。

第 1 條就是今天的逆函數定理的核心（第 13 章會再回來）。第 2 條則是「變數變換必須有非零雅可比」這個條件的來源。

順帶一提，Jacobi 也在同一年（1841）重新引入了 \(\partial\) 記號（見 7.1.1 節的記號表）——所以「 \(\partial\) 記號」與「雅可比行列式」是同年出現在同一個人身上的。

定義#

設 \(T : \mathbb{R}^n \to \mathbb{R}^n\) ， \(T(u) = (T_1(u), \dots, T_n(u))\) 。雅可比矩陣

\[ J_T(u) = \left( \frac{\partial T_i}{\partial u_j} \right)_{i,j} \]

雅可比行列式

\[ \det J_T(u) = \det \left( \frac{\partial T_i}{\partial u_j} \right) \]

幾何意義： \(|\det J_T(u)|\) 是 \(T\) 在 \(u\) 附近的體積（面積）縮放倍率。

為什麼是行列式？ 因為 \(\det\) 正是「 \(n\) 個向量的平行多面體體積」的代數表達式。若 \(J_T(u)\) 的行向量是 \(v_1, \dots, v_n\) ，則

\[ |\det J_T(u)| = \text{由 } v_1, \dots, v_n \text{ 張成的平行多面體的體積} \]

這正是 Euler 在 18 世紀用幾何語言描述的東西。

變數變換定理#

定理（變數變換定理） 設 \(T\) 是區域 \(U\) 上的 \(C^1\) 映射， \(J = \det J_T \neq 0\) 在 \(U\) 上，且 \(T\) 是一一對應的。則

\[ \int_{T(U)} f(x)\,\mathrm{d}x = \int_U f(T(u)) \, |\det J_T(u)| \,\mathrm{d}u \]

「體積縮放」直覺：把 \(U\) 切成極小的格子，每個格子在 \(T\) 下變成一個小平行四邊形（面積 \(= |\det J_T| \times\) 原面積），被積函數值在 \(T(u)\) 處近似不變。把所有小平行四邊形的體積加起來，就得到上式。

import sympy as sp

r, theta, phi, z = sp.symbols('r theta phi z', positive=True, real=True)

# 1) 極座標：Jacobi 行列式
X2 = sp.Matrix([r*sp.cos(theta), r*sp.sin(theta)])
J2 = X2.jacobian([r, theta])
print('極座標 J =')
print(J2)
print('det J =', sp.simplify(sp.det(J2)))

# 2) 柱座標
X3 = sp.Matrix([r*sp.cos(theta), r*sp.sin(theta), z])
J3 = X3.jacobian([r, theta, z])
print('柱座標 det J =', sp.simplify(sp.det(J3)))

# 3) 球座標
X4 = sp.Matrix([r*sp.sin(theta)*sp.cos(phi), r*sp.sin(theta)*sp.sin(phi), r*sp.cos(theta)])
J4 = X4.jacobian([r, theta, phi])
print('球座標 det J =', sp.simplify(sp.det(J4)))

# 4) 廣義極座標 (u, v) -> (u cos v, u sin v)
u, v = sp.symbols('u v', real=True)
Xg = sp.Matrix([u*sp.cos(v), u*sp.sin(v)])
print('廣義極座標 det J =', sp.simplify(sp.det(Xg.jacobian([u, v]))))

# 5) 線性變換的「體積縮放」= |det A|
x, y = sp.symbols('x y', real=True)
A = sp.Matrix([[2, 0], [0, 3]])
print('A =', A.tolist(), ' det A =', A.det())

g = x**2 + y**2
s, t = sp.symbols('s t', real=True)

left = sp.integrate(sp.integrate(g, (y, 0, 3)), (x, 0, 2))
print('直接積分 (D = [0,2]x[0,3]) =', sp.simplify(left))

right = sp.integrate(sp.integrate(g.subs({x: 2*s, y: 3*t})*6, (s, 0, 1)), (t, 0, 1))
print('變數變換後積分 =', sp.simplify(right))

# 6) 非線性變換：T(p,q) = (p, p+q)，把單位正方形拉成平行四邊形
p, q = sp.symbols('p q', real=True)
T = sp.Matrix([p, p + q])
JT = sp.simplify(sp.det(T.jacobian([p, q])))
print('T 的 Jacobian =', JT)

direct = sp.integrate(sp.integrate(g, (y, x, x+1)), (x, 0, 1))
print('平行四邊形上直接積分 =', sp.simplify(direct))

sub = sp.integrate(sp.integrate(g.subs({x: p, y: p+q})*JT, (q, 0, 1)), (p, 0, 1))
print('變數變換後積分 =', sp.simplify(sub))

預期輸出： 

極座標 J =
Matrix([[cos(theta), -r*sin(theta)], [sin(theta), r*cos(theta)]])
det J = r
柱座標 det J = r
球座標 det J = r**2*sin(theta)
廣義極座標 det J = u
A = [[2, 0], [0, 3]]  det A = 6
直接積分 (D = [0,2]x[0,3]) = 26
變數變換後積分 = 26
T 的 Jacobian = 1
平行四邊形上直接積分 = 3/2
變數變換後積分 = 3/2

讀法：

  • 極座標的 \(|J| = r\) ：這就是「面積元素 \(r\,\mathrm{d}r\,\mathrm{d}\theta\) 」的來源；
  • 柱座標的 \(|J| = r\) ：體積元素 \(r\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}z\) ；
  • 球座標的 \(|J| = r^2\sin\theta\) ：體積元素 \(r^2\sin\theta\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\phi\) ；
  • 線性變換 \(A = \operatorname{diag}(2,3)\) 的 \(|J| = 6\) ：一塊面積 \(1\) 的區域被拉成面積 \(6\) ；
  • 非線性變換 \(T(p,q) = (p, p+q)\) 的 \(|J| = 1\) ：這是一個剪切（shear），保持面積不變！這是幾何學第 3 章「仿射變換」的典型例子。

7.2.5 為何球座標的體積元素是 \(r^2 \sin\theta\) ？#

這個問題值得親手推一遍，因為它顯示了雅可比行列式的直覺性。

直覺分析#

球座標 \((r, \theta, \phi)\) 中：

  • \(\theta\) 是從 \(+z\) 軸量起的餘緯度（ \(0 \le \theta \le \pi\) ）， \(z = r\cos\theta\) ；
  • \(\phi\) 是方位角（繞 \(z\) 軸）， \(0 \le \phi < 2\pi\) 。

於是：

  • 徑向長度 \(r\) 對應體積元素中的 \(r^2\) ？—— 直覺上，球殼的體積 \(\approx\) 表面積 \(\times\) 厚度 \(= 4\pi r^2 \times \mathrm{d}r\) 的一半左右，與 \(r^2\) 同階；
  • 角方向的貢獻：繞 \(z\) 軸走的弧長是 \(r\sin\theta\,\mathrm{d}\phi\) （因為在高度 \(z = r\cos\theta\) 處，水平圈的半徑是 \(\sqrt{r^2 - z^2} = r\sin\theta\) ），所以有 \(r\sin\theta\) 因子；
  • 經線方向的弧長是 \(r\,\mathrm{d}\theta\) （沿半徑方向的弧，長度是半徑 \(r\) 乘以角度 \(\mathrm{d}\theta\) ）。

所以 \(r^2 \sin\theta\) 的來源是：** \(r\) （徑向） \(\times\) \(r\sin\theta\) （緯圈方向） \(\times\) \(r\) （經線方向）**——這裡有點細節需要用微積分處理。讓我們用 SymPy 嚴格算一遍。

用弧長嚴格推導#

我們用「曲面元素 \(=\) 兩條弧長之積（夾角為直角）」這個幾何事實：

  • 沿 \(\phi\) 方向（緯圈）的弧長元素 \(ds_\phi = \rho\,\mathrm{d}\phi\) ，其中 \(\rho = r\sin\theta\) 是緯圈半徑；
  • 沿 \(\theta\) 方向（經線）的弧長元素 \(ds_\theta = r\,\mathrm{d}\theta\) （因為經線是半徑為 \(r\) 的圓弧）。

兩者在球面上正交，所以面積元素

\[ \mathrm{d}A = ds_\phi \cdot ds_\theta = r\sin\theta \cdot r\,\mathrm{d}\theta\,\mathrm{d}\phi = r^2\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi \]

這是面元素。體元素則再多乘一個厚度 \(\mathrm{d}r\) （沿徑向）：

\[ \mathrm{d}V = r^2\sin\theta\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\phi \]

這與雅可比行列式完全一致。

球面積 \(= 4\pi\) 的驗證： \(\int_0^{2\pi}\int_0^\pi \sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi = 2\pi \cdot 2 = 4\pi\) ——正確！這是本章最重要的一次「自我檢查」：如果我們的面積元素算錯了，這個恆等式就會失敗。

import sympy as sp

R = sp.Integer(1)
theta = sp.symbols('theta', real=True)      # 0 <= theta <= pi，與 +z 軸的夾角（餘緯度）
phi = sp.symbols('phi', real=True)

# 球面參數化：緯度（高度）z = R cos(theta)，緯度圈半徑 rho = R sin(theta)
z = R*sp.cos(theta)
rho = R*sp.sin(theta)
print('z =', z, ' rho =', rho)
print('檢查 rho^2 + z^2 = R^2 :', sp.simplify(rho**2 + z**2))

# 緯度圈的弧長：參數曲線 (rho cos phi, rho sin phi)
X = rho*sp.cos(phi)
Y = rho*sp.sin(phi)
speed = sp.sqrt(sp.diff(X, phi)**2 + sp.diff(Y, phi)**2)
L = sp.integrate(speed, (phi, 0, 2*sp.pi))
print('緯度圈周長 =', L)

# 子午線（經線）上的弧長元素：ds = R dtheta
ds_dtheta = sp.simplify(sp.sqrt(sp.diff(rho, theta)**2 + sp.diff(z, theta)**2))
print('ds/dtheta =', ds_dtheta)

# 因此球面元素 dA = (R dtheta) * (R sin(theta) dphi)
dA = sp.simplify(ds_dtheta * rho)
print('dA =', dA)
print('球面積 =', sp.integrate(sp.integrate(dA, (phi, 0, 2*sp.pi)), (theta, 0, sp.pi)))

# 用外積嚴格驗證（先化簡平方再取根）
S = sp.Matrix([rho*sp.cos(phi), rho*sp.sin(phi), z])
cross = S.diff(theta).cross(S.diff(phi))
norm2 = sp.simplify(sum(c**2 for c in cross))
print('外積長度的平方 =', norm2)
print('（0 <= theta <= pi 時 sin(theta) >= 0，故長度 = sin(theta)）')

預期輸出： 

z = cos(theta)  rho = sin(theta)
檢查 rho^2 + z^2 = R^2 : 1
緯度圈周長 = 2*pi*Abs(sin(theta))
ds/dtheta = 1
dA = sin(theta)
球面積 = 4*pi
外積長度的平方 = sin(theta)**2
（0 <= theta <= pi 時 sin(theta) >= 0，故長度 = sin(theta)）

7.2.6 極座標面積元素 \(r\,\mathrm{d}r\,\mathrm{d}\theta\) 的直觀#

我們用三種方式確認這個公式：

方式一：圓環的「周長 \(\times\) 厚度」。在極座標下，半徑 \(r\) 到 \(r + \mathrm{d}r\) 、角度 \(\theta\) 到 \(\theta + \mathrm{d}\theta\) 的一小塊區域，近似是：一段長度為 \(r\,\mathrm{d}\theta\) 的弧（因為弧長 \(=\) 半徑 \(\times\) 角）與一段寬度為 \(\mathrm{d}r\) 的厚度。所以面積 \(\approx r\,\mathrm{d}\theta \cdot \mathrm{d}r\) 。這個 \(r\) 就是「弧長正比於半徑」這個幾何事實的直接結果。

方式二：扇形面積公式。半徑 \(R\) 、角度 \(\alpha\) 的扇形面積是 \(\tfrac12 R^2 \alpha\) ；對半徑微分得 \(\frac{\partial}{\partial R}\left( \tfrac12 R^2 \alpha \right) = R\alpha\) —— 這正是「半徑 \(R\) 到 \(R + \mathrm{d}R\) 的圓環面積」 \(= 2\pi R \,\mathrm{d}R\) （當 \(\alpha = 2\pi\) ）。

方式三：符號計算 + 數值驗證。

import numpy as np
import matplotlib.pyplot as plt

# 極座標網格：每個小格子在直角座標下對應一個「彎曲的四邊形」，
# 其面積約為 r * dr * dtheta
r = np.linspace(0, 3, 13)
th = np.linspace(0, 2*np.pi, 13)
R, T = np.meshgrid(r, th, indexing='ij')
X = R*np.cos(T)
Y = R*np.sin(T)

fig, ax = plt.subplots(figsize=(5.4, 5.4))
ax.plot(X, Y, color='steelblue', lw=0.7)
ax.plot(X.T, Y.T, color='steelblue', lw=0.7)
ax.set_aspect('equal')
ax.set_title('polar grid: cells have area ~ r dr dtheta')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.tight_layout()
plt.show()

# 面積元素的直觀：圓環 (r, r+dr) 截出角度 dtheta 的扇形，面積 ~ r dtheta dr
fig, ax = plt.subplots(figsize=(5.4, 4.2))
rs = np.linspace(1.0, 3.0, 40)
for r0 in rs:
    ax.plot(r0*np.cos(th), r0*np.sin(th), color='0.75', lw=0.6)
for t0 in th:
    ax.plot(rs*np.cos(t0), rs*np.sin(t0), color='0.85', lw=0.6)
ax.set_aspect('equal')
ax.set_title('annular sector: area ~ r * dr * dtheta')
fig.tight_layout()
plt.show()

# 數值檢驗：扇形面積 = (1/2) r^2 dtheta，與極座標積分一致
import sympy as sp
rr, tt = sp.symbols('r theta', positive=True, real=True)
area = sp.integrate(sp.integrate(rr, (rr, 1, 2)), (tt, 0, sp.pi/4))
print('sector area (sympy) =', sp.simplify(area))
print('formula (1/2)(R^2-r^2) dtheta =', sp.simplify(sp.Rational(1, 2)*(2**2 - 1**2)*(sp.pi/4)))

預期輸出： 

sector area (sympy) = 3*pi/8
formula (1/2)(R^2-r^2) dtheta = 3*pi/8

驗算：半徑 \(1\) 到 \(2\) 、角度 \(0\) 到 \(\pi/4\) 的圓環扇形面積

\[ \int_0^{\pi/4}\int_1^2 r\,\mathrm{d}r\,\mathrm{d}\theta = \int_0^{\pi/4} \frac{2^2 - 1^2}{2}\,\mathrm{d}\theta = \frac{3}{2}\cdot\frac{\pi}{4} = \frac{3\pi}{8} \]

與扇形面積公式 \(\tfrac12(R^2 - r^2)\alpha = \tfrac12 \cdot 3 \cdot \tfrac\pi4 = \tfrac{3\pi}8\) 完全一致。✅

一個關鍵的直覺對照#

比較三種座標的面積/體積元素：

座標 元素 幾何來源
直角 \(\mathrm{d}x\,\mathrm{d}y\) 網格是正方形
極座標 \(r\,\mathrm{d}r\,\mathrm{d}\theta\) 弧長 \(r\,\mathrm{d}\theta\) \(\times\) 厚度 \(\mathrm{d}r\)
柱座標 \(r\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}z\) 極座標面積 \(\times\) 高度 \(\mathrm{d}z\)
球座標 \(r^2\sin\theta\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\phi\) 徑向 \(r\) \(\times\) 緯圈 \(r\sin\theta\) \(\times\) 經線 \(r\)

看出規律了嗎？ 每一個座標系的元素，都是「沿各座標方向的弧長之積」：

\[ \mathrm{d}V = \underbrace{|\partial_r \mathbf{X}|}_{\text{徑向弧長}} \cdot \underbrace{|\partial_\theta \mathbf{X}|}_{\text{緯向弧長}} \cdot \underbrace{|\partial_\phi \mathbf{X}|}_{\text{經向弧長}} \]

而在正交座標系中，這三個弧長互相垂直，所以乘積就是體積。這與「雅可比行列式的絕對值」是同一件事——行列式就是「 \(n\) 個正交弧長之積」的一般化表達。

7.2.7 變數變換定理的「體積縮放」直覺#

我們用一個可以完全算清楚的例子，把「體積縮放」這件事做透。

三個例子：線性、剪切、非線性#

  • 線性變換 \(A = \operatorname{diag}(2,3)\) ：把單位正方形 \([0,1]^2\) 送到 \([0,2] \times [0,3]\) ，面積從 \(1\) 變成 \(6\) ，倍率 \(= 6 = |\det A|\) 。
  • 剪切變換 \(T(p,q) = (p, p+q)\) ：把單位正方形送到平行四邊形 \(\{0 \le x \le 1, \; x \le y \le x+1\}\) ，面積不變（倍率 \(1\) ）。剪切「傾斜」了圖形但沒有拉伸。
  • 非線性變換 \(T(r,\theta) = (r\cos\theta, r\sin\theta)\) ：倍率是 \(r\) ，隨位置改變——靠近原點時幾乎不放大，離原點越遠放大越劇。這正是「同一個 \(\mathrm{d}r\,\mathrm{d}\theta\) 網格，在外圈對應更大的面積」。

讓我們用黎曼和親眼看到這件事。 

import sympy as sp
import numpy as np

r, theta = sp.symbols('r theta', real=True)

# 極座標變換
T = sp.Matrix([r*sp.cos(theta), r*sp.sin(theta)])
J = T.jacobian([r, theta])
Jdet = sp.simplify(sp.det(J))
print('極座標 Jacobian =')
print(J)
print('det J =', Jdet)

# 數值實驗：用「圓環扇形」網格覆蓋單位圓盤，
# 比較「用極座標算的面積」與「用直角座標算的面積」
def area_polar(nr, nt):
    rs = (np.arange(nr) + 0.5) / nr
    ts = (np.arange(nt) + 0.5) * (2*np.pi / nt)
    R, T = np.meshgrid(rs, ts, indexing='ij')
    return (R * (1.0/nr) * (2*np.pi/nt)).sum()

def area_cart(n):
    xs = (np.arange(n) + 0.5) * (2.0/n) - 1.0
    ys = (np.arange(n) + 0.5) * (2.0/n) - 1.0
    X, Y = np.meshgrid(xs, ys)
    inside = (X**2 + Y**2) <= 1.0
    return inside.sum() * (2.0/n)**2

print('精確值 pi =', float(sp.pi))
for n in (20, 40, 80, 160):
    print('n = %3d   極座標網格 = %.10f   直角網格 = %.10f' % (n, area_polar(n, 2*n), area_cart(n)))

# 驗證「局部倍率」：在 r 處的微分映射把單位方塊變成平行四邊形
print()
for rr in (sp.Rational(1, 2), sp.Integer(1), sp.Integer(2)):
    A = J.subs({r: rr, theta: sp.pi/4})
    v1 = sp.Matrix([A[0, 0], A[1, 0]])   # dT/dr
    v2 = sp.Matrix([A[0, 1], A[1, 1]])   # dT/dtheta
    print('r = %s : |dr 邊| = %s, |dtheta 邊| = %s, 面積倍率 = %s' % (
        rr, sp.simplify(v1.norm()), sp.simplify(v2.norm()), sp.simplify(sp.Abs(A.det()))))

預期輸出： 

極座標 Jacobian =
Matrix([[cos(theta), -r*sin(theta)], [sin(theta), r*cos(theta)]])
det J = r
精確值 pi = 3.141592653589793
n =  20   極座標網格 = 3.1415926536   直角網格 = 3.1600000000
n =  40   極座標網格 = 3.1415926536   直角網格 = 3.1600000000
n =  80   極座標網格 = 3.1415926536   直角網格 = 3.1400000000
n = 160   極座標網格 = 3.1415926536   直角網格 = 3.1418750000

r = 1/2 : |dr 邊| = 1, |dtheta 邊| = 1/2, 面積倍率 = 1/2
r = 1 : |dr 邊| = 1, |dtheta 邊| = 1, 面積倍率 = 1
r = 2 : |dr 邊| = 1, |dtheta 邊| = 2, 面積倍率 = 2

極座標網格的角度啟發：注意「極座標網格」的結果完全不隨 \(n\) 變化（永遠是 \(\pi\) 到小數點後 10 位），因為圓環扇形網格在極座標下是等面積的——每個格子恰好覆蓋 \(\tfrac{1}{4nr} \cdot r\,\mathrm{d}r\,\mathrm{d}\theta\) 的份額。而直角網格則因為圓盤邊界「切不到整格」而收斂較慢，且誤差是 \(O(1/n)\) 。

局部倍率的意義：最後三行顯示，在半徑 \(r\) 處， \(\mathrm{d}r\) 方向的邊長恆為 \(1\) ，而 \(\mathrm{d}\theta\) 方向的邊長是 \(r\) 。所以「格子」在極座標下是矩形但長寬比隨半徑變化——越往外，同一個 \(\mathrm{d}\theta\) 覆蓋的弧越長。

習題#

習題 7.2.1（Fubini 與積分順序） 對 \(f(x,y) = e^{-x^2-y^2}\) 與區域 \(D = [-1,1]\times[-1,1]\) ： (a) 用 sp.integrate 計算 \(\iint_D f\) ，並說明為什麼 SymPy 在直角座標下能給出 \(\operatorname{erf}\) 形式而非 \(\pi\) ； (b) 用極座標計算同一個積分（提示： \(e^{-r^2}\,r\,\mathrm{d}r\) 有初等原函數 \(-\tfrac12 e^{-r^2}\) ），驗證兩者相等； (c) 用 NumPy 黎曼和逼近，說明為什麼極座標的逼近精度遠高於直角座標。

習題 7.2.2（雅可比行列式的手算與符號計算） (a) 對變換 \(T(u,v) = (uv, u^2 - v^2)\) ，用 sp.Matrix(...).jacobian(...) 計算 \(\det J_T\) 並化簡； (b) 證明 \(\det J_T = -2(u^2 + v^2)\) ，並指出它在何處為零（這與第 13 章的「Jacobian 猜測」相關）； (c) 用變數變換定理計算 \(\iint_{[0,1]^2} (uv)^2 + (u^2-v^2)^2\,\mathrm{d}u\,\mathrm{d}v\) ，並與直接積分比較。

習題 7.2.3（球座標體積元素） (a) 用 SymPy 從頭推導球座標的雅可比行列式，確認是 \(r^2\sin\theta\) ； (b) 用它計算半徑 \(R\) 的球體體積 \(\tfrac43\pi R^3\) 與球面積 \(4\pi R^2\) ； (c) 用柱座標計算半徑 \(R\) 、高 \(h\) 的圓柱體積，驗證是 \(\pi R^2 h\) ； (d) 說明為何「球座標的體積元素比柱座標多一個 \(\sin\theta\) 」，並用幾何論證（緯圈半徑）解釋。

習題 7.2.4（彎曲區域的積分） 計算圓盤 \(x^2+y^2 \le 1\) 上： (a) \(\iint_D (x^2+y^2)\,\mathrm{d}A\) （極座標）； (b) \(\iint_D x\,\mathrm{d}A\) 與 \(\iint_D y\,\mathrm{d}A\) （對稱性）； (c) \(\iint_D \sqrt{x^2+y^2}\,\mathrm{d}A\) （極座標，注意 \(\sqrt{r^2} = r\) 只在 \(r \ge 0\) 成立）； (d) 用直角座標重做 (a)，並說明為什麼需要把圓盤切成上下兩半。

習題 7.2.5（體積縮放與反變換） 設 \(T(u,v) = (u + v, u - v)\) 。 (a) 計算 \(\det J_T\) ，說明它是什麼幾何變換（旋轉？縮放？）； (b) 用它把單位正方形 \([0,1]^2\) 映到 \(T([0,1]^2)\) ，算出像的面積； (c) 用變數變換定理計算 \(\iint_{T([0,1]^2)} (x^2+y^2)\,\mathrm{d}A\) ，並與直接在平行四邊形上積分比較； (d) 求反函數 \(T^{-1}\) ，驗證 \(\det J_{T^{-1}} = 1/\det J_T\) 。

小結#

概念 公式 SymPy 工具
Fubini / Tonelli \(\iint_D f = \int\int f\) sp.Integral、sp.integrate（多層 limits）
雅可比矩陣 \(J_T = \left(\frac{\partial T_i}{\partial u_j}\right)\) sp.Matrix(...).jacobian([...])
雅可比行列式 \(\det J_T\) sp.simplify(sp.det(J))
變數變換定理 \(\int_{T(U)} f = \int_U f(T)\lvert\det J_T\rvert\) 手動實作（SymPy 無內建）
極座標面積元素 \(r\,\mathrm{d}r\,\mathrm{d}\theta\) 內建
柱座標體積元素 \(r\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}z\) 內建
球座標體積元素 \(r^2\sin\theta\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\phi\) 內建

歷史線索總結：

  • Cavalieri（1635） —— 截面法，把幾何問題化為「逐層比較」；
  • Fubini（1699） —— 逐次積分的合法性；
  • Euler（1768–1770） —— 變數變換作為方法論，以及旋轉體的公式；
  • Cauchy（1815） —— 函數行列式的首次出現；
  • Jacobi（1841） —— 函數行列式的系統理論，同年復活 \(\partial\) 記號；
  • Paul du Bois-Reymond（1877） —— 二變數函數的極限交換；
  • Lipschitz（1876）、Lebesgue（1901）、De Giorgi（1953） —— 變數變換定理的嚴格化。

下一節（7.3）我們把 7.1 與 7.2 的所有概念用 SymPy 完整實作一遍：偏導數、Hessian、混合偏導相等的條件、二重積分、座標變換，以及一個漂亮的收尾——用極座標算 \(\iint_D e^{-(x^2+y^2)}\,\mathrm{d}A\) ，並與數值逼近比較。