一道謎題：答案對了，步長卻像錯了
這裡有一幕小場景，困惑過一代又一代的初學者。你在解初值問題 y'(t) = -1000*(y - cos(t)) - sin(t)，從 y(0) = 0 出發，而真解很溫和：一眨眼之後它幾乎完全貼著 cos(t) 走，是一道平滑的波，變化速度從不超過 cos。於是你拿起第一篇那位可靠的主力，顯式尤拉法，挑一個像 h = 0.01 這種足以漂亮解析 cos(t) 的步長，然後按下執行。結果不是一點小誤差——而是一場災難。數字盪到一百、再到一萬，十幾步之內就溢位到無窮大。你追的那個解明明很平靜，方法卻爆炸了。
你的程式碼沒有錯，那個解也沒有錯。這場爆炸並不是你在第二、三篇遇過的那種精度失敗——那時步長 h 太大會給出一個可辨認、有界、能縮小的誤差。這是另一種野獸：一種與答案有多平滑毫無關係的不穩定。這條方程裡藏著一個快尺度——那個 -1000 的因子——它在最初一瞬就衰減殆盡、如今在物理上無關緊要，卻仍在號令你的顯式方法所被允許的步長。會這樣的方程有個名字，它是剛性的，而剛性正是這門學問始終在工具箱裡備著一整族隱式方法的頭號理由。
測試方程與絕對穩定區
要看清到底哪裡出錯，我們把問題剝到只剩骨架。把任何衰減行為線性化，就得到測試方程 y' = lambda*y，其中 lambda（常）是一個負實數——真解是 y(t) = y_0 * exp(lambda*t)，會衰減到零。每個方法真正在意的，是 lambda*h 這個數，也就是步長乘上速率。把顯式尤拉法套到測試方程：y_{n+1} = y_n + h*lambda*y_n = (1 + h*lambda)*y_n。所以每一步都把數值乘上因子 (1 + h*lambda)。真解每一步都在縮小，因此要讓數值解規矩，我們需要那個因子的大小至多為 1：|1 + h*lambda| <= 1。
讓某方法的成長因子留在單位圓盤之內的那一組 z = h*lambda 值（允許 z 為複數，因為真實系統也有振盪模態），就是該方法的絕對穩定區。對顯式尤拉法，條件 |1 + z| <= 1 是一個半徑為 1、圓心在 -1 的圓盤：一塊緊貼原點左側的、小而有界的區域。現在把我們的元兇代進去。lambda = -1000、h = 0.01，於是 h*lambda = -10，而 |1 + (-10)| = 9，遠在圓盤之外。每一步都把誤差乘上 9——這正是那場爆炸。要回到區域之內，我們需要 h*lambda >= -2，也就是 h <= 2/1000 = 0.002，硬被一個早已死去的瞬變所逼。
走向隱式：隱式尤拉法與它要付的代價
解圍之道在此。不要在你正要離開的點上計算斜率，而要在你正要抵達的那個點上計算——這就是隱式尤拉法：y_{n+1} = y_n + h*f(t_{n+1}, y_{n+1})。注意未知數 y_{n+1} 同時出現在等號兩邊；這個公式是隱式的，這一行就道盡了整個顯式對隱式的分別。把它套到測試方程 y' = lambda*y：y_{n+1} = y_n + h*lambda*y_{n+1}，重新整理成 y_{n+1} = y_n / (1 - h*lambda)。現在的成長因子是 1/(1 - h*lambda)。
看看這個因子在我們的剛性情形下做了什麼。h*lambda = -10 時，因子是 1/(1 - (-10)) = 1/11，舒舒服服地小於 1——這一步在衰減，不管步長有多大。事實上，對於每一個實部為負的 lambda，對所有 h > 0 都有 |1/(1 - h*lambda)| <= 1。它的絕對穩定區是整個複數平面的左半邊。穩定區包含整個左半平面的方法，稱為A 穩定，而 A 穩定正是剛性問題的黃金標準：它意味著步長上的穩定限制完全消失，你可以純粹為精度而挑 h。我們那道爆炸的問題，如今在 h = 0.01、甚至 h = 0.1 下都跑得開心。
沒有東西是免費的。隱式要付代價：每一步你都得解出 y_{n+1}，而不只是把它算出來。對我們那條線性測試方程，這是簡單的代數，但對真正的非線性系統 f(t, y)，它意味著每一步都要解一個（通常是非線性的）代數系統——而做這件事的工具，正是求根那一級的牛頓法推廣到系統的版本。所以一個剛性求解器，在每個時間步裡都藏著一輪牛頓迭代，而每個牛頓步都要解一個牽涉到 f 之雅可比矩陣的線性系統。這是實打實的工作量；一個隱式步可能要花顯式步的好幾倍成本。但這筆交易幾乎總是值得：你用一把昂貴卻被允許的大隱式步，換掉了成千上萬個被禁止的微小顯式步。
剛性究竟是什麼——以及如何辨認它
現在我們可以給剛性一個誠實的定義。一條剛性方程，是其解含有跨越天差地別尺度的衰減（或振盪）速率——一個幾乎瞬間消失的快瞬變，伴著一個你其實想追蹤的慢分量。對系統 y' = f(t, y)，相關的速率是雅可比矩陣的特徵值；剛性表現為最負與最不負的實部之間有一個很大的比值，即剛性比。我們那個純量例子假裝成一個系統：-1000 是快模態，cos(t) 是慢模態，比值約一千。化學動力學、同時含快與慢分量的電路，以及擴散型偏微分方程的空間離散化，都是經典的剛性來源。
有個常見誤解值得戳破：剛性並不等同於解變化得快。在最初的瞬變之後，我們的解平滑又緩慢到了極點，問題卻依然剛性。那個快尺度存在於方程之中，即便它已不在解之中——它支配著一個擾動會如何衰減，而顯式方法在每一步都得乖乖遵守它。這正是剛性殘酷的反諷：顯式方法被逼著邁出小到足以解析某個模態的步，而真解早已把那個模態忘得一乾二淨。一個好用的判症：如果顯式尤拉法或 RK4 要求一個荒謬地小的 h 才能保持有界，遠比答案的平滑程度所暗示的還小，那你幾乎肯定碰上了一個剛性問題。
BDF、更高階，與一個誠實的告誡
隱式尤拉法雖 A 穩定，卻只有一階，所以它太不準，當不了剛性問題的日常主力。實戰的答案來自第四篇那個多步族：後向差分公式（BDF）。一個 BDF 讓一條多項式穿過新點與若干個過去的點，再強迫它的導數在新的時間上匹配 f——和隱式尤拉法一樣是隱式的，且為剛性而生。BDF1 恰好就是隱式尤拉法；BDF2 在保持 A 穩定的同時達到二階；實務上會用到至多六階。這些就是著名剛性求解器內部的引擎，也是為何一個會自動偵測剛性的函式庫常式，一看見症狀就會悄悄切換到 BDF。
退一步，把整級放進同一個畫面看。第一篇給了你尤拉法與「跟著斜率走」的念頭；第二篇把局部誤差與整體誤差分開、並釘下了階；第三篇打造了每步高精度的 RK4；第四篇加上自適應步長控制與高效的多步法。這一篇補上了缺席的那個軸——穩定性——也補上了工具箱之所以是複數的理由。顯式方法（RK4、Adams-Bashforth）每步便宜，對非剛性問題堪稱完美；隱式方法（隱式尤拉法、BDF）每步要付一輪牛頓求解，卻能征服剛性。成熟的實踐者不挑最愛，而是診斷問題，再讓診斷去挑方法。也別忘了整門學問的根本誠實：求解器交回的每個數，都是有限精度算術裡的一個近似，所以穩定性為你買到的是一個有界、可信賴的答案——絕非一個精確的答案。