跳轉到

10.1 數值方法：Euler 與 Runge–Kutta#

本章導言#

絕大多數微分方程解不出顯式（Liouville 1841、2.3 節）；解析方法（級數、定性理論）給結構與行為，但工程需要的是數字：軌道的數值、預測的日期、控制的指令。數值方法用離散步進代替連續流動：

\[ y_{n+1} = y_n + h\, \Phi(y_n, h). \]

本節講述兩大經典——Euler 法（1768，最樸素）與 Runge–Kutta 法（1895–1901，最常用），以及它們的核心理論：局部截斷誤差、全局收斂階、穩定性與剛性。

Euler 法：直覺的原點#

Leonhard Euler 1768 年《積分法》 中提出：沿斜率場走折線。

\[ y_{n+1} = y_n + h\, f(t_n, y_n). \]

幾何：在 \((t_n, y_n)\) 處讀出斜率、直行一小步。誤差的來源：斜率在步內變化（用舊斜率走了新路）。

收斂理論（Grönwall 3.2 節的數值版）：

  • 局部截斷誤差：一步的理論誤差 \(\tau_n = O(h^2)\) （Taylor 展開：真解 \(y(t + h) = y + hf + \frac{h^2}{2}y'' + \cdots\) ，Euler 丟掉了 \(\frac{h^2}{2}y''\) 項）。
  • 全局誤差： \(N = \frac{T}{h}\) 步、每步誤差 \(O(h^2)\) 、且誤差會傳播（有界導數下每步傳播因子 \(\le 1 + Lh\) ）：總計
\[ \max_n |y_n - y(t_n)| \le C\, h \qquad \text{（一階收斂）}. \]

「局部 \(h^2\) 、全局 \(h\) 」——平方與一次之間的差就是「誤差傳播」的收費：Grönwall 不等式（3.2 節）的離散版是全局誤差分析的唯一工具（8.1 節場景、6.3 節的連續版在此會師）。

隱式 Euler（反向）： \(y_{n+1} = y_n + h f(t_{n+1}, y_{n+1})\) ——需要解方程（每步一次 Picard/Newton 疊代，3.1 節的戰場），但穩定性質大幅改善（見剛性段）。

Runge–Kutta：採樣斜率的藝術#

Martin Kutta（1867–1944）與 Carl Runge（1856–1927） 的想法：一個步長內多採幾個斜率、加權平均，使 Taylor 展開匹配到更高階。經典四階（RK4，1901 Kutta）：

\[ \begin{aligned} k_1 &= f(t_n,\ y_n), \\ k_2 &= f(t_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2} k_1), \\ k_3 &= f(t_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2} k_2), \\ k_4 &= f(t_n + h,\ y_n + h k_3), \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{aligned} \]

權重 \(\frac{1}{6}, \frac26, \frac26, \frac16\) ** 是 Simpson 求積公式（積分近似 \(\int \approx \frac{h}{6}(f(a) + 4f(\frac{a+b}{2}) + f(b))\) 的影印——RK 的本質是把 ODE 寫成積分方程（3.1 節！）後用高階求積）。局部誤差 \(O(h^5)\) 、全局 \(O(h^4)\) **——步長減半，誤差除以 16。

「每步四次函數調用換 16 倍精度」——RK4 是整個科學計算中使用最廣的單一演算法。歷史注記：Runge（物理學家、氪的發現者之一）1895 年提出二階版本；Kutta（應用數學家、空氣動力學 Kutta–Joukowski 定律的另一半）1901 年完成 \(s = 4\) 階分類。RK 方法的一般理論（Butcher 階條件、1987 Butcher 樹）到二十世紀後期才完備——「為什麼恰好是這些權重」的完整答案用了近百年。

收斂階對照表#

方法 每步 \(f\) 調用 全局階 半步長誤差比
Euler（顯式） 1 1 2
中點法（RK2） 2 2 4
RK4 4 4 16
隱式 Euler \(1 +\) 疊代 1 2（但 A-穩定）
Adams–Bashforth（多步） 1（複用） 可高階 —
Dormand–Prince（RK45, 1980） 6 5（內嵌 4 階誤差估計→自適應步長） —

現代實務：自適應步長（誤差估計 + 步長調整，Fehlberg 1969 的 RK45 是 MATLAB ode45 的內臟）——讓計算機自己決定步長，平滑處走大步、陡峭處走小步。

剛性：穩定性的深水區#

系統 \(\dot y = -\lambda_1 y + \cdots\) 、 \(\dot z = -\lambda_2 z + \cdots\) （ \(\lambda_1 = 0.1\) 、 \(\lambda_2 = 10000\) ）：快模態（ \(z\) ）幾微秒衰減殆盡，但顯式 Euler 的穩定條件（對 \(\dot y = \lambda y\) ， \(|1 + h\lambda| \le 1\) ⇒ \(h \le \frac{2}{|\lambda|}\) ）強制 \(h \le 0.0002\) ——為一個早已死亡的模態，全程用殺手步長。剛性 = 模態時間尺度懸殊（3.3 節場景四）。

A-穩定性（Dahlquist 1963）：方法的穩定域包含整個左半平面。隱式 Euler 是 A-穩定的： \(y_{n+1} = \frac{y_n}{1 - h\lambda}\) ， \(|1 - h\lambda|^{-1} \le 1\) 對一切 \(\operatorname{Re}\lambda < 0\) ✓——任意大步長都穩定（精度仍受限，但可以「正確地粗略」）。梯形法（隱式中點）同樣 A-穩定且二階。剛性問題（化學動力學、電路暫態、氣候模式）用隱式方法——每步多花的疊代成本，換來步長解放。

Stiff solvers 的工程地位：MATLAB 的 ode15s、化學反應器模擬（大氣化學的速率方程百餘種、時間尺度從納秒到年）、電力系統暫態——現代科學計算的一半精華是剛性處理（BDF 方法，Curtiss–Hirschfelder 1952、Gear 1971 標準化）。

一個具體例題#

例題 10.1.1：用 Euler 法與 RK4 解 \(y' = -2y + 1\) 、 \(y(0) = 0\) ， \(h = 0.1\) ，步進兩步至 \(t = 0.2\) ；與真解 \(y = \frac{1 - e^{-2t}}{2}\) 比較誤差，驗證階數關係。

解：

Euler： \(f(t, y) = -2y + 1\) 。

  • 步 1： \(y_1 = 0 + 0.1 \times (0 + 1) = 0.1\) 。
  • 步 2： \(y_2 = 0.1 + 0.1 \times (-0.2 + 1) = 0.1 + 0.08 = 0.18\) 。

RK4（步 1， \(t = 0, y = 0\) ）：

\[ k_1 = 1,\quad k_2 = f(0.05, 0.05) = -0.1 + 1 = 0.9,\quad k_3 = f(0.05, 0.045) = 0.91,\quad k_4 = f(0.1, 0.091) = 0.818. \]
\[ y_1 = 0 + \frac{0.1}{6}(1 + 1.8 + 1.82 + 0.818) = \frac{0.1}{6} \times 5.438 = 0.090633\ldots \]

（步 2 略，同法： \(y_2 \approx 0.16555\) 。）

真解： \(y(0.1) = \frac{1 - e^{-0.2}}{2} = 0.0906346\ldots\) ； \(y(0.2) = \frac{1 - e^{-0.4}}{2} = 0.164840\ldots\)

誤差對照：

Euler RK4
\(y(0.1)\) 誤差 \(0.1 - 0.09063 = 0.00937\) $
\(y(0.2)\) 誤差 \(0.18 - 0.16484 = 0.01516\) \(\approx 7 \times 10^{-7}\)

比例檢驗：RK4 誤差約為 Euler 的 \(\frac{1}{5000}\) ——四次函數調用的代價換來四個數量級的精度。若步長減半（ \(h = 0.05\) ）：Euler 誤差約減半（一階），RK4 誤差約除以 16（四階）——階數由「步長—誤差對數斜率」實測（數值實驗的標準驗收程序：取 \(h, h/2, h/4\) 畫對數圖，斜率即階）。∎

歷史啟示#

  • 離散即工程：Euler 折線法的誕生（1768）比嚴格存在性理論（1890）早一百二十年——計算的需求再次先行於理論（呼應 8.1 節 Heaviside）。
  • 階數的經濟學：從 1 階到 4 階的跨越，本質是在每步內部做更聰明的 Taylor 匹配——「多採樣、重加權」是所有數值學科的通用升級（求積、插值、譜方法）。
  • 剛性教會我們尊重時間尺度：顯式方法的失效不是計算不夠快，而是穩定域的幾何——理解「為什麼慢」，比「算得快」更重要（20 世紀數值分析的核心教訓）。

與其他章節及領域的關聯#

  • 3.1、3.2 節：Picard 積分方程 = 隱式方法每步的子問題；Grönwall = 全局誤差引擎。
  • 7.1 節：振盪系統的數值能量漂移（辛積分器 6.1 節場景五）。
  • 9.1、9.2 節：PDE 離散化 = 空間差分 + 時間步進（CFL 9.2 節場景四）；每一個 PDE 求解器都是一個 ODE 積分器。
  • 8.2 節：Laplace 域看剛性（極點模差）。
  • 10.2 節：應用模型的數值實現（流行病、藥動、生態）。
  • 航天：JPL 星曆的 Taylor 級數積分器（5.1 節場景四）；太陽系數值積分需保能量（辛/對稱方法）。
  • 機器學習：神經 ODE（2018）的反向傳播 = 隨伴法（adjoint method，1960s 控制論遺產）；訓練動態 = 剛性動力系統。

練習#

  1. 用 Euler 法（ \(h = 0.25\) ）解 \(y' = x + y\) 、 \(y(0) = 1\) 至 \(x = 0.75\) ，與真解 \(y = 2e^x - x - 1\) 比較。
  2. 證明：對 \(\dot y = \lambda y\) （ \(\lambda < 0\) ），顯式 Euler 的數值解為 \(y_n = (1 + h\lambda)^n y_0\) ，穩定需 \(h \le \frac{2}{|\lambda|}\) ；隱式 Euler 給 \(y_n = (1 - h\lambda)^{-n} y_0\) 恆穩。
  3. 對練習 1 的方程，用中點法（RK2： \(y_{n+1} = y_n + h f(t_n + \frac h2, y_n + \frac h2 f_n)\) ）重算兩步，比較誤差。
  4. 實驗設計題：寫出用 \(h, h/2, h/4\) 三步長測量某方法收斂階的程序（比值 \(\frac{e(h)}{e(h/2)} \approx 2^p\) ），並預測 RK4 的應有比值。
  5. 思考題：剛性系統 \(\dot y = -1000(y - \cos t) - \sin t\) 的真解光滑（ \(y \approx \cos t\) ）。請解釋：為什麼「解光滑」不意味「顯式方法好用」？這一點對「AI 訓練動態 = 微分方程」的觀點（10.3 節）有什麼警示？

進階研讀：數值 ODE 的五個深化場景#

場景一：辛積分器與長期保真。哈密頓系統的 Euler 顯式法使能量漂移；辛 Euler（ \(p_{n+1} = p_n - h H_q\) 、 \(q_{n+1} = q_n + h H_p(p_{n+1})\) ——交錯更新）保持辛結構，能量誤差有界（不漂移）。分子動力學（速度 Verlet，1967）與太陽系長期模擬（Laskar 1989–）的基石。「尊重幾何的步進」——數值方法的幾何覺醒（6.1 節場景五）。

場景二：誤差的統計視角。捨入誤差（浮點）與截斷誤差的競賽：步長過小時捨入誤差累積（ \(10^{-16}\) 每步 × \(N\) 步）反超截斷誤差——最優步長存在（總誤差 U 型曲線）。Wilkinson 1960s 的浮點誤差理論使「算錯」變成「可控」——數值分析是「誤差的物理學」。

場景三：延遲微分與時滯系統。 \(\dot y = f(y, y(t - \tau))\) （時滯 ODE，流行病潛伏期、控制回路的傳輸延遲）：需要歷史函數而非初始點；時滯使二維系統也能產生振盪與混沌（Mackey–Glass 1977 血細胞模型）。數值上 = 插值歷史 + 步進——記憶進入狀態，維數悄然無窮化。

場景四：驗證文化——守恆量檢查。工程模擬的驗收不是「圖好看」，而是守恆量漂移率（能量、動量、質量）與收斂階測試（Richardson 外推 1927：由兩層步長解外推精確解）。數值方法的「單元測試」——每個模擬團隊的日常儀式（CFD 的 benchmark、氣候模式的 Genesis 測試）。

場景五：從步進到學習。神經 ODE 把網絡層視為連續流 \(\dot z = f_\theta(z, t)\) ，訓練 = 對 \(\theta\) 的梯度流；adjoint 法（Pontryagin 1956 的伴隨方程）算梯度免去儲存全部中間層。強化學習的 policy 演化、擴散模型（diffusion models，2020s，生成式 AI 的骨幹）本質是學會反轉一個隨機微分方程的時間方向——數值積分器成為 AI 的推理引擎（每張 Stable Diffusion 圖片都是一次 RK/Heun 積分）。

補充：Runge 其人——化學家、光譜學家、數值先驅#

Carl Runge（1856–1927）是 Bunsen 的學生、與 Rayleigh 共事，1895 年與 Paschen 發現稀有氣體氪（光譜分析）。他的數值方法論文產生於光譜波長的最小平方擬合實務。1895 年論文同時給出「三係數求積」與「二階 RK」——實驗數據處理的需求催生了數值積分的規範化。他晚年主持哥廷根應用數學與航空研究（一戰的彈道計算）——數值分析的誕生始終伴隨著測量與工程。

補充：為什麼「捨入」毀掉了 Lorenz 的重跑（1961）#

Lorenz 用六位小數輸入、機器內部六位計算——重跑輸入 0.506 而原值 0.506127：初值誤差 \(1.3 \times 10^{-4}\) ，經 \(\lambda \approx 0.9\) 的指數放大（7.3 節），數十個模擬「天」後完全脫鉤。捨入誤差的命運由動力系統決定：穩定系統中捨入誤差安分（6.3 節），混沌系統中捨入誤差是「蝴蝶的翅膀」。數值精度與動力學性質的耦合——這是數值分析與動力系統的交界課題（shadowing theory，1970s：儘管單軌道漂移，數值軌道仍「影子跟蹤」某條真軌道——混沌模擬的救贖）。

名詞速查表（本節出現的核心術語）#

  • Euler 法：折線步進；一階。
  • 局部/全局截斷誤差： \(O(h^{p+1})\) / \(O(h^p)\) （ \(p\) 階方法）。
  • Runge–Kutta：多級斜率採樣；RK4 四階經典。
  • Simpson 求積聯繫：RK = 積分方程 + 高階求積。
  • 隱式方法：每步解方程；A-穩定。
  • 剛性：時間尺度懸殊；顯式步長被快模態綁架。
  • A-穩定性：穩定域含左半平面（Dahlquist 1963）。
  • 自適應步長：內嵌誤差估計 + 步長調控（RK45）。
  • CFL 條件：數值波速 ≥ 物理波速（PDE 離散化）。
  • 辛積分器：保幾何結構的長期穩定方法。
  • 伴隨法（adjoint）：梯度計算的反向積分；機器學習回傳的數學內臟。

年表速查（本節涉及的關鍵年份）#

  • 1768：Euler《積分法》折線法。
  • 1895：Runge 二階方法。
  • 1901：Kutta 四階與一般分類。
  • 1926：CFL 條件（1928 發表）。
  • 1927：Richardson 外推。
  • 1946–1950：von Neumann 穩定性分析（MANIAC 時代）。
  • 1952：BDF 思想（Curtiss–Hirschfelder）。
  • 1963：Dahlquist A-穩定性。
  • 1967：速度 Verlet（分子動力學）；Gear 標準軟體（1971）。
  • 1969：Fehlberg RK45 自適應。
  • 1980s：辛積分器理論（Hairer、Sanz-Serna 等）。
  • 2018：神經 ODE。

深化問答三則#

問：為什麼不直接用很小步長 + 簡單方法？ 答：三重障礙。（1）時間： \(10^9\) 步 × 每步 \(10^{-6}\) 秒 = 17 分鐘起跳，嵌在 PDE 與反覆試錯中不可承受；（2）捨入：累積誤差 \(\sim N \varepsilon_{\rm machine}\) 反噬（步長過小的 U 型陷阱）；（3）剛性：顯式方法在小步長下仍可能不穩（步長再小只要 \(> \frac{2}{|\lambda_{\max}|}\) 就炸）。高階方法不是奢侈品，是工程可行性。

問：RK4「無條件好用」嗎？ 答：不。三大禁區：（1）剛性——RK4 穩定域有限，剛性系統需隱式；（2）長期哈密頓模擬——能量漂移，需辛方法；（3）強非光滑——高階方法的 Taylor 假設需要 \(f\) 光滑，不連續輸入下 RK4 退化為低階。方法與問題的匹配度 > 方法本身的階數——數值分析的第一智慧。

問：數值解「正確」是什麼意思？ 答：現代標準三件套：（1）收斂性（步長→0 時解→真解，有理論階數）；（2）守恆性（不變量的漂移可控）；（3）驗證性（與已知解/實驗對照，Richardson 外推檢查）。注意「與實驗吻合」不保證正確（可能兩錯相抵）——V&V（驗證與確認）文化，從核工程到自動駕駛模擬，是計算科學的品質體系。數值解是假說，收斂階與守恆檢查是它的審判程序。

（本節完）