1983 - Ruth 辛積分器——保結構的數值方法#

案件摘要#

1983 年，Ruth 在 IEEE 論文中提出辛積分器（symplectic integrator）的第一個顯式三階方法：設計數值積分時，不只追求「單步精度」，而是要求每一步都保持辛結構（正則 2-形式）——

\[ \text{數值流 } \Phi_{\Delta t} \text{ 是辛同胚} \;\Longleftrightarrow\; \Phi_{\Delta t}^T J \Phi_{\Delta t} = J \]

傳統方法（如 RK4）單步誤差小，但長期模擬能量系統性漂移；辛積分器能量有微小振盪但永不漂移。天體力學、分子動力學、加速器物理的長期模擬從此改觀。

前因#

混沌與長期模擬的需求：Poincaré（1890）與 KAM（1954）揭示了相空間結構（見 1890-Poincare三體問題.md、1954-KAM定理.md），但傳統數值方法會破壞這些結構——KAM 環面在 RK4 的長期模擬中會「漏氣」。
劉維爾定理（1838）：正則流不可壓縮（見 1843-Jacobi正則變換.md）——數值流若可壓縮，長期統計必然失真。
FPUT 實驗（1955）的教訓：電腦模擬是第三種偵探方法（見 1955-FPUT實驗.md），但方法本身必須保真。
太陽系數值積分的困境：1970s 的長期模擬（百萬年尺度）顯示能量緩慢漂移——是物理還是數值誤差？無法分辨。

線索與推理#

偵查線索 1：傳統方法的漂移病因#

以諧振子 \(H = \frac{1}{2}(p^2 + q^2)\) 為例，Euler 顯式法：

\[ q_{n+1} = q_n + \Delta t\, p_n, \qquad p_{n+1} = p_n - \Delta t\, q_n \]

一步變換的矩陣 \(M = \begin{pmatrix} 1 & \Delta t \\ -\Delta t & 1 \end{pmatrix}\) ， \(\det M = 1 + \Delta t^2 > 1\) —— 面積放大，且 \(|M| > 1\) ⇒ 能量 \(E_n \sim (1 + \Delta t^2)^{n} E_0\) 指數漂移。病因：Euler 法不是辛的。

偵查線索 2：保辛的「手術」——半隱式 Euler#

只調整一步順序（用新算的 \(q\) 更新 \(p\) ）：

\[ q_{n+1} = q_n + \Delta t\, p_n, \qquad p_{n+1} = p_n - \Delta t\, q_{n+1} \]

矩陣 \(M = \begin{pmatrix} 1 & \Delta t \\ -\Delta t & 1-\Delta t^2 \end{pmatrix}\) ， \(\det M = 1\) —— 面積守恆（辛）！且可證明此方法精確求解另一個哈密頓系統

\[ \tilde H = H + \Delta t\, H_1 + \Delta t^2 H_2 + \cdots \quad(\text{後向誤差分析}) \]

辛積分器 = 精確求解一個「鄰近」的哈密頓系統——能量 \(\tilde H\) 精確守恆，故 \(H\) 只振盪不漂移。

偵查線索 3：Ruth 的三階構造與分裂法#

Ruth 把可分離哈密頓 \(H = T(p) + V(q)\) 的每一步拆成「踢–漂–踢」的組合：

\[ \Phi_{\Delta t} = e^{\Delta t\, c_1 A_V}\, e^{\Delta t\, d_1 A_T}\, e^{\Delta t\, c_2 A_V}\, \cdots \]

（ \(A_T, A_V\) 分別是 \(T\) 、 \(V\) 生成的流，各自精確可解。）用 BCH 展開（Baker–Campbell–Hausdorff）匹配 \(\dot z = J\nabla H\) 到三階，解出係數 \(c_i, d_i\) 。 分裂法 + BCH 展開 = 任意階辛積分器的構造機器。

對照表：數值方法的長期行為#

方法	單步精度	長期能量	相空間結構
顯式 Euler	一階	指數漂移	破壞
RK4	四階	緩慢漂移	破壞
半隱式 Euler（辛）	一階	振盪不漂移	保留（微形變）
Ruth 三階（辛）	三階	振盪不漂移	保留

Python 驗證：Euler vs 半隱式 Euler 的長期能量#

import numpy as np

dt, N = 0.05, 100000
q = p = 1.0
E_euler = [0.5*(q*q + p*p)]
qe, pe = q, p
for _ in range(N):                       # 顯式 Euler（不辛）
    qe, pe = qe + dt*pe, pe - dt*qe
    E_euler.append(0.5*(qe*qe + pe*pe))
E_simpl = [0.5*(q*q + p*p)]
qs, ps = q, p
for _ in range(N):                       # 半隱式 Euler（辛）
    qs = qs + dt*ps
    ps = ps - dt*qs
    E_simpl.append(0.5*(qs*qs + ps*ps))
print(f"Euler 漂移: {E_euler[-1]:.3e}（爆炸）")     # >> 1e10
print(f"辛法振盪: {min(E_simpl):.4f} ~ {max(E_simpl):.4f}")  # 恆定在 1.0 附近

結案#

「數值精度 ≠ 數值保真」的教訓確立：長期模擬必須保辛。
1988 年 Forest–Ruth 四階方法普及，太陽系億年級模擬成為可能（見 1988-幾何數值積分.md）。
連鎖反應：分子動力學（Verlet 算法 1967 年其實就是辛的，遲至 1980s 才被認出）、加速器束流物理、數值相對論——幾何數值積分成為獨立學科。