為什麼偏微分方程塞不進電腦
偏微分方程把一個量在時間與空間中的變化綁在一起。最乾淨的例子,也是我們整級的夥伴,就是一維熱傳導方程 u_t = a * u_xx:一根細桿上的溫度 u(x, t),其中 u_t 是隨時間的變化率、u_xx 是空間上的二階導數,而 a > 0 是熱擴散的快慢。直覺上它說:當一個點比它鄰居的平均還冷時就會升溫——熱往低處流。麻煩在於 u(x, t) 是連續的 x 與連續的 t 的函數,所以它帶著無窮多個值。沒有任何電腦能儲存一個在不可數無窮多個點上都有定義的函數。
於是我們做數值偏微分方程裡最古老的交易:我們同意只在有限的支架點上知道 u,其餘地方就退而求其近似。這就是貫穿本階梯每一級的先離散再求解策略——我們在把積分變成求和、把導數變成差商時都見過它。偏微分方程只是同時往兩個方向做這件事:空間與時間。本級的功夫就在於把這個替換做得讓離散配方真的能追上連續的物理,而非(我們將會看到)某個看似合理、卻悄悄爆掉的冒牌貨。
鋪設網格
挑一個空間步長 h 與一個時間步長 k,鋪下一張計算網格:點 x_j = j*h(j = 0, 1, ..., M)橫越細桿,時間層 t_n = n*k(n = 0, 1, 2, ...)往上攀升。把它想成方格紙:桿子由左到右、時間由下而上。每個交叉點坐著一個數,寫作 u_j^n,是我們對真實溫度 u(x_j, t_n) 的近似。上標 n 是時間層,不是次方——這個記號衝突每個人都會被絆一次,所以把 u_j^n 讀作「第 j 行、第 n 列的網格值」。我們全部的工作,就是一列接一列把這格點陣填滿,從一個已知的底列(初始溫度)與兩條邊緣行(邊界條件)出發。
把這格點陣具體地想像出來。底列 n = 0 完全已知——它就是在每一行上取樣的初始溫度 u(x, 0)。最左與最右兩行(j = 0 與 j = M)也對所有時間都被釘住:這些是邊界條件,比方說桿子兩端被冰固定住。內部的一切都是未知,而一個格式不過是一條規則,一次填一整列新值,從下方已經確定的列,算出第 n+1 列上的每個值。我們從不一次儲存超過兩三列;當我們往上攀,過去就從底部捲走了。
用模板取代導數
現在是關鍵的一步。導數是差商的極限;在間距固定的網格上我們無法取極限,所以就保留那個商。這正是微分那篇裡的有限差分公式想法,如今搬到空間上。對空間的二階導數 u_xx,我們重用二階中心差分:節點 j 處的 u_xx 約為 (u_{j-1} - 2*u_j + u_{j+1}) / h^2。三個相鄰值,加權 +1、-2、+1,再除以 h^2。這個小小的加權節點花樣,就是一個有限差分模板——一個可重複使用的形狀,我們在每個內部網格點蓋上它,純靠算術就讀出一個導數。
那個「中間是 -2」的花樣從哪來?又有多好?又是泰勒定理。把 u(x + h) 與 u(x - h) 的展開相加:一階導數項因對稱而消掉,u(x) 項合併,除以 h^2 之後存活下來的是 u''(x) + (h^2 / 12)u''''(x) + ...。所以這個模板的截斷誤差是 O(h^2):它是資料的精確*二階導數,加上一個大小正比於 h^2 的殘餘。把 h 減半,那誤差掉為四分之一。這個 O(h^2) 標籤就是模板在空間上的精度階數,也是我們追問一個格式好不好時,第一個要查的數字。
組裝格式——並向前推進
現在把 u_t = a * u_xx 的兩邊都在網格節點 (j, n) 處離散化。對時間導數,用時間上最簡單的前向差分:(j, n) 處的 u_t 約為 (u_j^{n+1} - u_j^n) / k。對空間導數,在當前列 n 上用那個三點模板。令兩個近似相等,得到 (u_j^{n+1} - u_j^n) / k = a * (u_{j-1}^n - 2*u_j^n + u_{j+1}^n) / h^2。因為未來那列上只出現一個未知數 u_j^{n+1},我們可以直接把它解出來。這就是FTCS 格式——前向時間、中心空間——也是顯式方法的原型。
let r = a * k / h^2 (the mesh ratio)
u_j^{n+1} = u_j^n + r * ( u_{j-1}^n - 2*u_j^n + u_{j+1}^n )
row n+1 (new) <- built entirely from row n (old)
stencil shape: u_{j-1}^n u_j^n u_{j+1}^n
\ | /
+----- u_j^{n+1}看看我們造出了什麼。要得到一整列新值,我們讓 j 掃過內部、套用那一條公式——每個未來值都顯式地用我們已知的值寫出來。沒有線性方程組要解,每節點只要一次乘加。用網格比 r = a*k / h^2,更新式讀作 u_j^{n+1} = r*u_{j-1}^n + (1 - 2*r)*u_j^n + r*u_{j+1}^n,是三個鄰居的加權平均。當 0 < r <= 1/2 時,三個權重全為非負,新值就是一個貨真價實的平均,令人安心地貼近物理:熱會被抹平。這是穩定性的第一聲低語,而且絕非巧合——下一篇會把這個暗示講得精確。
這條離散方程,意思對嗎?
我們把連續的偏微分方程換成了一條代數更新式,但我們換到的是同一條方程嗎?誠實的檢驗就是相容性:把真解 u(x, t) 代進離散公式,問它違反這條式子的程度有多糟。那份不吻合就是格式的局部截斷誤差,而它正好是我們已經算過的兩個模板誤差之和——前向時間差分的 O(k),加上中心空間差分的 O(h^2)。當 h 與 k 都縮到零,那誤差消失,所以這條離散方程確實收斂到 u_t = a*u_xx。當一個格式的截斷誤差隨網格間距趨於零時,它就是相容的;FTCS 是相容的,階數為 O(k) + O(h^2)。
你已掌握的,以及這一級往哪走
你現在掌握了有限差分法的完整流程:在空間與時間上鋪一張網格,把每個導數換成由鄰近網格值組成的模板,在每個內部節點令兩邊相等,再從初始條件開始一列接一列地把解往上推進。我們是對熱傳導方程做的,但這份配方是整門學問的模板——波動方程 u_tt = c^2 * u_xx 只是在時間上也用一個二階差分模板,給出蛙跳格式;對流與完整的流體模型也照同一套劇本走,只是換不同的模板。我們算出的每個值當然都是近似:模板的截斷誤差立下一道地板,我們負擔得起的網格決定解析度,而浮點捨入則潛伏在這一切底下。
兩個懸念帶進本級接下來的篇章。第一,FTCS 把空間模板放在舊列上,使更新顯式而便宜——但我們暗示過它只在 r <= 1/2 時安全。改把模板放在新列上,每個節點就與它未知的鄰居耦合,逼你每一步都得解一個線性方程組:這就是隱式格式,也就是第二篇要拆解的顯式對隱式抉擇。第二,那個 r <= 1/2 的限制絕非民間傳說;第四篇會從馮諾伊曼分析與著名的 CFL 條件乾淨地推導出它。夾在中間的第三篇,釘下把這一切綁在一起的深刻結果——對一個相容的格式而言,穩定性恰好就是收斂所需要的全部。網格與模板是容易的部分;讓它們說真話,才是這趟攀登的其餘路程。