從連續的積分到一個大矩陣
到現在為止，概念上最難的工作你都已經做完了。前幾篇指南把一條 PDE 重寫成弱形式，挑了一張由小元素組成的網格，並選了形狀函數——一個個帳篷形狀的帽子函數，每個節點一個——作為近似解的基底。接著Galerkin 法說：把未知的 u 寫成「未知係數乘上那些帽子」的和，要求弱形式方程對每個帽子作為測試函數都成立，於是掉出一個有限的線性方程組 A u = b。本篇要談的，就是那個不起眼卻絕對核心的機器，它把那份承諾變成真正的數字：定義 A 與 b 的積分如何被算出來、塞進矩陣裡，然後你又如何解它。
矩陣 A 就是剛度矩陣：元素 A_{ij} 是「帽子 i 的梯度點乘帽子 j 的梯度」在整個區域上的積分（以 -u'' = f 這類問題為例）。向量 b 是負載項：b_i 是源項 f 乘上帽子 i 的積分。建造 A 的天真做法，是對所有節點對 (i, j) 跑一個雙重迴圈，每對都算一個全域積分——但這既慢又浪費，因為幾乎每一對都給出零。兩個帳篷不重疊的帽子毫無貢獻。漂亮的修補是把迴圈裡外翻轉：不要對節點對跑迴圈再去拜訪元素，而是對元素跑迴圈，再去拜訪每個元素碰到的那少數幾個節點。這個重新組織，就是組裝的全部精神。
局部元素矩陣
把鏡頭拉近到單一個元素——比方說一維裡的某個小區間 [x_k, x_{k+1}]，或二維裡的一個三角形。只有以「這個元素自己的節點」為中心的帽子，在它內部才非零：一段線段有兩個節點，一個線性三角形有三個。所以這個元素對那個巨大剛度矩陣的貢獻，可以用一個微小的元素矩陣來捕捉——一維裡是個 2×2 的方塊，三角形則是個 3×3 的方塊。你只在這一個元素上對「梯度點乘梯度」積分，就能算出它的各個元素；在這裡帽子只是簡單的直線片段，積分很容易。對一張間距為 h 的均勻一維網格，每個內部元素都給出同樣的小剛度方塊 (1/h) 乘以 [[1, -1], [-1, 1]]——這個結果你動手一分鐘就能推出來。
在真正的程式裡，一般元素上的積分並不是用封閉形式算的；它是用數值積分法——通常是高斯積分——在少數幾個點上求值，而且是在一個固定的參考元素（單位區間、單位三角形）上做，再映射到實際的形狀。那個從參考到實體的映射，正是等參元素所形式化的東西，也是用同一個小核心去處理彎曲邊界與扭曲格子的方法。要抓住的重點是：建造局部矩陣既便宜、又是局部的，而且對給定類型的每個元素，結構都一模一樣。網格只是一次遞給你一個元素而已。
灑加：組裝迴圈
現在是核心訣竅。每個元素是用局部編號來認得自己的節點的（這段線段的節點 1 與節點 2），但這些節點在完整網格裡有全域編號（比方說節點 7 與節點 8）。組裝就是把局部元素矩陣的每個元素，依照全域節點編號所指的位置，加進全域剛度矩陣裡的動作。是「加」，不是「寫」：一個內部節點屬於兩個元素（或好幾個三角形），而每一個都貢獻那個節點對角線元素的一部分。這些貢獻會累積起來。正是這個「灑開再相加」，使得均勻一維網格的全域 A 對角線最後是 2/h 而非 1/h——相鄰的兩個元素各自在那裡存入 1/h。
A = zeros(N, N); b = zeros(N) // global, all N nodes
for each element e in the mesh:
Ke = local_stiffness(e) // small dense block, e.g. 2x2 or 3x3
fe = local_load(e) // small local load vector
g = global_node_ids(e) // map: local index -> global index
for local a in nodes(e):
b[g[a]] += fe[a] // scatter-add the load
for local c in nodes(e):
A[g[a], g[c]] += Ke[a, c] // scatter-add the stiffness
// A is now the assembled global stiffness matrix; b the global load有兩件事讓這個迴圈成為瑰寶。第一，它的代價是 O(元素數)，而非 O(N^2)——你只碰每個元素一次，做固定且極少量的工作。第二，組裝出來的 A 是稀疏的：每個節點只和它的網格鄰居耦合，所以無論網格長到多大，一列都只有寥寥幾個非零元素。你絕不該把 A 存成一個滿的 N×N 陣列；你要用稀疏格式來建它（一串你最後加總起來的 (列, 行, 值) 三元組，或一個壓縮列的結構），好讓百萬等級的 N 還塞得進記憶體。稀疏性在這裡不是個小小的最佳化——它是「一個解得動的模型」與「一個需要一兆位元組記憶體的模型」之間的分界。
釘住邊界
如果盲目地對所有節點組裝，矩陣 A 其實是奇異的——它沒有唯一解，因為到目前為止還沒有任何東西指明解被錨定在哪裡。一個純粹的 Neumann（通量）問題，會讓解可以整體上下平移一個常數；那份自由表現為一個零特徵值。補救之道是施加邊界條件，而弱形式早已把它們誠實地分成兩類。自然邊界條件（指定通量，Neumann）不需要特別處理——它們本就烘焙在弱形式的邊界項裡，直接進入 b。本質邊界條件（指定值，Dirichlet）則必須直接施加在方程組上，因為它們釘死的是未知量本身。
施加 Dirichlet 值 u_j = g 最乾淨的方式是「抬升並消去」法：把那個已知值的影響搬到右邊（對其他每一列，從 b 減去 g 乘上第 j 行），然後整列整行第 j 刪掉，把方程組縮小到真正未知的那些節點。一個比較偷懶但常見的捷徑，是把第 j 列覆寫成「對角線放 1、b 裡放 g」——寫起來快，雖然它會破壞 A 的對稱性，而下一節會告訴你你想保住那份對稱性。無論哪種做法，本質條件一旦施加，方程組就變成非奇異，答案也被釘住了。如果你是在時間上推進一條熱傳導或波動方程，這裡也正是你把剛度矩陣在時間上的搭檔——質量矩陣——併進來的地方。
求解 A u = b，別丟掉它的天賦
現在你手上握著一個真正的線性方程組 A u = b——而你在前幾級學過的直接法與迭代法求解器，全都回家了。新手的第一直覺，計算 u = A 的反矩陣乘 b，正是不該做的事：求反矩陣既昂貴又稠密（稀疏矩陣的反矩陣通常是滿的），而且在數值上比直接解方程組更糟。有限元素的 A 通常是對稱且正定的，所以最自然的直接求解器是 Cholesky 分解，A = L L^T——只需一般 LU 一半的工作量，而且不必選主元就保證穩定。再搭配一個用重排序來讓因子保持稀疏的稀疏直接求解器，這就是中等規模問題的首選。
對於三維模擬裡真正龐大的網格，連稀疏的 Cholesky 因子也填入得太多，於是你改用迭代求解器。由於 A 是對稱正定的，完美的搭配是共軛梯度法——它只需要「矩陣乘向量」的乘積，而這在稀疏的有限元素矩陣上代價是 O(非零元素數)，並且它完全不形成任何因子。共軛梯度法的收斂速度由 A 的條件數主宰，而這裡有個誠實的陷阱：把網格細化、縮小 h，並不只是多加未知量，還會讓 A 的條件變差（對二階問題大致像 1/h^2）。所以一張更細、更準確的網格，同時也是一個更難解的方程組——這正是為什麼一個好的前置條件子，例如多重網格，在大規模下不是可有可無，而是不可或缺。