別一口氣把所有東西都離散化
這一級到目前為止,都把空間與時間放在同等地位:你在兩者上都鋪了格點,你把每一個導數——時間導數與空間導數一視同仁——都換成一個樣板,於是冒出像 FTCS 或克蘭克-尼科爾森這樣單一的配方,告訴你如何從當前時間層得到整個下一個時間層。這行得通,但它把兩件很不一樣的工作捆進同一個公式裡,並逼你為每一個新的時間步進想法重新推導一套新格式。線方法(MOL) 是一種把這兩件工作解開的態度轉變。它的口號很簡單:離散化空間,但暫時別碰時間。
拿一根桿上的熱方程 u_t = alpha * u_xx。只鋪空間格點:點 x_1, x_2, ..., x_M 間距為 h,兩端邊界固定住。在每個內部點定義 u_i(t),即那個固定位置的溫度作為時間的函數。現在做半離散化。把時間導數 u_t 原封不動留著——它是一個如假包換的對 t 的導數。但把空間導數 u_xx 換成中心二階差分樣板 (u_{i-1} - 2*u_i + u_{i+1}) / h^2。點 x_i 處的偏微分方程就變成 du_i/dt = (alpha/h^2) * (u_{i-1} - 2*u_i + u_{i+1})。那不再是偏微分方程了。它是一個對 t 的常微分方程,每個格點一個。
想想這名字為何貼切。每個內部格點坐落於一個固定的 x,並把它的溫度沿時間往前描出——一條橫貫時空平面的水平線。偏微分方程被切成一束這樣的線,每個節點一條,而沿著每條線,規則如今是一個常微分方程。這些線並非彼此獨立:x_i 處的方程伸手抓向它的鄰居 x_{i-1} 與 x_{i+1},所以整束是一個巨大的、必須一起推進的耦合系統常微分方程。但那是上一級的問題,不是一個新問題。
偏微分方程變成一個巨大的常微分方程系統
把未知數疊成一個向量 u(t) = (u_1(t), u_2(t), ..., u_M(t))——時刻 t 時整根桿的溫度分佈,凍結成一串數字。中心差分樣板在每個內部點的作用方式相同,於是那束常微分方程塌縮成一個工整的矩陣方程:du/dt = A u,其中 A 是熟悉的三對角二階差分矩陣,對角線是 -2、兩側各是 +1,再乘上 alpha/h^2(一維的離散拉普拉斯算子)。這恰恰是一個系統常微分方程的初值問題——與常微分方程那一級的 y' = f(t, y) 同樣的形狀,只不過 y 是一個長向量、f 是線性映射 A。初始條件是起始溫度分佈 u(0);邊界值作為已知常數,餵進第一個與最後一個方程。
回報就在這裡,而且很大。你再也不需要為偏微分方程發明一個時間步進格式了——你伸手進常微分方程那一級,抽出任何你喜歡的求解器。把前向尤拉法用在 du/dt = A u 上,你會原樣復現 FTCS 格式;那不是巧合,那是 FTCS 真正的身分。用後向尤拉法你得到 BTCS。用梯形法則——步首與步尾斜率的平均——你得到克蘭克-尼科爾森。你當作三個分開的發明來認識的那三個格式,其實只是三個不同的常微分方程積分器,栓到同一個空間離散化上。而且沒什麼能阻止你栓一個更好的上去:一個龍格-庫塔法、一個自適應求解器、一個多步法,這些你都已經懂了。
PDE: u_t = alpha * u_xx on 0 < x < 1
step 1 (space): grid x_i = i*h, replace u_xx by central stencil
du_i/dt = (alpha/h^2) * (u_{i-1} - 2 u_i + u_{i+1})
stack into a vector u = (u_1, ..., u_M):
du/dt = A u A = (alpha/h^2) * tridiag(1, -2, 1)
step 2 (time): hand du/dt = A u to ANY ODE solver
forward Euler -> FTCS (explicit, cheap, CFL-limited)
backward Euler -> BTCS (implicit, solve A x = b each step)
trapezoidal -> Crank-Nicolson(implicit, 2nd order in time)為何 CFL 限制再度現身——以及它其實就是剛性
線方法不只是省力:它從單一個制高點解釋了前面幾篇的穩定性結果。一個線性系統 du/dt = A u 在某個常微分方程求解器下穩定,恰好當 A 的每一個特徵值 lambda,其縮放後的值 h_t * lambda(時間步長乘以特徵值)落在那個求解器的絕對穩定域之內——正是你在常微分方程那一級研究過的那個域。這裡 A 是實對稱的,一個對稱特徵值問題,所以它的特徵值全是實數,這讓圖像保持在實軸上。於是「我的偏微分方程格式穩定嗎?」這個問題,變成了「A 的特徵值落在哪裡,而我的時間步長有沒有把它們縮進穩定域裡?」同一張圖,用兩次。
現在替熱方程算一下。二階差分矩陣 A 的特徵值範圍從接近 0 一直到大約 -4*alpha/h^2——全是實的負數,而最負的那個隨你細化格點像 1/h^2 那樣增長。前向尤拉法的穩定域沿負實軸只伸到距離 2,所以穩定性要求 h_t * (4*alpha/h^2) <= 2,也就是 h_t <= h^2 / (2*alpha)。FTCS 那一篇的顯式時間步上限,就在這裡,從特徵值直接掉了出來——而那個懲罰人的 h^2 一點也不意外,一旦你看出它來自 A 裡的那個 1/h^2。後向尤拉法與克蘭克-尼科爾森的穩定域覆蓋整個左半平面,所以無論 h_t 多大,那些負特徵值都逃不出去:無條件穩定,正如隱式格式那一篇所宣稱的。
誠實的細則
線方法很強大,卻不是魔法,而這座學習階梯其餘部分的紀律同樣適用。首先,空間離散化帶著它自己的截斷誤差:中心二階差分只有 O(h^2) 的準確度,所以即使配上完美無瑕的時間積分器,你的答案在空間上也不會好過那個樣板所允許的。把時間上的高階常微分方程求解器配上空間上的低階樣板是白費力氣——準確度由兩者中較弱的那個決定。其次,每一步仍在浮點數裡運行,那裡 0.1 沒有精確的二進位形式,所以永遠縮小 h 並不會一直有幫助:過了某個點,1/h^2 樣板裡的捨入誤差會占主導,這與限制任何有限差分導數的那道捨入誤差地板是同一回事。
第三,選擇隱式時間積分器並不會讓工作消失——它把工作搬了個地方。後向尤拉法與克蘭克-尼科爾森要求在每個時間步解一個線性系統 A x = b,而你很少靠形成矩陣的逆來做這件事。在一維裡 A 是三對角的,所以 托馬斯演算法以 O(M) 次浮點運算解它,划算得很。在二維或三維裡,離散拉普拉斯算子是稀疏的、卻不再是三對角的,於是你會去取一個配上好預條件子的迭代求解器,或一個多重網格法——迭代法那一級的工具——正是因為在大規模下,形成或求逆那個矩陣是想都別想的。
再來一個誠實的但書。乾淨的特徵值故事,對簡單格點上的線性、常係數問題妥妥地成立;對非線性偏微分方程或狂野的係數,A 變成一個依賴狀態的映射,分析就是局部且近似的,是個嚮導而非保證。而對純對流問題如 u_t + a*u_x = 0,半離散的特徵值是純虛數而非負實數,這改變了哪些常微分方程求解器能保持穩定,並重新引入數值色散——提醒你線方法組織了分析、卻沒讓你免於去做分析。
你帶出這一級的東西
退一步,看這些拼塊如何彼此咬合。前面幾篇教你鋪格點、把導數變成樣板;權衡顯式對隱式;檢驗一個格式的相容性、穩定性與收斂性,並信賴那個「相容性加穩定性等於收斂性」的拉克斯等價定理;以及用馮諾依曼分析與 CFL 條件找出穩定的步長。線方法是那個統一的想法,讓這一切感覺像一個學科、而非四個花招:離散化空間得到一個常微分方程系統,然後把整套常微分方程工具箱——求解器、階、穩定域、剛性——搬來對付時間。
它也指向外頭。有限差分是離散化空間最溫和的方式,卻不是唯一的:下一級的有限元素(伽遼金)法與有限體積法產生一個不同的 A——以及左邊一個質量矩陣 M,給出 M*du/dt = A u——然而配方的後半截毫無改變,因為一旦你有了一個時間上的常微分方程系統,你早就知道該怎麼辦。譜方法用一個更稠密、更準確的 A 做同樣的事。線方法是「我該如何在格點上表示一個函數?」與「我該如何讓它隨時間前進?」之間的那道鉸鏈,而世上幾乎每一個正經的偏微分方程程式,都是繞著那道鉸鏈造出來的。