數值微分

 

2. 太小真的不好

先前的討論我們處理了控制 truncation error 的問題,無論是在哪一個階數的 truncation error 乍看之下我們希望讓 \(h\) 越小越好,然而當我們考慮到round-off error 時情況就不是如此。

假設我們要計算 \(f'(x)\), 我們可能會利用到

\(f(x+2h), f(x+h), f(x), f(x-h), f(x-2h)\)

這幾個點。而這些理想中的數字到了電腦裡會有 round-off error ,實際上會存在電腦裡的值會是

\(y_2 = f(x+2h)+e_2, y_1 = f(x+h)+e_1, y_0 = f(x)+e_0, y_{-1} = f(x-h)+e_{-1}, y_{-2} = f(x-2h)+e_{-2}\)

我們以 \(D_{f1}(x,h)\) 為例,

\(D_{f1}(x,h) = \frac{y_1-y_0}{h} = \frac{f(x+h)+e_1-f(x)-e_0}{h} = f'(x)+\frac{K_1}{2}h +\frac{e_1-e_0}{h}\)

這裡的數值誤差就是 \(|D_{f1}(x,h)-f'(x)| = \frac{|K_1|}{2}h +\frac{|e_1-e_0|}{h}\)。前面 \(K_1 = f^{(2)}(\xi)\) 的部份是 truncation error 貢獻的,後面 \(e_i\) 的部份是 round-off error 貢獻的。很明顯這誤差是一個和 \(h\) 有關的函數,並且有一項是 \(\frac{1}{h}\) ,這讓 \(h\) 靠近0 時會使得 round-off error 放大,但也會讓 truncation error 變小,這稱為 step size dilemma。雖然說不准 \(e_i\) 到底是多少,但配合不同位元的電腦會有一個對應的上限。以32 位元的電腦為例,雙精度的數值及算 \(e_i\) 的值會小於 \(10^{-16}/2\),我們令這個值為 \(\epsilon\)。因此,我們計算 \(|D_{f1}(x,h)-f'(x)| \leq \frac{2\epsilon}{h} + \frac{|K_1|}{2}h\) 的最小值發生在,\(h_0 = 2\sqrt{\epsilon/|K_1|}\)。

這樣的分析還是讓人無法決定 \(h\) 的最適大小,因為我們還不知道 \(K_1\) 的影響到底有多大?當我們限制自己不去對函數 \(f (x)\) 進行解析的情況下,對於任意電腦可以表達的 \(h\) 與整數 \(k\) ,我們有 \(f(x+kh)\) 的值,我們也有 \(D_{f1}(x,h)\) 的值。令 \(h_i = 10^{-i}\) 時,我們還可以觀察 \(D_{f1}(x,h_i)\) 的行為,特別是當我們觀察 \(D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})\) 時,該式可以整理成

\(D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})= 2\epsilon (9)10^{i-1} - 9K10^{-i}/2\)

我們知道 \(\epsilon\) 大約是多少 \(10^{-16}/2\),我們也可以調整 \(i\) 使得 \(D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})\) 最靠近零,意思是當我們找到一個 \(i'\) 使得 \(D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})\) 最靠近零時,我們可以預估 \(K\) 的值是多少?例如,當 \(i' = 9\) 時 \(|D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})|\) 有最小值,我們就可以預估 \(K\) 的值大約是 \(2(10^{2i'-17-i})=20\)。我們可以把這個 \(K\) 值當作是 \(K_1\) 的估計值,因此我們再看一下對應的 \(|D_{f1}(x,h_i) - f'(x)|\) 大約是 \(10^{i-16}+10^{2i'-17-i}\),因此 \(h_i\) 在 \((10^{-9},10^{-8})\) 之間時會使得我們的差分值最靠近微分值。從以上關係,我們得到一個重要的結論,當最適合的 \(i\) 會落在 \((i'-1,i')\) 這個區間。這樣我們就可以利用觀察 \(|D_{f1}(x,h_i)-D_{f1}(x,h_{i-1})|\),讓我們在沒有真解的情況下決定最適合的 \(h_i\)。可以看見在不考慮\(K_1\) 的影響下,最好的 \(h\) 大約是 \(10^{-8}\) 等級左右。