數值微分

1. Richardson's extrapolation

如果我們把剛才的差分式中的 \(h\) 換成換成 \(2h\) 得到

\(D_{f1}(x,2h) = f'(x) + \frac{2h}{2}f^{(2)}(x) + \frac{4h^2}{3!}f^{(3)}(x)+\dots = f'(x) + \dots\)

這裡有一個有趣的結果就是我們可以利用 \( D_{f1}(x,h) \) 和 \( D_{f1}(x,2h) \)  的線性組合得到 \( f'(x) \) 但是卻有更好的誤差。

\(2D_{f1}(x,h)-D_{f1}(x,2h)=f'(x)-\frac{2h^2}{3!}f^{(3)}(x)+\dots=f'(x)+O(h^2)\)

為甚麼說這是更好的誤差?我們假設 \(h\) 是小於1 且很靠近0 的數字,在數值微分的領域,因為我們無法知道 \(f^{(2)}\) 和 \(f^{(3)}\) 到底等於多少,我們只能假設他們具有一定程度的連續性。因此在某一個區間內,他們會靠近一個常數。接著 \(h^2\) 當然比 \(h\) 更靠近0 。因此,我們期待 \(O(h^2)\) 的誤差會比 \(O(h)\) 來的更小。

這種由forward method 做出來的高階逼近,我們記作 \(D_{fn}(x,h)\) ,下標的\(f\) 代表forward, \(n\) 代表誤差項的次數。同理,我們可以定義 \(D_{bn}(x,h)\),其中 \(b\) 代表backward

\(D_{b1}(x,h) = \frac{f(x)-f(x-h)}{h}\)

若我們把 \(f(x+h)\) 和 \(f(x-h)\) 對 \(f(x)\) 的泰勒展示同時寫開,我們得到

\(D_{c2}(x,h) = \frac{f(x+h)-f(x-h)}{2h} = f'(x) + \frac{h^2}{3!}f^{(3)}(x) + \dots = f'(x) + O(h^2) \)

這裡雖然只用到一次的差分,但是用到的資訊分別是 \(x\) 前一個 \((x+h)\) 和後一個 \((x-h)\),因此我們稱為 central difference,並且用下標符號c 代表。

既然我們可以把一階的誤差項提高為二階,我們當然可以繼續仿造相同的方法,把階數繼續提高。我們有Richardson's extrapolation 公式:

  • \(D_{f,n+1}(x,h) = \frac{2^nD_{f,n}(x,h)-D_{f,n}(x,2h)}{2^n-1}\)
  • \(D_{b,n+1}(x,h) = \frac{2^nD_{b,n}(x,h)-D_{b,n}(x,2h)}{2^n-1}\)
  • \(D_{c,2(n+1)}(x,h) = \frac{2^{2n}D_{c,2n}(x,h)-D_{c,2n}(x,2h)}{2^{2n}-1}\)

這種透過調整間距以及適當的線性組合達到提高逼近能力的方式並不是只出現在數值微分,若我們用數值方法 \(N_p(f,h)\) 去逼近一個數學上關於某個平滑函數 \(f (x)\)的計算,我們稱實際值為 \(M(f)\),並且這一個數值方法的 truncation error 階數是 \(p\),亦即

(式一):\(M(f) - N_p(f,h) = c h^p + d h^q + O(h^r)\)

其中常數 \(c\),\(d\) 與 \(f (x)\) 的高階導數有關,\(p<q<r\)。我們利用兩倍的間距 \(2h\) 來計算數值解,得到

(式二):\(M(f) - N_p(f,2h) = c (2h)^p + d (2h)^q + O((2h)^r)\)

我們的目標是計算 \(M(f)\),並且讓 truncation error term 的階數提昇。我們可以將式一乘上 \(2^p\) 然後減去式二,便可以將 \(h^p\) 次方項消去。新的逼近式為

\(\frac{2^pN_p(f,h)-N_p(f,2h)}{2^p-1} = M(f) + O(h^q)\)

這就是上述公式的由來。

看起來我們可以利用這樣的公式不斷地把階數提高,那麼階數是否越高越好呢?請別忘記,當階數提高時,我們需要用到的差分資訊也越多,並且關於誤差的分析都建築在一個假設下,就是這些點很靠近要計算的點,並且在這個點附近函數是平滑的。這裡有太多的未知,因此考量到計算成本與實際需求的情況下,我們並不會用很高的階數來計算數值微分。