用電腦算微分

前情提要: 用電腦算極限

這裡我們要講的是用數值計算來算函數的微分值.

已知一個函數 $f(x)$ 在某個點 $a$ 的微分值定義是 $$ f'(a) = \lim_{h\to 0} \frac{f(a+h)-f(a)}{h}. $$

我們用一個簡單的例子試試看. 假設我們想求 $f(x)=x^2$ 在 $x=\pi$ 的微分. 根據定義我們有

$$ f'(\pi) = \lim_{h\to 0} \frac{(\pi+h)^2-\pi^2}{h}. $$

接著我們將 $h$ 取靠近 $0$ 的 $11$ 的點並帶入上列這個式子試著來算其極限值. Julia 程式如下

h = range(0, 0.01, length=11); h = h[11: -1 : 2];
fp = ((pi.+h).^2 .- pi^2)./h;
hcat(h, fp, fp.-2*pi)

結果如下:

上列數字最左邊是 $h$ 值, 中間為估計的微分值. 我們發現的確這個數字會越來越接近真實的解, 也就是 $2\pi$, 約等於 $6.283185307179586$. 最右邊為這個估算值與真實值 $2\pi$ 之間的差. 的確, 當 $h$ 越接近零, 這個估計出來的微分值離 $2\pi$ 的距離越來越小.

Question: 用數值計算微分能有多精確? 這個誤差能不能一直遞減下去?

接著我們取更多靠近 $0$ 的點來計算微分的極限值, 我們列出其與真實值 $2\pi$ 之間的差, 並且把它畫出來.

h = range(1, 16, length=16); h = 10.0.^(-h[1: 16]);
fp = ((pi.+h).^2 .- pi^2)./h;
hcat(h, fp, abs.(fp.-2*pi))

結果如下:

觀察最後一下發現, 當 $h$ 很小時微分竟然與真實質差更多, 更不準了!! 比如說當 $h=10^{-14}$ 時, 誤差竟然大到約是 $10^{-1}$.

我們把它畫出來看看, Julia code 如下:

using Plots
plot(log10.(h), log10.(abs.(fp.-2*pi)),label="Error")
xlabel!("log10(h)")
ylabel!("log10(Error)")

發現雖然誤差在 $h$ 大的時候遞減, 不過當 $h$ 接近零的時候卻又遞增上去了.

那誤差最小值出現在什麼時候呢?

我們發現當 當 $h=1.0*10^{-8}$ 時, 其估計出來的微分值離真實值誤差最小, 其誤差為 $4.26 *10^{-8}$.

不過, WHY?? 為什麼誤差不會一直往下遞減? 其實這也是因為 捨入誤差(rounding-error) 的關係.

觀察一下我們的式子 $$ \frac{f(a+h)-f(a)}{h} $$ 當我們在用數值計算這個式子的時候其實並不完全是這樣子, 在分子應該會有捨入誤差在, 也就是說, 其實我們看到的數字應該是以下這個式子算出來的 $$ \frac{f(a+h)-f(a) + \epsilon}{h} $$ 其中的 $\epsilon$ 就是捨入誤差. 所以, 我們計算的時候會多出了 $\frac{\epsilon}{h}$ 這麼多.

泰勒展開式

再深一點來說, 我們可以利用泰勒展開式知道 $$ \frac{f(a+h)-f(a)}{h} = f'(a) + \frac{h}{2}f''(\xi), \quad a\leq \xi \leq a+h. $$ 這個式子告訴我們, 用這方式算微分誤差應該是 $\frac{h}{2}f''(\xi) = O(h)$, 誤差會隨著 $h$ 減少而線性變小.

數學上我們有以上這個等式, 而數值計算上則是有以下這個等式 $$ \frac{f(a+h)-f(a) + \epsilon}{h} = f'(a) + \frac{h}{2}f''(\xi) + \frac{\epsilon}{h}, \quad a\leq \xi \leq a+h. $$ 也就是說, 真正的誤差公式應該是 $$ \frac{h}{2}f''(\xi) + \frac{\epsilon}{h}, $$ 當 $h$ 非常小的時候 $\frac{\epsilon}{h}$ 這項就會變很大.

比如說, 依我們之前所算的 $\epsilon\approx 10^{-16}$, 那當 $h=10^{-8}$ 時, 算出來的數字會多了大約 $\frac{10^{-16}}{10^{-8}} = 10^{-8}$.

而當 $h=10^{-14}$ 時, 算出來的數字會多了大約 $\frac{10^{-16}}{10^{-14}} = 10^{-2}$. 跟我們之前所發現的完全吻合!! 而這也就是為什麼當 $h$ 很靠近零的時候誤差會上升的原因.

Optimal $h$

那給定一個微分公式, 要怎麼知道 $h$ 能小到什麼程度呢? 一個簡單的感覺是這樣的, 由於誤差的第一項 $\frac{h}{2}f''(\xi)$ 會隨著 $h$ 變小而變小, 第二項 $\frac{\epsilon}{h}$ 則會變大, 因此整體最小值約會發生在兩項交叉時, 也就是當 $$ h \sim \frac{\epsilon}{h}, $$ (由於我們不知道 $f''(\xi)/2$ 是多少, 簡單起見設成 1). 上式稍微計算一下發現誤差最小值約發生在 $h=10^{-8}$, 誤差最小值則約為 $10^{-8}$.

Summary:

最後總結一下:

• 我們可以用數值計算來估計一個函數在某點的微分值 $$ f'(a) \approx \frac{f(a+h)-f(a)}{h}. $$ 這樣的做法稱為有限差分法 (finite difference method).

• 不過計算時 $h$ 值不能無限取小, 需考慮到捨入誤差的影響.