從無窮菜單到有限菜單
上一篇把我們的模型問題——比方說在區間 [0, 1] 上、配上 u(0) = u(1) = 0 的 -u''(x) = f(x)——改寫成它的弱形式:找一個函數 u,使得 u'(x)*v'(x) 的積分等於 f(x)v(x) 的積分,對每一個*合格的測試函數 v 都成立。這一行是精確的,但它要的 u 活在一個無窮維的函數空間裡——大到根本存不下。整個有限元素法就是一個有紀律的妥協:別再搜尋那個無窮空間,只搜尋一個我們真能用數字寫下來的有限維子空間。
這就是上一篇的 Galerkin 法,如今具體化了。我們挑一組有限的建構積木函數 phi_1, phi_2, ..., phi_N——也就是基底——並堅持近似解必須是一個組合 u_h(x) = c_1*phi_1(x) + ... + c_N*phi_N(x)。未知數不再是一整個函數,而只是 N 個數,也就是係數 c_1, ..., c_N。接著 Galerkin 說:要求弱形式成立的對象不是每一個測試函數,而只是依次取每一個基底函數 phi_i 當作 v。這就是 N 條方程配 N 個未知數——一個有限的線性系統,可以交給前幾級的求解器。
網格:把區域切成元素
在能建構局部函數之前,我們得先有個地方掛它們。把區間 [0, 1] 在節點 0 = x_0 < x_1 < ... < x_N < x_{N+1} = 1 處切開。每一小段子區間 [x_{j}, x_{j+1}] 就是一個元素,而所有節點與元素的整體,就是有限元素網格。在一維裡,網格只是一串珠子;在二維裡,它通常是鋪滿一塊區域的三角形拼布,在三維裡則是四面體的堆砌。網格不必均勻——這份自由正是本法的超能力。你可以在解很尖銳的地方(裂縫尖端、邊界層)擠進細小的元素,在它平緩的地方讓元素伸展開來,貼合規則網格永遠追不上的彆扭形狀。
現在來談基底。在這張網格上,給每個內部節點 x_j 掛上一個形狀函數 phi_j:一個連續的帳篷(或帽子)函數,它在自己的節點 x_j 處等於 1,朝兩個相鄰節點線性下降到 0,其餘地方則恰好為 0。想像一頂搭在 x_j、釘樁在 x_{j-1} 與 x_{j+1} 的馬戲團帳篷。因為每頂帳篷只在它相鄰的兩個元素上不為零,這些基底函數幾乎不重疊——而這份局部性就是整盤棋的關鍵。用這些帳篷組出來的函數 u_h 自動是一個分段線性內插:節點之間是直線段,而 c_j 就是節點 j 處的高度(也方便地恰好是 u_h(x_j) 的值)。
矩陣從哪裡來
把 u_h = c_1*phi_1 + ... + c_N*phi_N 代入弱形式,讓 v 跑遍每一個 phi_i。u_h'*phi_i' 的積分變成對 j 求和、每一項是 c_j 乘上 phi_j'*phi_i' 的積分。把這些積分收進一個矩陣 A,元素為 A_ij = phi_i'(x)*phi_j'(x) dx 的積分,把右端收進一個向量 b,b_i = f(x)*phi_i(x) dx 的積分。無窮維的弱問題就坍縮成熟悉的有限系統 A x = b,其中 x = (c_1, ..., c_N) 裝著各節點的值。這個矩陣 A 就是剛度矩陣——名字借自它在結構工程裡的出身,在那裡 A_ij 真的衡量節點 i 對節點 j 處推力抵抗得有多硬。
局部帳篷在這裡帶來驚人的回報。元素 A_ij 是 phi_i'*phi_j' 的積分,而這個乘積只在兩頂帳篷都不為零的地方才不為零——也就是只有當節點 i 與 j 同在一個元素裡(即彼此是鄰居)時。所以除非 |i - j| <= 1,否則 A_ij = 0。在一維裡,剛度矩陣是三對角的;在二維裡每個節點只碰到它那少數幾個網格鄰居,給出一個漂亮的稀疏矩陣,每列在 N 個位置中只有寥寥幾個非零。我們從不儲存、也不想要一個稠密的 N×N 陣列。對間距 h 的均勻網格,簡短計算給出對角線 2/h、每個次對角線 -1/h——再整體除以 h,你或許會注意到它正好就是上一級那個 (-1, 2, -1) 有限差分模板。兩種方法在這個簡單問題上相遇了;一般情形它們不同,但這份親緣是真的。
weak form, with u_h = sum_j c_j phi_j and v = phi_i :
A_ij = integral_0^1 phi_i'(x) phi_j'(x) dx
b_i = integral_0^1 f(x) phi_i(x) dx
A x = b , x = (c_1, ..., c_N)
uniform mesh, spacing h, hat functions:
[ 2 -1 ] tridiagonal,
A = 1 [ -1 2 -1 ] symmetric,
--- [ -1 2 -1 ] positive definite
h [ . . .] (sparse: |i-j|<=1)組裝,以及矩陣的好教養
實務上沒人用一個龐大的全域積分去算 A_ij。我們改採逐元素組裝矩陣。對每個元素跑迴圈,只在那個元素上活動的形狀函數之間,算一個小小的局部 2×2(一維)或 3×3(三角形)元素矩陣的重疊積分,再把這些小數字加進全域 A 的對應格子裡。每個元素貢獻一點,貢獻在元素共用節點處堆疊起來。這個逐元素的局部迴圈是每一份真實有限元素程式的引擎:它能一致地處理非結構化網格與彎曲元素,也正是等參元素映射(把參考三角形扭曲到彎曲三角形上)施展身手之處。
由對稱弱形式得到的剛度矩陣是對稱的(A_ij = A_ji,因為被積函數 phi_i'*phi_j' 不在乎順序)且正定的(它源自一個確實被極小化的能量)。這是份禮物:一個對稱正定系統正是線性代數那幾級裡最舒服的情形——直接解可用Cholesky 分解,大型稀疏的可用共軛梯度法,都不必為主元換位傷腦筋。但仍有兩點得誠實交代。細化網格會讓 A 變大也變得更病態——它的條件數約以 O(1/h^2) 增長,所以更細的網格悄悄吃掉你一些雙精度位數,並拖慢迭代求解器。而且在施加邊界條件之前,赤裸的 A 是奇異的,這正是下一段的主題。
邊界條件、質量矩陣,以及我們的落腳處
並非所有邊界條件都以同樣方式進入——這是個確實微妙、值得釘清的點。一個固定 u 值的條件,例如 u(0) = 0,是本質邊界條件:它無法用弱形式的積分表達,所以我們直接強加它,把那個節點的未知數移除(它不自由,係數就是 0)。一個關於通量的條件,例如 u'(1) = g,是自然邊界條件:它從建構弱形式的分部積分裡掉出來,就直接出現在右端 b——你什麼特別的事都不必做,它自己會處理好。這個本質對自然的區分常令初學者吃驚,而把它搞反是有限元素法的經典臭蟲。
有一個姊妹矩陣值得一個名字。若問題還帶著一個未經微分的 u 項——一個含時的熱傳導方程 u_t = u_xx,或一個特徵值問題——同樣的組裝會產生第二個矩陣 M,元素為 M_ij = phi_i(x)*phi_j(x) dx 的積分(是形狀函數本身,不是它們的導數)。這就是質量矩陣,同樣稀疏且對稱正定,它正是把空間上的有限元素法與常微分方程那幾級的時間步進方法接起來的東西:用元素離散空間,你就剩下一個關於節點係數的常微分方程組 M*c'(t) = -A*c(t)——又是線法,如今時間導數上多了一個有品味的質量矩陣。
退一步,看看你已跨過的這座橋。上一篇抽象的弱形式,已變成一個具體、稀疏、對稱正定的系統 A x = b,其未知數是一個分段線性解的各節點值。它裡頭的每個數仍是近似:我們把搜尋限制在帳篷函數上(這個建模選擇,誤差理論與 Cea 引理會告訴我們,在恰當的範數下代價為 O(h^2)),積分帶著求積誤差,而條件數 O(1/h^2) 在我們細化時侵蝕掉幾個位數。下一篇真正去求解這個系統並讀出它的答案;再下一篇則衡量它的準確度,並在痛處細化網格。把微積分變成矩陣這一道困難的概念跳躍,如今已在你身後。