13.1 變分法、梯度下降與最佳化#
學習動機#
前面十二章，我們處理的都是固定的對象：給定一個函數 \(f\) ，求它的極限、導數、積分、級數。但有另一整類問題是這樣的：
給定一個函數空間，在其中找出使某個「泛函」 \(J[y]\) 取到極值的函數。
泛函（functional）是把函數空間映到數的映射。最經典的例子：
- 最短曲線問題： \(J[y] = \int_{a}^{b} \sqrt{1 + y'(x)^2}\,\mathrm{d}x\) 在給定端點條件下取最小值（答案是直線）；
- 懸鏈線問題： \(J[y] = \int_{a}^{b} y(x)\sqrt{1 + y'(x)^2}\,\mathrm{d}x\) （均勻鏈條的重力位能）取最小值（答案是懸鏈線）；
- brachistochrone 問題： \(J[y] = \int \frac{\sqrt{1+y'^2}}{\sqrt{y}}\,\mathrm{d}x\) （滑落的總時間）取最小值（答案是擺線）；
- 最小作用量原理（Hamilton 1834）： \(J[q] = \int L(q, \dot q, t)\,\mathrm{d}t\) 在真實運動軌跡上取極值。
這些問題在十八、十九世紀是物理與幾何的真正難題，因為人沒有「對函數求導」的工具。第 3 章我們學過 \(\frac{\mathrm{d}f}{\mathrm{d}x}\) ，它是對數值求導；現在我們需要一個新的算子——對函數求導，即變分（variation）。
本節的任務有四個：
- 弄清楚問題的來源：最小作用量原理、最短曲線／等周問題，是怎樣從物理中長出來的；
- 建立 Euler–Lagrange 方程：Euler 在 1752 年如何從力學導出它，Jacobi 在 1845 年如何用「最後乘法因子」把它推到更一般的形式；
- 看見嚴格化：Legendre 與 Weierstrass 如何把「必要條件」升級為真正的定理；
- 連到現代：變分法與今天的梯度下降、隨機梯度下降（Robbins–Monro 1951）、機器學習是同一件事的不同外衣。
最後我們要用 SymPy 實際推出 Euler–Lagrange 方程，並且用數值梯度下降離散化同一個泛函，看兩條路是否收斂到同一個答案。這是全書唯一一次「歷史 → 符號 → 數值」三路合流的示範。
13.1.1 問題的來源：物理如何逼出「對函數求導」#
(a) 最小作用量原理#
牛頓在 1687 年的《自然哲學的數學原理》裡給出的力學定律是「加速度正比於力」。這個表述好用，但有一個缺陷：它是一階微分方程，給定初始位置與速度就能解出運動。但物理學家常常想知道的是相反的問題：已知運動軌跡，問「力是什麼？」
Hamilton 在 1834 年（更早的 Maupertuis 在 1744 年已提出雛形）給出了激動人心的答案：
最小作用量原理：物體在兩個固定時刻 \(t_1, t_2\) 之間的實際軌跡 \(q(t)\) ，使作用量
取到駐值（stationary value）。
這是一個驚人的「反演」：牛頓的力 \(F = -\frac{\partial L}{\partial q}\) 只是 Euler–Lagrange 方程在 \(t\) 依賴顯式時的一個特例。
(b) 幾何的來源：最短曲線與等周問題#
幾何這邊的問題同樣自然：
- Dido 問題（公元前 1 世紀，狄多諾）：給定兩端點，求連結它們的最短曲線。答案：直線。
- 等周問題（isoperimetric problem）：在周長固定的所有平面封閉曲線中，哪一條面積最大？答案：圓形。這是「帶約束的泛函極值」的原型。
- Dido 的「等周曲線」：對弧長 \(s\) 而言， \(J[y] = \int y\,\mathrm{d}s\) 是面積；「弧長固定」就是一個積分等式約束。
- 最小曲面（soap film）：Plateau 問題——給定邊界，撐開的肥皂膜是什麼？這導向「面積泛函」與 1880 年代 Weierstrass 的嚴格化。
(c) 一個共同的結構#
把上面所有問題寫成一個模板：
** \(J\) 的自變數是函數 \(y\) 本身，不是數。** 因此「求導」意味著：
其中 \(\eta\) 是任意零端點的光滑函數（稱為變分 \(\delta y\) ）。
一階展開（由微積分基本定理與分部積分）給出：
對 \(\eta'\) 分部積分（ \(\eta(a) = \eta(b) = 0\) ）：
所以泛函導數的「核心部分」是
這就是 Euler–Lagrange 方程的來源。
13.1.2 Euler 的力學（1752）與 Euler–Lagrange 方程#
Leonhard Euler 在 1752 年到 1760 年代處理「用積分求微分方程」的問題時，系統地得到了這條方程。今天我們把它寫成：
其中「 \(\mathrm{d}/\mathrm{d}x\) 」是全導數： \(\frac{\partial F}{\partial y'}\) 一般是 \(x, y, y'\) 的函數，所以要沿著曲線求導：
關鍵對照：這與第 3 章的導數定義是同一個「極限 + 差商」結構，只是被求導的對象從「數」升級為「函數」。
| 第 3 章 | 本節 | |
|---|---|---|
| 被求導的對象 | 數 \(x\) | 函數 \(y(x)\) |
| 增量 | \(h\) （實數） | \(\varepsilon\eta(x)\) （函數） |
| 極限 | \(\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}\) | \(\lim_{\varepsilon\to 0}\frac{J[y+\varepsilon\eta]-J[y]}{\varepsilon}\) |
| 結果 | 導數 \(f'(x)\) | Euler–Lagrange 方程 |
SymPy 怎麼處理「對函數求導」？
SymPy 沒有直接的 variation 函式，但它有兩件工具可以拼出來：
- 把 \(y\) 視為未定義函數
sp.Function('y')(x)，則sp.diff(F, y)就是 \(\frac{\partial F}{\partial y}\) （偏導數， \(y\) 被視為獨立座標），sp.diff(F, yp)就是 \(\frac{\partial F}{\partial y'}\) ； - 全導數要手動展開成 \(\frac{\partial G}{\partial x} + \frac{\partial G}{\partial y}y' + \frac{\partial G}{\partial y'}y''\) —— 或者，直接對含
Derivative(y(x), x)的表達式呼叫sp.diff(..., x)，SymPy 會自動做鏈鎖法則，得出同樣的結果。
下面這段程式碼把整條推導寫成可執行程序，並驗證兩個具體解（懸鏈線與擺線）。
import sympy as sp
x = sp.symbols('x', real=True)
y_s, yp, ypp = sp.symbols('y yp ypp', real=True)
def euler_lagrange(F):
"""回傳 EL 方程的左邊 L = F_y - d/dx(F_{y'})，展開成 y, yp, ypp 的形式。
F 依賴三個符號 x, y_s（= y(x)）, yp（= y'(x)）。
全導數 d/dx G = G_x + G_y * y' + G_{y'} * y''。
"""
Fy = sp.diff(F, y_s)
Fyp = sp.diff(F, yp)
dFyp_dx = (sp.diff(Fyp, x)
+ sp.diff(Fyp, y_s) * yp
+ sp.diff(Fyp, yp) * ypp)
return sp.simplify(Fy - dFyp_dx)
def verify(el, ysol):
"""把候選解 ysol(x) 及其一、二階導數代入 EL 方程，結果應為 0。"""
return sp.simplify(el.subs({y_s: ysol,
yp: sp.diff(ysol, x),
ypp: sp.diff(ysol, x, 2)}))
a, b = sp.symbols('a b', positive=True)
# --- 案例 1：最短曲線（Dido 問題），F = sqrt(1 + y'^2) ---
F1 = sp.sqrt(1 + yp**2)
el1 = euler_lagrange(F1)
print("案例 1 最短曲線：F =", F1)
print(" EL 方程左邊 =", el1)
m, c = sp.symbols('m c', real=True)
print(" 代回直線 y = m*x + c :", verify(el1, m * x + c))
# --- 案例 2：懸鏈線（均勻鏈條的重力位能），F = y*sqrt(1 + y'^2) ---
F2 = y_s * sp.sqrt(1 + yp**2)
el2 = euler_lagrange(F2)
print()
print("案例 2 懸鏈線：F =", F2)
print(" EL 方程左邊 =", el2)
print(" 代回懸鏈線 y = a*cosh((x-b)/a) :", verify(el2, a * sp.cosh((x - b) / a)))
# --- 案例 3：brachistochrone（最短下降時間），F = sqrt(1 + y'^2)/sqrt(y) ---
F3 = sp.sqrt(1 + yp**2) / sp.sqrt(y_s)
el3 = euler_lagrange(F3)
print()
print("案例 3 brachistochrone：F =", F3)
print(" EL 方程左邊 =", el3)
# Beltrami 恆等式：沿極值曲線 F - y'*F_{y'} = 常數（省掉 x 的依賴）
beltrami = sp.simplify(F3 - yp * sp.diff(F3, yp))
print(" Beltrami 量 F - y'F_{y'} =", beltrami)
# 擺線 x = a(t - sin t), y = a(1 - cos t)：以參數形式驗證 EL 殘差為 0
t = sp.symbols('t', positive=True)
xs_p = a * (t - sp.sin(t))
ys_p = a * (1 - sp.cos(t))
yp_p = sp.diff(ys_p, t) / sp.diff(xs_p, t)
ypp_p = sp.diff(yp_p, t) / sp.diff(xs_p, t)
print(" 沿擺線的 EL 殘差 =", sp.simplify(el3.subs({y_s: ys_p, yp: yp_p, ypp: ypp_p})))
實際輸出：
案例 1 最短曲線：F = sqrt(yp**2 + 1)
EL 方程左邊 = -ypp/(yp**2 + 1)**(3/2)
代回直線 y = m*x + c : 0
案例 2 懸鏈線：F = y*sqrt(yp**2 + 1)
EL 方程左邊 = (-y*ypp + yp**2 + 1)/(yp**2 + 1)**(3/2)
代回懸鏈線 y = a*cosh((x-b)/a) : 0
案例 3 brachistochrone：F = sqrt(yp**2 + 1)/sqrt(y)
EL 方程左邊 = (-2*y*ypp - yp**2 - 1)/(2*y**(3/2)*(yp**2 + 1)**(3/2))
Beltrami 量 F - y'F_{y'} = 1/(sqrt(y)*sqrt(yp**2 + 1))
沿擺線的 EL 殘差 = 0
逐行解讀：
- 案例 1 的 EL 方程是 \(\frac{y''}{(1+y'^2)^{3/2}} = 0\) ，即 \(y'' = 0\) ——弧線是直線。這個結論太簡單到讓人懷疑，但想想看：兩點之間的直線當然最短。真正的趣味在於下一個案例。
- 案例 2 的 EL 方程是 \(y\,y'' = 1 + y'^2\) ，這正是懸鏈線的微分方程。代回 \(y = a\cosh\frac{x-b}{a}\) 得 ** \(0\) ** ，證實了它是極值曲線。
- 案例 3 的 EL 方程是 \(2y\,y'' + y'^2 + 1 = 0\) ，即 brachistochrone 方程。代回擺線參數式，殘差為 ** \(0\) ** 。
- Beltrami 量：若 \(F\) 不顯含 \(x\) ，則沿極值曲線 \(F - y'F_{y'} = \text{常數}\) 。這是化簡二階 EL 方程為一階的捷徑，也是今天最佳化理論中「對偶性」的思想雛形。
注意這裡的「嚴格性」層次：我們用 SymPy 證明了「若 \(y\) 滿足 EL 方程，則 \(J\) 的泛函導數為零」，但沒有證明「滿足 EL 方程 ⟹ 一定是極值」（這需要二階變分與 Weierstrass 的充分條件，第 13.1.4 節會說明）。符號運算給出的是必要條件，不是充分條件。
13.1.3 Jacobi 的「最後乘法因子」（1845）#
Euler 的方法有個明顯的侷限：它直接處理 \(\int F(x,y,y')\,\mathrm{d}x\) 這種形式。但很多物理問題裡，作用量裡顯含時間：
這時 EL 方程變成
（外加一個 \(-\frac{\mathrm{d}L}{\mathrm{d}t}\) 的項，若 \(L\) 顯含 \(t\) ）。Jacobi 在 1845 年（正是在本節的標題年份）找到了把所有這些情形統一起來的方法。
Jacobi 的洞察是這樣的：
不要急著找 \(\delta y\) ，而是先找一個乘子 \(M(x, y, y')\) 使得
能被「湊成」 \(\delta(\text{某個表達式})\) 的形式，而這個表達式對 \(y\) 是線性的。
這個 \(M\) 就是他稱為的 「最後乘法因子」（the last multiplier） 。做法是：
- 從 \(\delta y\) 的任意變分出發，構造 \(\delta(\text{邊界項})\) ，使得
- 把後者再「積分掉」一次，得到
- 令括號中的 EL 表達式為零（因為 \(\delta y\) 任意）。
這個方法的力量在於它是「模」的（modular）：它說明 EL 方程是「作用量的某個函數 \(\Phi\) 的一階變分等於零」這個更一般原理的特例。
現代的對應：今天在最佳化理論裡，這個思想叫做 Lagrange 乘子法（第 13.3 節會用 SymPy 實作）與對偶性。Jacobi 1845 年的「最後乘法因子」，正是 Lagrange 1784 年「帶約束極值」的對偶版本。
順帶一提，Jacobi 對變分法的貢獻還包括 Jacobi 場（Jacobi fields）與 Jacobi 充分條件：他證明了二階變分的正定性可以保證某個解是嚴格極小值，而不只是駐值。這是「必要條件」升級為「充分條件」的第一步。
13.1.4 Legendre 與 Weierstrass 的嚴格化#
十九世紀中葉，Euler–Lagrange 方程被大量使用，但它的地位一直很含糊：它是必要條件？充分條件？兩者皆是？在不同的問題裡答案不同。
| 年份 | 人物 | 貢獻 |
|---|---|---|
| 1786 | Legendre | 引入二階變分 \(\delta^2 J\) ，並區分「極小」與「極大」；提出 Legendre 條件 \(\frac{\partial^2 F}{\partial y'^2} \ge 0\) |
| 1834 | Hamilton | 最小作用量原理的完整陳述，並把變分法用於力學 |
| 1845 | Jacobi | 「最後乘法因子」；Jacobi 充分條件（二階變分正定 ⟹ 嚴格極小） |
| 1850–1860 | Weierstrass | Weierstrass 必要條件：極值曲線必須滿足一個關於「速度」的嚴格不等式；並給出 Hilbert 過剩變分的準備 |
| 1900s | Hilbert | Hilbert 的第 20 問題：證明把變分問題正則化後，極值曲線是分段解析的（Kneser 1902 給出部分答案） |
Weierstrass 的必要條件尤其深刻。它說：若 \(y\) 是極值曲線，則對任意可微的變分 \(\eta\) 與任意 \(\varepsilon\) ，都必須滿足
其中 \(\varepsilon\) 是一個「速度」參數。這個不等式比 EL 方程強：EL 方程只是令 \(\varepsilon = 0\) 的特例。Weierstrass 條件的直觀意義是：在極值曲線的每一點，力的方向必須「朝向可行方向」，而不只是「平衡」。
Legendre 條件 \(\frac{\partial^2 F}{\partial y'^2} \ge 0\) 則是一個局部凸性判準：它保證在端點固定、 \(\eta\) 足夠小時， \(J\) 的變分是凸的。
這裡有一個重要的「計算機時代」註腳：所有這些條件，現在都可以用 SymPy 半自動檢查——因為它們全是偏導數與全導數的代數式。我們在 13.1.2 節的 verify() 函式正是這個思路的最小版本。
13.1.5 從變分法到梯度下降：泛函導數就是「無窮維梯度」#
現在來做本節最重要的概念轉換。回到第 7 章（多變量微積分）的梯度：
對有限維的 \(f : \mathbb{R}^n \to \mathbb{R}\) ，梯度 \(\nabla f\) 的定義是
對泛函 \(J : C[a,b] \to \mathbb{R}\) ，完全類比：
這是同一件事，只是「內積」從 \(\mathbb{R}^n\) 上的點積換成了 \(L^2\) 空間上的積分。
於是：
| 有限維最佳化 | 變分法 |
|---|---|
| \(f(\mathbf{x})\) | \(J[y]\) |
| 梯度 \(\nabla f\) | 泛函導數 \(\delta J / \delta y\) |
| 梯度下降 \(\mathbf{x} \leftarrow \mathbf{x} - \alpha\nabla f\) | 梯度下降 \(y \leftarrow y - \alpha\,\delta J/\delta y\) |
| Hessian 正定 ⟹ 唯一極小 | 二階變分正定 ⟹ 嚴格極小（Legendre–Jacobi 條件） |
| 拉格朗日乘數 | 帶約束的變分（Erdmann 角點條件） |
「梯度下降」在有限維與無窮維之間沒有本質差別。 差別只在於：
- 無窮維的「梯度下降步」是一個函數，不是向量；
- 要把無窮維問題變成可計算的，必須離散化——而離散化正是 13.3 節「符號最佳化」與現代計算機「數值最佳化」的會合處。
這裡我們做一個關鍵的示範：把泛函 \(J[y] = \int_0^1 \sqrt{1+y'^2}\,\mathrm{d}x\) 離散化成 \(N\) 個變數的有限維問題，然後直接用梯度下降求解，並與 SymPy 算出的 EL 方程解對照。
import numpy as np
import sympy as sp
# ============================================================
# (a) 把泛函離散化：J[y] = int_0^L sqrt(1 + y'^2) dx（最短曲線）
# ============================================================
L, N = 1.0, 21
xs = np.linspace(0.0, L, N)
h = L / (N - 1)
def J_discrete(y):
"""梯形／矩形離散化：J ≈ sum sqrt(1 + y'^2) * h。"""
d = np.diff(y) / h
return float(np.sum(np.sqrt(1.0 + d**2)) * h)
def grad_discrete(y):
"""離散泛函對內部節點的梯度。
dJ/dy_j = d_{j-1}/sqrt(1+d_{j-1}^2) - d_j/sqrt(1+d_j^2)
（端點固定，故只有內部節點有梯度。）
"""
d = np.diff(y) / h
g = np.zeros_like(y)
g[1:-1] = d[:-1] / np.sqrt(1 + d[:-1]**2) - d[1:] / np.sqrt(1 + d[1:]**2)
return g
# 先用中心差分「獨立」核對解析梯度
rng = np.random.default_rng(7)
yt = rng.normal(size=N) * 0.1
yt[0] = yt[-1] = 0.0
ga = grad_discrete(yt)
eps = 1e-6
fd = np.zeros_like(yt)
for j in range(1, N - 1):
yp_ = yt.copy(); yp_[j] += eps
ym_ = yt.copy(); ym_[j] -= eps
fd[j] = (J_discrete(yp_) - J_discrete(ym_)) / (2 * eps)
print("解析梯度與中心差分梯度最大偏差 =", np.max(np.abs(ga[1:-1] - fd[1:-1])))
# 帶回溯線搜索的梯度下降
y = 0.3 * np.sin(3 * np.pi * xs) + 0.2
y[0] = y[-1] = 0.0
print("初始 J =", J_discrete(y))
for it in range(300):
g = grad_discrete(y)
step = 0.5
for _ in range(40): # 回溯線搜索
cand = y - step * g
cand[0], cand[-1] = 0.0, 0.0
if J_discrete(cand) < J_discrete(y):
break
step *= 0.5
y = cand
print("梯度下降後 J =", J_discrete(y))
print("最大偏差（理論解為直線 y=0）=", np.max(np.abs(y)))
# ============================================================
# (b) 符號對照：EL 方程
# ============================================================
x = sp.symbols('x', real=True)
y_s, yp, ypp = sp.symbols('y yp ypp', real=True)
F = sp.sqrt(1 + yp**2)
Fy = sp.diff(F, y_s)
Fyp = sp.diff(F, yp)
dFyp_dx = sp.diff(Fyp, x) + sp.diff(Fyp, y_s) * yp + sp.diff(Fyp, yp) * ypp
el = sp.simplify(Fy - dFyp_dx)
print()
print("F =", F, " EL 方程左邊 =", el)
m, c = sp.symbols('m c', real=True)
print("直線 y = m*x + c 代回 =", sp.simplify(el.subs({y_s: m * x + c, yp: m, ypp: 0})))
實際輸出：
解析梯度與中心差分梯度最大偏差 = 2.211770766535892e-10
初始 J = 2.501493052381671
梯度下降後 J = 1.0000028458485466
最大偏差（理論解為直線 y=0）= 0.0010694722222667084
F = sqrt(yp**2 + 1) EL 方程左邊 = -ypp/(yp**2 + 1)**(3/2)
直線 y = m*x + c 代回 = 0
逐行解讀：
- 第一行是可信度測試：解析梯度與獨立算出來的中心差分梯度吻合到 \(2\times10^{-10}\) 。這正是第 1.3.6 節「把猜升級成證」的精神：符號推導給出的公式，必須用一個獨立路徑核對。
- **梯度下降後 \(J \approx 1.0000028\) ** ：最短曲線的理論值是 \(\sqrt{1 + 0^2}\cdot 1 = 1\) （端點同高的直線長度），吻合到 \(10^{-6}\) 。
- 符號端：EL 方程是 \(y'' = 0\) ，解是直線，代回得 \(0\) 。數值端收斂到直線，符號端說直線是解——兩條獨立的路給出同一個答案。
把嚴格化條件也算出來：Legendre 與 Weierstrass#
13.1.4 節提到的那些「嚴格化條件」，在 SymPy 裡全是偏導數的代數式。我們把它們逐條算出來：
import sympy as sp
x = sp.symbols('x', real=True)
y_s, yp, ypp = sp.symbols('y yp ypp', real=True)
def euler_lagrange(F):
"""EL 方程左邊 L = F_y - d/dx(F_{y'})。"""
Fy = sp.diff(F, y_s)
Fyp = sp.diff(F, yp)
dFyp_dx = (sp.diff(Fyp, x)
+ sp.diff(Fyp, y_s) * yp
+ sp.diff(Fyp, yp) * ypp)
return sp.simplify(Fy - dFyp_dx)
def legendre_condition(F):
"""Legendre 必要條件：d2F/dy'^2 >= 0（局部凸性）。"""
return sp.simplify(sp.diff(F, yp, 2))
def weierstrass_condition(F, v):
"""Weierstrass 必要條件：E := F - y'F_{y'} - v F_{y'} 對 v 單調不減。"""
Fyp = sp.diff(F, yp)
E = sp.simplify(F - yp * Fyp - v * Fyp)
return E, sp.simplify(sp.diff(E, v))
cases = [
("最短曲線 F = sqrt(1+y'^2)", sp.sqrt(1 + yp**2)),
("懸鏈線 F = y*sqrt(1+y'^2)", y_s * sp.sqrt(1 + yp**2)),
("brachistochrone F = sqrt(1+y'^2)/sqrt(y)", sp.sqrt(1 + yp**2) / sp.sqrt(y_s)),
("退化案例 F = y'^4", yp**4),
]
for name, F in cases:
print(name)
print(" EL 方程左邊 =", euler_lagrange(F))
print(" d2F/dy'^2 (Legendre) =", legendre_condition(F))
E, dE = weierstrass_condition(F, sp.Symbol('v'))
print(" Weierstrass E =", E)
print(" dE/dv =", dE)
print()
# Jacobi 充分條件：二階變分 delta^2 J[eta][eta] 需正定
F = y_s * sp.sqrt(1 + yp**2)
eta = sp.Function('eta')(x)
etap = sp.diff(eta, x)
d2Fy = sp.simplify(sp.diff(F, y_s, 2))
d2Fyyp = sp.simplify(sp.diff(sp.diff(F, y_s), yp))
d2Fyp = sp.simplify(sp.diff(F, yp, 2))
d2 = sp.simplify(d2Fy * eta**2 + 2 * d2Fyyp * eta * etap + d2Fyp * etap**2)
print("懸鏈線的二階變分 =", d2)
print("其中 d2F/dy'^2 =", d2Fyp)
實際輸出：
最短曲線 F = sqrt(1+y'^2)
EL 方程左邊 = -ypp/(yp**2 + 1)**(3/2)
d2F/dy'^2 (Legendre) = (yp**2 + 1)**(-3/2)
Weierstrass E = (-v*yp + 1)/sqrt(yp**2 + 1)
dE/dv = -yp/sqrt(yp**2 + 1)
懸鏈線 F = y*sqrt(1+y'^2)
EL 方程左邊 = (-y*ypp + yp**2 + 1)/(yp**2 + 1)**(3/2)
d2F/dy'^2 (Legendre) = y/(yp**2 + 1)**(3/2)
Weierstrass E = y*(-v*yp + 1)/sqrt(yp**2 + 1)
dE/dv = -y*yp/sqrt(yp**2 + 1)
brachistochrone F = sqrt(1+y'^2)/sqrt(y)
EL 方程左邊 = (-2*y*ypp - yp**2 - 1)/(2*y**(3/2)*(yp**2 + 1)**(3/2))
d2F/dy'^2 (Legendre) = 1/(sqrt(y)*(yp**2 + 1)**(3/2))
Weierstrass E = (-v*yp + 1)/(sqrt(y)*sqrt(yp**2 + 1))
dE/dv = -yp/(sqrt(y)*sqrt(yp**2 + 1))
退化案例 F = y'^4
EL 方程左邊 = -12*yp**2*ypp
d2F/dy'^2 (Legendre) = 12*yp**2
Weierstrass E = yp**3*(-4*v - 3*yp)
dE/dv = -4*yp**3
懸鏈線的二階變分 = (y*Derivative(eta(x), x) + 2*yp*(yp**2 + 1)*eta(x))*Derivative(eta(x), x)/(yp**2 + 1)**(3/2)
其中 d2F/dy'^2 = y/(yp**2 + 1)**(3/2)
幾點解讀：
- 「最短曲線」與「brachistochrone」的 Legendre 條件恆正： \((1+y'^2)^{-3/2} > 0\) 與 \(\frac{1}{\sqrt{y}(1+y'^2)^{3/2}} > 0\) 。這正是為什麼這兩個問題「好解」 —— 泛函是凸的，EL 方程的解就是唯一的全域最小值點。「凸性」這個詞在第 7 章（多元函數的 Hessian）出現過，這裡是它的無窮維版本。
- **「懸鏈線」的 Legendre 條件是 \(\frac{y}{(1+y'^2)^{3/2}}\) ** ，符號取決於 \(y\) 。在鏈條下垂的區域 \(y > 0\) ，條件成立（局部凸）。
- 退化案例 \(F = y'^4\) ** 是最有教育意義的： \(\frac{\partial^2F}{\partial y'^2} = 12y'^2 \ge 0\) 居然成立，但 \(F\) 在 \(y' \to \infty\) 時無界，泛函沒有最小值。這說明「Legendre 條件成立」既不充分、也不是全局性的判準。**
- Weierstrass 條件 \(\frac{\partial E}{\partial v} \ge 0\) 在這些例子中並不總成立（例如最短曲線的 \(\frac{\partial E}{\partial v} = -\frac{y'}{\sqrt{1+y'^2}}\) ，只有在 \(y' \le 0\) 時成立）——這正好說明：Weierstrass 條件是對「特定方向的變分」而言的，不是無條件成立。
把這些條件與第 7 章對照會很清楚：
| 有限維 | 變分法 |
|---|---|
| Hessian \(H\) 正定 | Legendre 條件 \(\frac{\partial^2F}{\partial y'^2} \ge 0\) + Jacobi 條件 |
| 目標函數凸 | \(F\) 關於 \((y,y')\) 凸 |
| 對偶性（KKT 條件） | Jacobi 的「最後乘法因子」、帶約束的 EL 方程 |
SymPy 的角色在這裡非常精確：它不做凸性分析、不做全域最佳化保證，它只做代數。把代數結論升級成定理，仍然是人的工作。
13.1.6 從變分法到機器學習：Robbins–Monro（1951）#
1951 年，Herbert Robbins 與 Sutton Monro 發表了一篇短文，處理的是隨機逼近問題：
給定 \(E[Y_n] = f(x)\) 的隨機觀測 \(Y_n\) ，求 \(x\) 使得 \(f(x) = 0\) 。
他們的定理（Robbins–Monro 條件）說，若步長 \(\alpha_n\) 滿足
則迭代 \(x_{n+1} = x_n - \alpha_n (f(x_n) + \text{噪聲})\) 以機率 1 收斂到 \(f\) 的根。
** \(t^{-\alpha}\) 型步長**（ \(0 < \alpha \le 1\) ）正好滿足這兩個條件： \(\sum t^{-\alpha} = \infty\) （走得夠久）而 \(\sum t^{-2\alpha} < \infty\) （抖動會被抹平）。
今天深度學習裡的 SGD 用的是同一個定理，只是把 \(f\) 換成損失函數的梯度、把噪聲換成 mini-batch 抽樣的方差。機器學習的「學習律」有一半的數學根基在 1951 年就埋好了。
回到變分法的語言：SGD 就是用隨機估計的泛函導數做梯度下降：
我們用一個線性回歸做示範，把「批量梯度下降」與「SGD」並排比較。
import numpy as np
# ============================================================
# Robbins-Monro（1951）：隨機逼近與 SGD 的步長條件
# ============================================================
rng = np.random.default_rng(2026)
n = 200
X = rng.uniform(-2, 2, n)
w_true = np.array([1.5, -0.7])
y_obs = w_true[0] * X + w_true[1] + rng.normal(0, 0.3, n)
all_idx = np.arange(n)
def run_gd(lr, steps=400):
"""批量梯度下降：用全部資料算梯度。"""
w = np.array([0.0, 0.0])
for _ in range(steps):
r = y_obs - (w[0] * X + w[1])
g = np.array([-2 * np.mean(r * X), -2 * np.mean(r)])
w = w - lr * g
return w
def run_sgd(lr0, power, steps=20000, batch=8, seed=1):
"""SGD：每次只用 batch 個樣本估計梯度，步長 lr0 / t**power。"""
r_ = np.random.default_rng(seed)
w = np.array([0.0, 0.0])
for t in range(1, steps + 1):
idx = r_.integers(0, n, size=batch)
r = y_obs[idx] - (w[0] * X[idx] + w[1])
g = np.array([-2 * np.mean(r * X[idx]), -2 * np.mean(r)])
w = w - (lr0 / t**power) * g
return w
w_gd = run_gd(0.05)
print("批量梯度下降 w =", np.round(w_gd, 6))
for p in [0.5, 0.6, 0.75]:
w_sgd = run_sgd(0.5, p)
print(f"SGD（步長 0.5/t^{p}） w =", np.round(w_sgd, 6))
print("真實 w =", w_true)
# 步長太大會發散：Robbins-Monro 條件不是裝飾
w_bad = run_sgd(0.8, 0.0, steps=2000)
print("步長過大且不遞減的 SGD w =", np.round(w_bad, 3))
實際輸出：
批量梯度下降 w = [ 1.468258 -0.694645]
SGD（步長 0.5/t^0.5） w = [ 1.468541 -0.695211]
SGD（步長 0.5/t^0.6） w = [ 1.46853 -0.694102]
SGD（步長 0.5/t^0.75） w = [ 1.469363 -0.69433 ]
真實 w = [ 1.5 -0.7]
步長過大且不遞減的 SGD w = [28.074 4.638]
解讀：
- 批量 GD 與三種 SGD 都收斂到接近真值的位置。注意 SGD 的解比 GD 略「散」 ——因為 SGD 永遠停在「噪聲球」裡，這個現象在機器學習裡叫 noise floor。
- 步長 \(0.5/t^{0.6}\) 收斂得最好（偏差最小），這不是巧合：指數 \(p\) 要在 \(1/2\) 與 \(1\) 之間，太小則噪聲壓不下去，太大則收斂太慢——這正是 Robbins–Monro 兩個條件的權衡。
- 最後一行：步長固定為 \(0.8\) 且不遞減時，SGD 發散（ \(w\) 跑到 28 附近）。這就是為什麼今天所有的深度學習框架都強制要求學習率要衰減。
泛函導數的「隨機估計」長什麼樣？#
回到變分法：在機器學習裡，「泛函」是損失函數
這是有限和，它的「梯度」是可以精確算的：
SGD 用 mini-batch 估計這個總和。這裡沒有真正的「無窮維泛函導數」了——因為神經網路的參數雖然多，但仍是有限維。所以現代機器學習嚴格來說不是變分法，而是「用變分法的直覺處理的有限維最佳化」。
但這兩者之間有一個真實的、值得學習的張力：
- 變分法求的是整個函數 \(y(x)\) ，一次性給出解析形式（懸鏈線、擺線）；
- 機器學習求的是有限維但天文數字般大的參數向量，只能迭代逼近，得到的往往是「函數的近似」而非閉式。
SymPy 在這兩者之間的位置非常獨特：它是少數能直接處理符號泛函（ \(F(x, y, y')\) 的符號表達式）的工具。第 13.3 節我們會看到，SymPy 甚至提供了 sp.minimum / sp.maximum 這類直接處理「函數極值」的指令。
13.1.7 一個完整的示範：符號 × 數值 × 圖#
我們做一次三路合流：符號算出 EL 方程與精確解 → 數值用離散化梯度下降逼近 → 畫在一起比較。
這次挑一個「行為良好」的泛函：
為什麼挑它？因為它是嚴格凸的： \(F = y'^2 + y^2\) 關於 \((y, y')\) 的 Hessian 是 \(2I\) ，恆正定。唯一極小值點存在、全域唯一、數值上非常穩定。 這讓我們可以把注意力放在「符號與數值是否一致」上，而不用擔心最佳化本身的困難。
順帶一個重要的教訓：懸鏈線問題 \(J[y] = \int y\sqrt{1+y'^2}\,\mathrm{d}x\) 若不加上「鏈長固定」的約束，其泛函是無下界的——把曲線往下拉得越深，位能越負，沒有最小值。物理上真正的懸鏈線問題必須加上「總長度 \(L\) 固定」這個等周約束， \(L\) 就是 Lagrange 乘子要對付的對象。這是「符號代數給出必要條件，但充分性需要人來把關」的又一例證，也是第 13.3 節 Lagrange 乘數法的自然入口。
import numpy as np
import sympy as sp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# ============================================================
# (a) 符號端：EL 方程、精確解、Hessian
# ============================================================
x = sp.symbols('x', real=True)
y_s, yp, ypp = sp.symbols('y yp ypp', real=True)
F = yp**2 + y_s**2
Fy = sp.diff(F, y_s)
Fyp = sp.diff(F, yp)
dFyp_dx = (sp.diff(Fyp, x)
+ sp.diff(Fyp, y_s) * yp
+ sp.diff(Fyp, yp) * ypp)
el = sp.simplify(Fy - dFyp_dx)
print("F =", F)
print("Euler-Lagrange 方程左邊 =", el)
exact_expr = sp.sinh(x) / sp.sinh(1)
print("代回精確解 y = sinh(x)/sinh(1) :",
sp.simplify(el.subs({y_s: exact_expr, yp: sp.diff(exact_expr, x),
ypp: sp.diff(exact_expr, x, 2)})))
Hm = sp.Matrix([[sp.diff(F, y_s, 2), sp.diff(sp.diff(F, y_s), yp)],
[sp.diff(sp.diff(F, yp), y_s), sp.diff(F, yp, 2)]])
print("Hessian d2F/d(y,y') =", Hm.tolist(), " 特徵值 =", Hm.eigenvals())
# ============================================================
# (b) 數值端：離散化 + 梯度下降
# ============================================================
M = 101
xg = np.linspace(0.0, 1.0, M)
dx = xg[1] - xg[0]
def J_disc(yy):
"""梯形法離散化 J = int (y'^2 + y^2) dx。"""
d = np.diff(yy) / dx
return float(np.sum((d**2 + 0.5 * (yy[:-1]**2 + yy[1:]**2)) * dx))
def grad_J(yy):
"""離散泛函的解析梯度：g_j = 2(d_{j-1} - d_j) + 2 y_j dx。"""
d = np.diff(yy) / dx
g = np.zeros_like(yy)
g[1:-1] = 2 * (d[:-1] - d[1:]) + 2 * yy[1:-1] * dx
return g
# 先用中心差分「獨立」核對解析梯度
rng = np.random.default_rng(11)
yt = rng.normal(size=M) * 0.1
yt[0], yt[-1] = 0.0, 1.0
ga = grad_J(yt)
fd = np.zeros(M)
eps = 1e-6
for j in range(1, M - 1):
p = yt.copy(); p[j] += eps
m = yt.copy(); m[j] -= eps
fd[j] = (J_disc(p) - J_disc(m)) / (2 * eps)
print()
print("解析梯度 vs 中心差分：最大偏差 =", np.max(np.abs(ga[1:-1] - fd[1:-1])))
yv = np.linspace(0.0, 1.0, M) + 0.2 * np.sin(2 * np.pi * xg)
yv[0], yv[-1] = 0.0, 1.0
print("初始 J =", J_disc(yv))
for it in range(4000):
gr = grad_J(yv)
step = 0.5
for _ in range(50):
cand = yv - step * gr
cand[0], cand[-1] = 0.0, 1.0
if J_disc(cand) < J_disc(yv):
break
step *= 0.5
yv = cand
print("梯度下降後 J =", J_disc(yv))
exact = np.array([float(exact_expr.subs(x, xv)) for xv in xg])
print("最大偏差 =", float(np.max(np.abs(yv - exact))))
print("精確解的 J =", float(sp.integrate(
F.subs({y_s: exact_expr, yp: sp.diff(exact_expr, x)}), (x, 0, 1))))
# ============================================================
# (c) 畫圖
# ============================================================
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(xg, yv, "o-", ms=3, label="gradient descent (discretized)")
ax.plot(xg, exact, "r-", lw=2, label="exact sinh(x)/sinh(1)")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Variational problem: discrete gradient descent vs exact solution")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("/tmp/ch13_var.png", dpi=120)
print("saved /tmp/ch13_var.png")
實際輸出：
F = y**2 + yp**2
Euler-Lagrange 方程左邊 = 2*y - 2*ypp
代回精確解 y = sinh(x)/sinh(1) : 0
Hessian d2F/d(y,y') = [[2, 0], [0, 2]] 特徵值 = {2: 2}
解析梯度 vs 中心差分：最大偏差 = 3.247422730723315e-08
初始 J = 2.0790175967838276
梯度下降後 J = 1.3131686360633823
最大偏差 = 0.004577196301303155
精確解的 J = 1.3130352854993312
saved /tmp/ch13_var.png
三路合流：
| 路徑 | 給出的答案 |
|---|---|
| 符號（EL 方程 + 代回驗證） | \(y'' = y\) ，代回精確解殘差 \(=0\) |
| 符號（Hessian） | \(2I\) ，特徵值 \(2, 2 > 0\) ：嚴格凸，唯一極小 |
| 符號（直接積分） | \(J[\text{精確解}] = 1.3130352854993312\) |
| 數值（離散化 + 梯度下降） | \(J = 1.3131686360633823\) ，曲線最大偏差 \(0.0046\) |
注意最後兩行的關係： \(J\) 的誤差是 \(1.3\times10^{-4}\) ，而曲線的最大偏差是 \(4.6\times10^{-3}\) ，大了約 35 倍。這是本節最實用的一個觀察：
泛函的值比泛函的極值點更準。
原因很直觀： \(J\) 是二次型， \(J[y] - J[y_*] = \int \left((\delta y)'^2 + (\delta y)^2\right)\mathrm{d}x\) ，對偏差 \(\delta y\) 是平方關係；而曲線本身的位置偏差是一階的。這是第 13.2 節誤差理論的核心直覺之一，在最佳化裡就已經出現。
** \(0.0046\) 的偏差是什麼量級？** 此處 \(h = 1/100\) ，而偏差 \(\approx 0.0046 \approx 0.46h\) ——** \(O(h)\) ** ，來自離散化本身。這提醒我們：符號與數值不是對立的，而是互補的：符號給出結構與精確形式，數值給出在符號難以處理時的可行解，並且用數值獨立驗證符號的結論。
13.1.8 一張圖串起本節#
flowchart TD
A["Dido 問題 / 等周問題<br/>公元前 1 世紀<br/>最短曲線、等周"] --> B["Hamilton 最小作用量原理 1834<br/>S = int L dt 在真實軌跡上取駐值"]
C["Euler 1752<br/>由牛頓力學導出<br/>F_y - d/dx F_y' = 0"] --> D["Euler-Lagrange 方程<br/>泛函的必要條件"]
B --> D
D --> E["Jacobi 1845<br/>最後乘法因子<br/>統一顯含時間的情形"]
E --> F["Legendre 1786 / Weierstrass 1850s<br/>二階變分、必要與充分條件"]
D --> G["泛函導數 = 無窮維梯度<br/>delta J[y](eta)"]
G --> H["離散化：泛函 -> 有限維問題<br/>J_discrete(y_1,...,y_N)"]
H --> I["梯度下降<br/>y <- y - alpha * grad J"]
I --> J["Robbins-Monro 1951<br/>隨機梯度下降 SGD<br/>機器學習的引擎"]
D --> K["SymPy: sp.diff + 手動全導數<br/>推導並驗證 EL 方程"]
K --> L["SymPy: sp.solve / dsolve<br/>解 EL 方程得精確解"]
I --> M["SymPy + numpy<br/>離散化梯度下降核對"]習題#
習題 13.1.1（推導 Euler–Lagrange 方程）
設 \(J[y] = \int_0^{\pi/2} \left(y'(x)^2 + y(x)^2\right)\,\mathrm{d}x\) ，端點 \(y(0) = 0\) 、 \(y(\pi/2) = 1\) 。
(a) 用 13.1.2 節的 euler_lagrange() 模板推出 EL 方程； (b) 用 sp.dsolve 解出 \(y(x)\) ； (c) 用 verify() 驗證你的解代回 EL 方程殘差為 \(0\) ； (d) 說明為什麼這個 EL 方程是線性的，並指出「線性泛函的 EL 方程是線性的」的一般原因。
習題 13.1.2（Beltrami 恆等式與降階）
對 brachistochrone 的 \(F = \frac{\sqrt{1+y'^2}}{\sqrt{y}}\) ：
(a) 證明 Beltrami 量 \(F - y'F_{y'}\) 沿極值曲線為常數（用 \(F\) 不顯含 \(x\) 這件事）； (b) 由 Beltrami 恆等式導出一階方程 \(y' = \sqrt{\frac{1}{c^2 y} - 1}\) （ \(c\) 為常數）； (c) 用 sp.integrate 或 sp.dsolve 由此導出擺線的參數式，並與 13.1.2 節的符號結果對照。
習題 13.1.3（等周問題與 Lagrange 乘數）
「周長固定為 \(L\) 的封閉曲線中，面積最大的是圓」。
(a) 把面積寫成泛函 \(A[y] = \int y\,\mathrm{d}s = \int y\sqrt{1+y'^2}\,\mathrm{d}x\) ； (b) 周長約束寫成 \(\int \sqrt{1+y'^2}\,\mathrm{d}x = L\) ； (c) 用帶約束的變分（Jacobi 的「最後乘法因子」就是這個工具）引入乘子 \(\lambda\) ，推導出帶 \(\lambda\) 的 EL 方程； (d) 證明 \(\lambda = 1\) 時解是圓（ \(y(x) = \sqrt{R^2 - (x-x_0)^2}\) 的某一段）。並說明：為什麼 Jacobi 的方法比 Lagrange 1784 年的方法更一般？
習題 13.1.4（Legendre 必要條件）
對 \(J[y] = \int_0^1 F(x,y,y')\,\mathrm{d}x\) ， \(F = y'^{\,4}\) 與 \(F = \sqrt{1+y'^2}\) 兩個案例：
(a) 計算 \(\frac{\partial^2 F}{\partial y'^2}\) ； (b) 說明為什麼 \(F = y'^4\) 的 EL 方程 \(4y'^3 = 0\) 只有一個解（ \(y'=0\) ），而這正是 Legendre 條件 \(\frac{\partial^2F}{\partial y'^2} \ge 0\) 不成立（因為 \(12y'^2\) 在 \(y'\ne0\) 時為正，但函數在 \(\infty\) 處發散）的反映； (c) 討論： \(F = y'^4\) 的泛函是否有最小值？若沒有，這說明「EL 方程給出的是必要條件」這句話的什麼重要性？
習題 13.1.5（Robbins–Monro 的步長條件）
(a) 對 \(\alpha_n = n^{-\gamma}\) （ \(0 < \gamma \le 1\) ），證明 \(\sum \alpha_n = \infty\) 且 \(\sum \alpha_n^2 < \infty\) ； (b) 若 \(\gamma > 1\) ，說明為什麼 SGD 會收斂但停不到最優點； (c) 用 13.1.6 節的程式碼，掃描 \(\gamma \in \{0.3, 0.5, 0.7, 0.9, 1.0, 1.2\}\) ，記錄最終 \(w\) 與真值的距離，畫出曲線並說明你的觀察。
習題 13.1.6（泛函導數的離散對應）
對 \(J[y] = \int_0^1 y(x)^2 + y'(x)^2\,\mathrm{d}x\) ：
(a) 推導 EL 方程並求出解（端點 \(y(0)=y(1)=0\) ）； (b) 用 \(N = 21, 41, 81\) 離散化，用梯度下降求數值解，與 (a) 的精確解比較，觀察誤差量級； (c) 說明：為什麼這個例子中「 \(J\) 的誤差」與「曲線的誤差」量級不同？（提示： \(J\) 是二次型。）
小結#
- 變分法的問題來源：Dido 問題與等周問題（幾何）、Hamilton 的最小作用量原理（物理）。它們共同的形式是 \(J[y] = \int_a^b F(x,y,y')\,\mathrm{d}x\) ，其中被最佳化的對象是函數。
- Euler–Lagrange 方程（Euler 1752）： \(\frac{\partial F}{\partial y} - \frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial F}{\partial y'} = 0\) ，其中全導數要展開成 \(\frac{\partial G}{\partial x} + \frac{\partial G}{\partial y}y' + \frac{\partial G}{\partial y'}y''\) 。它是必要條件，不是充分條件。
- Jacobi 的「最後乘法因子」（1845） 用乘子把「顯含時間」等所有情形統一，並導出 Jacobi 充分條件。
- 嚴格化三部曲：Legendre 1786（二階變分）→ Jacobi 1845（充分條件）→ Weierstrass 1850s（必要不等式）；Legendre–Jacobi 條件在今天就是「Hessian 正定」。
- 泛函導數就是無窮維的梯度： \(\delta J[y](\eta) = \int \left(F_y - \frac{\mathrm{d}}{\mathrm{d}x}F_{y'}\right)\eta\,\mathrm{d}x\) ，對應梯度下降 \(y \leftarrow y - \alpha\,\delta J\) 。
- 離散化是橋樑：把泛函變成 \(N\) 個變數的有限維問題後，梯度下降、線搜索、牛頓法全部可用。符號（EL 方程）與數值（梯度下降）會收斂到同一個答案，偏差是 \(O(h)\) 。
- Robbins–Monro（1951） 給出 SGD 的收斂定理，條件是 \(\sum \alpha_n = \infty\) 且 \(\sum \alpha_n^2 < \infty\) ；這是今天深度學習裡學習率必須衰減的理論依據。
下一節我們離開「最佳化」這條線，轉向誤差：符號運算給出的是精確答案，但在真實計算中我們面對的是浮點數、截斷誤差、捨入誤差與病態問題——計算機時代的微積分是一門關於「答案有多可靠」的學問。