當小步長是被「錯誤的理由」逼出來的
到現在，你已經會用好幾種方式把解往前推：尤拉法粗糙的切線一跳、RK4四次取樣的精修、預測—修正法回收歷史的精打細算。它們全都帶著同一個承諾——把步長 h 縮小，局部截斷誤差縮得更快，所以較小的 h 換來更高的精度。這篇要談的是更棘手的情況：明明真解平滑無趣、大步長就綽綽有餘地準確，你卻被迫用荒謬的極小步長。這時步長付的錢不是買精度，而是買「讓方法別爆掉」。
這裡有個最乾淨的試管。取線性衰減 y' = -50 y，y(0) = 1。精確解是 y = e^(-50 x)，一條先急速下墜、再幾乎貼著零拉平的曲線——溫馴至極。現在對它跑前向尤拉。每一步做 y_new = y + h * (-50 y) = (1 - 50 h) y。這個數 1 - 50 h 是一個放大因子：每走一步就把答案乘上它一次。若 h = 0.1，則 1 - 50 h = -4，於是你的數值解變成 1, -4, 16, -64, 256, ……——一邊振盪一邊爆向無窮，而真正的答案卻在悄悄歸零。方程本身毫無問題。慌張發作的是方法。
y' = -50 y, y(0) = 1 true: y = e^(-50x) -> 0 forward Euler: y_(n+1) = (1 - 50 h) y_n h = 0.1 factor = -4.0 1, -4, 16, -64, ... EXPLODES h = 0.05 factor = -1.5 1, -1.5, 2.25, ... EXPLODES h = 0.04 factor = -1.0 1, -1, 1, -1, ... oscillates forever h = 0.02 factor = 0.0 1, 0, 0, 0, ... (lucky) decays h = 0.01 factor = 0.5 1, 0.5, 0.25, ... decays, stable
穩定區域：一個方法的安全地帶
替這個模式起個名字。把一個單步法用在模型方程 y' = lambda y 上，其中 lambda 是常數（此處為負，因為解在衰減）。走一步就把 y 乘上某個放大因子 R，而 R 只依賴於組合 z = h * lambda。對前向尤拉，R(z) = 1 + z。當且僅當 |R(z)| <= 1 時，方法保持有界——不會爆掉。使其成立的那組複數 z，就是該方法的絕對穩定區域，即穩定區域。它是 z 平面上一個固定的形狀，是演算法本身的性質，與你手上的特定方程無關。
至此整個故事豁然開朗。要得到衰減的真解，你需要 lambda < 0；要讓方法保持有界，你需要 h * lambda 落進那個固定的區域裡。對前向尤拉，|1 + z| <= 1 是一個圓心在 z = -1、半徑為 1 的小圓盤，因此在實軸上它要求 -2 <= h * lambda <= 0，亦即 h <= 2 / |lambda|。代入 lambda = -50 就是 h <= 0.04——恰好是我們看到的門檻。RK4 的穩定區域更胖、有瓣狀外凸，因此能容忍稍大一點的 h，但它依然是個有界區域：每一個顯式方法都有一塊有限的安全地帶，而一個很大的負 lambda 終究會把 h 逼到趨近於零。
剛性：快的成分早已消失，卻仍在發號施令
單獨一個衰減指數很少會困住你；真正的痛苦來自混合。設想一個系統的解大致長成 y = e^(-1 x) + e^(-1000 x)。第二項在 x = 0.01 之後就消失在雜訊裡；其餘區間的答案基本上只剩溫和的 e^(-x)。然而要讓顯式方法保持穩定，你卻必須讓 h * (-1000) 留在穩定區域內，把 h 逼到 0.001 左右甚至更小——而且是整段都得如此，遠在那個快速項死透之後。你正用上百萬個小碎步去描一條口袋計算機就能勾勒的曲線。這種失衡就是剛性：時間尺度相差懸殊，由最快的（早已衰減的）模態替所有人決定步長。
剛性很難下精確定義——並沒有單一數值會「啪」地一聲跳起來告訴你。一個好用的工作描述是：若在某區間上，純粹為了穩定，顯式方法被迫採用遠比「解的平滑程度所能正當化」更小的步長，這問題就是剛性的。化學家與工程師舉目所及皆是剛性系統：含快慢反應的化學動力學、小電容緊鄰大電感的電路、快速致動器配上緩慢動態的控制系統。誠實的結論是：在剛性問題上，RK4 的精度白白浪費了——卡住你的是它有限的穩定區域，而不是它的誤差。
隱式解藥：後向尤拉與 BDF
如果問題出在穩定區域太小，解藥就是換一個穩定區域超大的方法。把斜率從「在舊點取值」改成「在新點取值」：後向（隱式）尤拉寫成 y_(n+1) = y_n + h * f(x_(n+1), y_(n+1))。在模型 y' = lambda y 上，這給出 y_(n+1) = y_n / (1 - h * lambda)，於是放大因子是 R(z) = 1 / (1 - z)。對任何實部為負的 z，|R(z)| < 1 自動成立——穩定區域是整個左半平面。這個方法是 A-穩定 的：無論負 lambda 多大、步長 h 多大，計算都保持有界且衰減，正如真解一般。那個剛性項再也無法挾持你。
天下沒有白吃的午餐，只有換一張帳單。由於 y_(n+1) 同時出現在等號兩邊，隱式方法的每一步都是一道「要解的方程」，而非一條「拿來算的公式」——對非線性的 f，你每一步都得跑一輪牛頓迭代。這使得每一步更昂貴。但在剛性問題上，這是一筆極划算的交易：後向尤拉也許能走 h = 0.1，而顯式 RK4 卻被釘在 h = 0.001，於是縱使每步要幾輪牛頓迭代，它仍以數量級之差勝出。實用的剛性求解器採用後向差分公式（BDF），亦即後向尤拉的高階表親，既保有大穩定區域，又把精度提升到一階以上。
自適應步長：讓求解器自己挑選步調
還有第二個讓固定步長浪費的理由，它和穩定性毫無關係。真實的解在平緩處慢吞吞地晃，到了急轉彎處又猛地一甩；小到足以應付急彎的步長，用在平緩處是殺雞用牛刀，而按平緩處設定的步長到了急彎處又會搞砸。解法是自適應步長控制：在每一步估計局部誤差，若它舒舒服服地低於你的容許值，就接受這一步並把 h 放大；若它超標，就拒絕這一步、把 h 縮小，再重試。求解器這樣一路摸索著走，在曲線偷懶處大步流星，在它忙碌處則碎步而行。
可是，在不知道真答案的情況下，你要怎麼估計誤差呢？那個漂亮的技巧是一對嵌入式公式：跑兩個不同階的龍格—庫塔公式，並刻意讓它們共用斜率的取值，於是第二個答案幾乎是白送的。它倆的差就是局部誤差一個廉價又可靠的估計。RKF45（Runge-Kutta-Fehlberg）做的正是這件事——它把一個 4 階與一個 5 階的公式配成一對、重用同一批取樣的斜率，再用兩者之間的落差來驅動步長控制器。MATLAB 裡著名的 ode45 與 SciPy 預設的 RK 方法，都是這個想法的現代後裔。
- 從當前點出發，用嵌入式公式對（一個 p 階、一個 p+1 階）算出下一個值的兩個估計。
- 把兩者的差當作局部誤差估計 e，並與你的容許值 tol 比較。
- 若 e <= tol，接受這一步（採用高階的那個值）並把下一次的 h 放大；若 e > tol，拒絕它、縮小 h，再重做這一步。
- 用形如 h_new = h * (tol / e)^(1/(p+1)) 的公式（再乘一個安全係數）來縮放新步長，讓控制器平緩地朝著目標誤差靠攏。
把它們串起來——以及還有什麼會咬人
這兩個想法——穩定性與自適應——並非對手；最好的求解器把它們編織在一起。自適應控制器藉由在尖銳特徵附近縮小 h，保護你免於精度的災難；而一個隱式核心（BDF、A-穩定）則在剛性區段保護你免於穩定性的災難。這正是成熟的函式庫為何同時提供一個剛性求解器與一個非剛性求解器，有些甚至會嗅探問題並在兩者間自動切換。數值常微分方程的功夫，大半在於把方法配上方程的脾性，而不是抓著某個心頭好不放。
對剩下的陷阱也要誠實。自適應控制約束的是「每一步」的局部誤差，而非累積的全域誤差——微小的每步誤差在長程積分中仍會累加，再緊的容許值也不保證一條全域忠實的曲線。把容許值壓得太低，你會撞到地板：有限精度算術的捨入誤差開始主導，於是再縮 h 只會把事情弄得更糟而非更好。而且沒有任何步進法能拯救一個不適定的問題——若方程本身具有敏感依賴（混沌），或有一個在有限時間內爆破的解，那些數字會忠實地追蹤那份不穩定，一路追到它再也追不下去為止。這些方法是誠實的鏡子；它們映照出方程真實的性格，連瑕疵也一併映出。