JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

顯式與隱式格式

同一條熱傳導方程,有兩種向前推進時間的方式:一種每步便宜卻被一個極小的步長上限綁死,另一種每步較貴卻掙脫了那條鎖鏈。認識 FTCS、BTCS 與 Crank-Nicolson,並理解為什麼「解一個線性方程組」反而可能是更快的選擇。

岔路口

在上一篇指南裡,你把一張網格鋪在一根金屬棒上,並用有限差分模板取代導數:某一點的空間二階導數 u_xx 變成 (u_{i-1} - 2 u_i + u_{i+1})/h^2,也就是一個值與它兩個鄰居的差,再除以間距的平方。現在我們把時間加進來。熱傳導方程說 u_t = alpha u_xx——溫度隨時間的變化率,等於常數 alpha 乘上它在空間中的曲率。我們已經知道怎麼離散化右邊。本篇的整個問題,就在一個看似無辜的選擇上:我們要在哪個時間層上計算那個右邊——是我們正要離開的時間層,還是我們正要抵達的時間層?

這正是你在上一級、處理常微分方程時遇過的同一個岔路。在那裡,前向歐拉法用的是當前點的斜率——一個顯式、代進去就走的步——而後向歐拉法用的是終點的斜率,把未知量同時困在方程式的兩邊。PDE 的情形完全繼承了這套性格分裂,只是現在「未知量」不再是一個數字,而是棒子上整整一排的網格值。就是這一個差別——一個未知量,相對於一整排彼此耦合的未知量——把一次平凡的更新變成了一個線性方程組,而它正是接下來一切的核心。

FTCS:顯式的走法

把右邊取在「舊」的時間層 n。時間上用前向差分、空間上用中央差分,你就得到 FTCS 格式——Forward-Time, Centred-Space(前時、中央空間),也就是 FTCS 格式。把它解出那個唯一的新值 u_i^{n+1},右邊的一切都早已從你手上那一排舊值得知。每個新溫度,不過是舊溫度被鄰居的加權平均輕輕推了一下。沒有方程組要解,沒有矩陣要反;你橫掃過棒子一次就完成了。這就是顯式的意思:新值被明明白白地用舊值寫出來。

let r = alpha * dt / h^2     (the key ratio)

FTCS update for each interior i:
    u_i^{n+1} = u_i^n + r * (u_{i-1}^n - 2 u_i^n + u_{i+1}^n)

= (1 - 2r) u_i^n + r u_{i-1}^n + r u_{i+1}^n
FTCS:一次顯式橫掃。新值是三個舊鄰居的加權混合,由比值 r 控制。

便宜又簡單——那陷阱在哪?看看中心點的係數 (1 - 2r)。如果比值 r = alpha dt/h^2 大於 1/2,那個權重就變成負的。從物理上說,熱會把一個分佈抹平;但在數值上,一個負的中心權重做的恰恰相反——它放大擾動,一個微小的凸起會長成一片正負來回擺盪、爆炸的棋盤格。FTCS 只是條件穩定的:唯有當 r 至多為 1/2 時它才管用。完整的論證要等到兩篇之後的 von Neumann 分析,但你現在已經能在骨子裡感受到它的後果。

BTCS:隱式的走法

現在改把空間模板取在「新」的時間層 n+1。這就是 BTCS 格式——Backward-Time, Centred-Space(後時、中央空間),也就是 BTCS 格式。u_i^{n+1} 的更新式現在牽涉到新層上的三個未知量:u_{i-1}^{n+1}、u_i^{n+1} 與 u_{i+1}^{n+1}。你沒辦法孤立地解出單一個網格點,因為每個新值都倚靠它的新鄰居,而那些鄰居又倚靠它們的鄰居,一路橫貫整根棒子。新的一排只被隱式地定義出來——作為一組耦合方程式的解,每個內部點對應一條方程式。

把所有這些方程式收攏起來,它們就構成一個線性方程組 A u^{n+1} = u^n(邊界項已併入)。這裡的 A 不是一個稠密、嚇人的矩陣——每條方程式只牽涉一個點與它最近的兩個鄰居,所以 A 只有在主對角線與緊鄰的兩條對角線上才有非零元素。那是一個三對角矩陣,而它是個禮物。一般的 n×n 求解要花 O(n^3),但三對角的情形可以交給 Thomas 演算法——只碰那三條活躍帶的高斯消去法——僅以 O(n) 個運算解決,與網格點數成線性。一個隱式步的代價,只比一次顯式 FTCS 橫掃多出一個小小的常數倍。

而回報很大:BTCS 是無條件穩定的。無論你把 dt 取得多大,這個格式都不會爆炸——這和後向歐拉法在常微分方程那一級馴服一條剛性方程是同一個道理。FTCS 那個負權重陷阱在這裡根本不存在,因為隱式耦合在每一步都抑制擾動,而不是放大它們。你可以單純為了準確度來選 dt,而不必去討好穩定性警察。對於空間網格必須很細(h 很小)的問題,這份自由可能意味著少走數千倍的步數,於是即使每步裡藏著一次線性求解,隱式方法仍然全面勝出。

Crank-Nicolson:折衷取中

FTCS 與 BTCS 共有一個 von Neumann 那篇會量化的缺點:兩者在時間上都只有一階準確度,誤差 O(dt)。把時間步長減半,時間誤差也只減半——當空間誤差已經是 O(h^2) 時,這是個糟糕的交易。修補方法很優雅:把空間模板取在兩個時間層的正中間,當作新舊兩個模板的平均。那就是 Crank-Nicolson 格式,也就是 Crank-Nicolson 格式——FTCS 與 BTCS 各半的混合。它仍然是隱式的(新層依舊出現),所以你每步仍要解一個三對角方程組,但這個對稱平均把時間差分置中在中點,將準確度提升到 O(dt^2)。

Crank-Nicolson 既無條件穩定,又在空間與時間上皆為二階——這正是它成為熱傳導方程主力、並在無數求解器中當預設的原因。但誠實要求我們讀小字。它是 A-穩定但非 L-穩定:當 dt 非常大時,它的放大因子趨近 -1 而非 0,於是初始資料裡的一個尖銳跳變不會被消滅,只是每步被翻一次正負號,產生像被敲響的鐘一樣緩慢衰減的振盪。能量永遠不會爆炸(它保持穩定),但答案看起來可能很醜。實務上的解法是先用一兩個 BTCS 步把尖銳的暫態壓垮,再切換成 Crank-Nicolson 處理之後平滑的部分。

誠實地選一個格式

那麼你該伸手抓哪一個?沒有放諸四海皆準的贏家,只有幾條誠實的經驗法則。這個決定取決於統治整門學科的同一個三重奏:你需要的準確度、每步的代價,以及穩定性上限咬得有多兇。

  1. 如果網格很粗、而 dt 自然滿足 r <= 1/2,FTCS 就是能用的最簡單選擇——幾行程式碼,不需要線性求解器。原型與教學住在這裡。
  2. 如果空間網格必須很細(h 很小),顯式的 dt <= h^2/(2 alpha) 上限就會令人動彈不得。改用隱式:BTCS 或 Crank-Nicolson 讓 dt 為準確度而選,而 O(n) 的三對角求解讓每步維持便宜。
  3. 如果你想要二階準確度且資料平滑,Crank-Nicolson 是預設首選。如果初始資料有尖銳跳變或激波,先用幾個 BTCS 步把它壓掉,以避免那種振鈴。
  4. 永遠記得答案仍是近似的。每一步都在浮點數裡捨入,所以即使是在良態問題上的無條件穩定格式,也只交付一個近似值;而它的階 O(dt^2) 或 O(h^2) 是關於趨勢的承諾,不是對位數的保證。

最後帶一個警告往前走。「無條件穩定」是「不會爆炸」的承諾,不是準確度的承諾。用 BTCS 配上巨大的 dt,解會保持有界、行為良好,卻可能把一個移動的鋒面抹成一團糊——你拿一場爆炸換來了一片模糊。一個收斂到正確答案的格式,與一個只是維持有限的格式,這兩者之間的區別,正是下一篇談一致性、穩定性與收斂性的指南要去解開的東西。你在這裡做的顯式對隱式的選擇搭好了舞台;而 Lax 等價定理會告訴你它何時真正划算。