問題:穿過每一點的曲線
假設一個氣象站每小時記錄一次溫度,而你需要對下午 2:37——一個沒人量測過的時刻——給出合理的估計。又或者一張物理常數表給了某材料在 0、100、200 度時的密度,而你想要 150 度時的值。這個反覆出現的任務是同一個:有人交給你幾個點 (x_0, y_0)、(x_1, y_1)、...、(x_n, y_n),要你找一條「恰好」穿過它們全部的平滑曲線,好讓你讀出中間的值。這就是插值,而最簡單可用的曲線是多項式——形如 a_0 + a_1 x + a_2 x^2 + ... 的冪次和。多項式無比方便:計算、微分、積分都易如反掌,這正是為何多項式插值坐落在如此多數值工作的基石上——數值積分與常微分方程求解器都悄悄建立在它之上。
在我們去獵捕曲線之前,有一個事實讓整個主題變得堅實:它是「唯一」的。給定 n + 1 個 x 值「互異」的點,恰好存在唯一一個次數至多為 n 的多項式穿過它們全部。兩點決定一條直線(一次),三點決定一條拋物線(二次),四點決定一條三次曲線,依此類推——n + 1 個條件釘住一個 n 次多項式的 n + 1 個係數。這不是好幾個彼此競爭、任你挑選的答案;而是同一個物件,我們即將用三種不同的筆跡把它寫下來。附帶條款裡那句「次數『至多』為 n」很要緊:三個共線的點還你一條直線,而非硬湊出一條拋物線。
誠實卻凶險的做法:直接解係數
最直白的想法是直接找出係數。寫下 p(x) = a_0 + a_1 x + ... + a_n x^n,並在 n + 1 個資料點上各要求 p(x_i) = y_i。這給出 n + 1 條關於 n + 1 個未知係數的線性方程式——一個方程組 A a = y,你可以用本階梯前面學過的高斯消去法來解。這裡的矩陣,其第 i 列為 (1, x_i, x_i^2, ..., x_i^n),就是著名的范德蒙矩陣。它的行列式是所有 i < j 的 (x_j - x_i) 之乘積,恰在節點互異時非零——所以這個觀點為我們剛宣稱的唯一性給出一個乾淨的線性代數「證明」。
觀念上可愛,計算上可怕。范德蒙矩陣隨次數增長變得嚴重「病態」:對等距節點,其條件數大致以 2^n 爆炸。回想條件數那一級的經驗法則——一個接近 10^8 的條件數,已經要花掉你大約 16 個雙精度數字裡的 8 個。在等距格點上到了約 20 次,求得的係數就成了數值垃圾,儘管真正的插值多項式完全沒問題。這就是「問題」(插值,本身乖巧)與解它的一個「壞演算法」(建立並求解 A a = y)之間的差別。所以我們把范德蒙矩陣留作解釋,而伸手去拿更好的食譜來真正計算——這正是「寧取後向穩定的方法、勿用樸素方法」的精神。
拉格朗日的食譜:直接調配 y 值
這裡有一種寫下答案、完全不必解任何方程組的方式。拉格朗日插值把多項式建構成 y 值的加權混合,而每個權重都是一個量身訂做的小多項式。對第 i 個節點,你打造一個基底多項式 L_i(x),讓它在 x_i 處等於 1、在其他每一個節點處等於 0。建構很直接:把所有 j 不等於 i 的因式 (x - x_j) 相乘(這讓它在每個其他節點為零),再除以同一乘積在 x_i 代入後的值(這把它在 x_i 正規化為 1)。於是整個插值多項式就是 p(x) = y_0 L_0(x) + y_1 L_1(x) + ... + y_n L_n(x)。
它為何穿過每一個點?代入任意節點 x_k。每個 i 不等於 k 的基底多項式 L_i 都被造成在 x_k 處消失,所以那些項全都不見了;只有 L_k 存活,且它等於 1,留下 p(x_k) = y_k。它依建構自動命中每一個資料點——沒有方程式、沒有消去、沒有范德蒙矩陣。兩點時這就重現穿過它們的直線;三點時則是拋物線。它與先前是同一個唯一多項式,只是用「把工作分攤到各 y 值上」的筆跡寫出來。
Through (0,1), (1,3), (2,2):
(x-1)(x-2) (x-0)(x-2) (x-0)(x-1)
L_0 = -------------- L_1 = -------------- L_2 = --------------
(0-1)(0-2) (1-0)(1-2) (2-0)(2-1)
p(x) = 1*L_0(x) + 3*L_1(x) + 2*L_2(x)
Check x = 1: L_0(1)=0, L_1(1)=1, L_2(1)=0 -> p(1) = 3拉格朗日形式在理論與證明上很優美——答案被明確寫出,看不到一絲線性代數。它誠實的缺點是實務上的。在這個樸素形式下,在一個查詢點計算 p(x) 要花 O(n^2) 的工作量,更糟的是,「加入」一個新資料點會迫使你從頭重建「每一個」基底多項式,因為所有 (x - x_j) 乘積都變了。牛頓形式修正了第二個抱怨,而一個叫做重心公式的巧妙改寫則一次修正了兩者——但那個打磨過、又快又穩定的版本,是本級後面切比雪夫那篇的主題。
牛頓的形式:一次一個點地長出曲線
現在想像你要逐步地建構曲線,讓每個新資料點只「加上」一個修正,而不擾動你已有的擬合——就像為級數附加一個更細的項。牛頓差商形式正是如此。它把「同一個」唯一多項式寫成一個巢狀和:p(x) = c_0 + c_1 (x - x_0) + c_2 (x - x_0)(x - x_1) + c_3 (x - x_0)(x - x_1)(x - x_2) + ...。看這個結構:每個新因式 (x - x_j) 在所有「較早」的節點處都為零,所以附加一個新項不可能破壞較早項已釘好的點。這正是為何加入一個資料點只需附加一個係數與一個因式——很便宜,而拉格朗日卻要求整個重建。
係數 c_k 就是差商,在一張小小的三角形表中計算。零階差商就是 y 值:f[x_i] = y_i。一階差商是相鄰點之間的斜率:f[x_i, x_{i+1}] = (f[x_{i+1}] - f[x_i]) / (x_{i+1} - x_i)。更高階遞迴而得——每一個都是兩個較低階差商之差,除以它們所涵蓋的 x 值跨度:f[x_i, ..., x_{i+k}] = (f[x_{i+1}, ..., x_{i+k}] - f[x_i, ..., x_{i+k-1}]) / (x_{i+k} - x_i)。這張表的領先邊緣——f[x_0]、f[x_0,x_1]、f[x_0,x_1,x_2]、...——直接讀作 c_0、c_1、c_2、...。
- 把 y 值寫進表的第一欄:f[x_0] = y_0、f[x_1] = y_1、...、f[x_n] = y_n。
- 由每一對相鄰項形成一階差商:f[x_i, x_{i+1}] = (f[x_{i+1}] - f[x_i]) / (x_{i+1} - x_i)。
- 向更高階遞迴:每個較高階差商是(右下那個減左下那個)除以 (x_{i+k} - x_i),填出越來越短的欄。
- 從表的頂部對角線上直接讀出係數 c_0、c_1、...、c_n。
- 用巢狀(霍納式)乘法計算 p(x),從最深的項向外運算——每個查詢點 O(n) 的工作量。
用同樣三個點 (0, 1)、(1, 3)、(2, 2) 做個小小的計算:y 值給出 f[0] = 1、f[1] = 3、f[2] = 2。一階差商:f[0,1] = (3-1)/1 = 2、f[1,2] = (2-3)/1 = -1。二階差商:f[0,1,2] = (-1 - 2)/(2 - 0) = -1.5。所以 p(x) = 1 + 2(x - 0) - 1.5(x - 0)(x - 1)。把它展開,你會得到與拉格朗日形式所產生「完全相同」的那條拋物線,也與范德蒙求解會給出的同一條。一條曲線,三種筆跡。建表是一次性的 O(n^2) 成本;之後每次計算都是 O(n)。
同一條曲線,三種筆跡——以及一個誠實的警告
這值得重申,因為它常令人意外:范德蒙、拉格朗日與牛頓「不」給出三個不同的插值結果。依唯一性定理它們不可能——穿過 n + 1 個互異點、次數至多 n 的多項式只有一個,所以三者都必須產生它。不同的只有記帳方式:設定的成本、計算的成本、數值穩定性,以及加入新點的優雅程度。范德蒙是最清楚的解釋卻是最糟的計算;拉格朗日在理論上最乾淨;牛頓則是資料逐步到來時、或你想要下一篇所倚靠的誤差估計時的實用選擇。
有一座更深的橋值得留意。領先差商與一個導數相連:對資料範圍內某點 xi,有 f[x_0, ..., x_n] = f^(n)(xi) / n!。這並非巧合——它是插值誤差公式的種子,其中下一個被略去的差商,表現得恰如多項式與真實函數之間的差距。換句話說,牛頓形式在它的下一個係數裡就攜帶著自己的誤差估計,這也是它在自適應方案中仍是主力的部分原因。