「剩下的誤差」這個想法
在上一篇我們認識了中心差分 D(h) = (f(x+h) - f(x-h)) / (2h)，把它當作 f'(x) 的一個近似，也看到了它的截斷誤差：泰勒展開告訴我們 D(h) = f'(x) + c2*h^2 + c4*h^4 + ...，其中 c2 = f'''(x)/6，更高次的項則騎在 f''''' 、 f''''''' 等等之上。重點標題是誤差為 O(h^2)：把 h 減半，誤差縮成四分之一。但仔細看它的結構。這個誤差並不是一團沒有形狀的模糊——它是一個關於 h 的已知多項式，由單一的領頭項 c2*h^2 主宰，那個係數我們也許不知道，但它的形狀我們是知道的。
這單單一個事實，就是整齣戲的開場。如果領頭誤差的行為恰好像 h^2，那我就可以刻意算出兩個估計，一個用步長 h、一個用步長 h/2，再讓它們互相較勁。兩者都錯，但它們錯得可預測、成比例——而兩個你知道其比值的錯，可以組合成一個對。這就是理查森外推法的核心：不要跟領頭誤差硬碰，把它量出來，然後減掉它。同樣的想法適用於任何「準確度的階是 h 的已知次方」的方法，而那幾乎就是整個數值分析。
兩個錯湊成一個對：那次消掉
用 A 表示我們想要的真值，A(h) 表示步長為 h 時的估計，並假設 A(h) = A + c*h^p + （更高次項），其中階 p 是已知的——對中心差分而言 p = 2。把步長減半得到 A(h/2) = A + c*(h/2)^p + ... = A + c*h^p / 2^p + ... 。如今 A(h) 的領頭誤差，恰好是 A(h/2) 領頭誤差的 2^p 倍。於是組出一個能殺掉那一項的組合：取 2^p 份較細的估計，減掉一份較粗的，再除以 (2^p - 1)。未知的 c 乾乾淨淨地消掉，A 存活下來，而誤差剩下的部分是下一項——對於像中心差分這種「只含偶數次方」的級數，現在是 O(h^{p+2})。
Improved = ( 2^p * A(h/2) - A(h) ) / ( 2^p - 1 ) # central difference: error is O(h^2), so p = 2, 2^p = 4 Improved = ( 4 * D(h/2) - D(h) ) / 3 # error now O(h^4)
舉個具體例子。假設 D(h) = 1.040、D(h/2) = 1.012，兩者都是對某個真值為 1.000 的導數所做的粗略猜測。理查森外推給出 (4*1.012 - 1.040)/3 = (4.048 - 1.040)/3 = 3.008/3 = 1.0027——已經比任何一個輸入都接近得多，而我們花的也不過就是那兩個估計本來就需要的函數求值。兩個誤差 +0.040 與 +0.012，大致呈理論所預測的 4 比 1 之比，而那個加權差把當中大部分的誤差擠了出去。
一次一階，往上爬那張表
為什麼只走一步就停？改良後的估計有它自己的領頭誤差，現在是 p+2 階，帶著它自己的未知係數——和之前是同一種局面，只是高了一階。那就再外推一次，這回用為新的階所選的權重。建一張三角形的表：第一欄放原始估計 A(h)、A(h/2)、A(h/4)、...；之後每一欄，都由左邊那一欄的兩個相鄰元素以理查森組合而成，用的權重是 2^q，q 是那一欄所要殺掉的階。每多一欄就抹掉下一個誤差項，於是準確度能以 O(h^2)、O(h^4)、O(h^6)、... 的速度快得驚人地往上爬——對光滑函數而言，減半個幾次就能逼近機器精度。
龍貝格法：理查森與梯形法則的聯姻
最有名的回報落在積分上。複合梯形法則是那匹謙遜的主力馬——用一個個直線梯形去近似曲線下的面積——它本身只有 O(h^2) 的準確度，聽起來平庸。但對一個光滑的被積函數，一個深刻的結果（歐拉-麥克勞林公式）說：它的誤差是一個關於 h 的乾淨偶數冪級數：誤差 = c2*h^2 + c4*h^4 + c6*h^6 + ... 。這恰恰就是理查森所狼吞虎嚥的那種結構。於是在 h、h/2、h/4、... 處算出梯形估計，再讓它們沿著外推表往上跑。結果就是龍貝格積分，而梯形法則的第一次外推，竟然恰好就是喬裝過的辛普森法則——一個令人滿意的暗示：我們重新發現了下一篇會正面處理的東西。
- 用一個區間算出複合梯形估計 R(0,0)，然後反覆把步長減半，得到 R(1,0)、R(2,0)、...。重複利用先前的求值——每次減半只新增那些新的中點，所以加密很便宜。
- 用 R(i,j) = ( 4^j * R(i, j-1) - R(i-1, j-1) ) / ( 4^j - 1 ) 填入每一欄外推值。這裡的 4^j = 2^p 扮演那個權重，殺掉第 2j 階的誤差項。
- 對角線上的 R(i,i) 就是你最好的估計；盯著相鄰的對角元在你在乎的位數上不再變動，就在那裡停手。
實務上的魔法在於：每次梯形加密都重複利用較粗網格上的每一個函數值，所以主要成本只是那些新的點，而建立外推欄不過是幾次便宜的算術組合。你大致只付了最細梯形網格的代價，卻得到任何可行大小的梯形網格都達不到的準確度。對於一個有限區間上的光滑被積函數，龍貝格法常常從一張只有五、六列高的表裡，就釘下十位以上正確的數字。
它救不了你的那些事
理查森外推攻擊的是截斷誤差——公式本身的誤差——而且只攻擊它。它對捨入誤差無能為力，而對微分來說這個區別格外尖銳。回想上一篇的步長取捨：當 h 縮小，截斷誤差下降，但 (f(x+h) - f(x-h)) 裡的相減抵消使捨入誤差增長，於是存在一個最佳 h，比它更小的步長反而把事情弄糟。用非常小的 h 去外推一個有限差分公式並不能逃離那道地板；它甚至可能放大它，因為權重 4、-1 會把已經污染輸入的雜訊一同放大。在外推一個導數時，請用適中、而非微小的步長。
這張表能爬多高，也有它誠實的極限。外推法假設存在無窮多個光滑的誤差項，但在浮點數裡，每一欄都是由「本身只準到約 16 位」的輸入算出的，所以一旦某欄達到那個準確度，之後的欄就只是在搬弄捨入雜訊——表停滯了，而更後面的對角元甚至可能漂離真值。讀龍貝格表或微分表的正確方式，是盯著對角線：相信「相鄰元首次在你所需精度上一致」的那個點，並忽略它之後的翻攪。這個方法是一台漂亮的加速器，不是一台永動機。