用電腦算極限

這裡我們要介紹如何用電腦算極限, 以及我們來看一下當我們真的這樣實作的時候有可能會發生什麼問題. 我們以 $sinc$ 函數為例來做說明.

sinc function

首先我們要介紹一個特別的函數, $sinc(x)$, 定義如下: $$ sinc(x) = \frac{\sin(x)}{x}, \quad x\ne 0. $$ 很明顯可以看出來當 $x=0$ 的時候分母會等於零, 在數學定義上是不確定的, 所以把 $x=0$ 這個點先拿掉.

為了觀察這個函數的行為我們可以把這個函數畫出來. 首先我們在 $[-20, 20]$ 這個區間取 $1000$ 個點, 然後帶入上面 $sinc$ 函數的定義, 再把所有點連起來看看.

Remark: 有件事需要先說明一下, 由於我們是在 $[-20, 20]$ 這個區間均勻的取偶數個點, 所以會有 $500$ 個正數以及 $500$ 個負數, 重點是保證不會取到 $x=0$ 這個點, 所以沒有問題, 避免了程式出現 DivisionByZero 的錯誤. 相反的, 如果取奇數個點就一定會取到 $x=0$, 那就會有函數無定義的問題了.

我們很輕易可以看出來, 連起來的線還蠻"光滑"的. 函數值在 $x=0$ 附近似乎不會趨近正無窮大或負無窮大, 也沒有跳躍的現象. 這暗示了雖然 $x=0$ 處於分母, 但該處可能存在一個有限的極限值.

接著我們試著在 $x=0$ 附近放大一點看看, 我們在 $[-0.1, 0.1]$ 這個區間取 $1000$ 個點, 然後帶入 $sinc$ 函數的定義再把它畫出來:

看起來真的很光滑!! 而且似乎當 $x$ 很靠近 $0$ 時, $sinc(x)$ 的值很靠近 $1$.

接著我們取一個會越來越靠近 $0$ 的數列, 然後看一下當把 $sinc$ 函數在這個數列的點上取值時, 其值會不會越來越靠近 $1$. 我們取以下數列:

 1.670170079024566e-5
 6.14421235332821e-6
 2.2603294069810542e-6
 8.315287191035679e-7
 3.059023205018258e-7
 1.1253517471925912e-7
 4.139937718785167e-8
 1.522997974471263e-8
 5.602796437537268e-9
 2.061153622438558e-9

算一下 $sinc$ 函數在這些點上面的值, 並觀察他的趨勢:

 0.9999999999535089
 0.9999999999937081
 0.9999999999991485
 0.9999999999998848
 0.9999999999999845
 0.9999999999999979
 0.9999999999999997
 1.0
 1.0
 1.0

赫然發現算到後來就等於 $1$ 了!! 所以我們發現

$$ \lim_{x\to 0} sinc(x) = \lim_{x\to 0} \frac{\sin(x)}{x} = 1. $$

不過有一點點詭異的是, 在剛剛的計算裡我們最多也只是取到離 $0$ 很近的點而已, 但是算出來的結果卻是 $1$. 難道不只是 $sinc(0)=1$, 我們也有 $sinc(2.061153622438558\times 10^{-9})=1$ 嗎?

事實上並不是這樣.

在數學上, $sinc(x)=1$ 只會在極限意義下成立(或是定義擴充後), 但在電腦的輸出中, 我們看到 $sinc(2.06\times 10^{-9})$ 竟然精確地等於 $1$. 這並不是因為數學公式變了, 而是因為電腦的精度限制. 電腦有所謂的捨入誤差(rounding error), 因為電腦需要用有限位數來表達無窮小數, 所以一定要捨棄後面的位數. 我們把算出來的數字減去 $1$ 看看:

 -4.649114426769074e-11
 -6.291855925155687e-12
 -8.515410598874951e-13
 -1.1524114995609125e-13
 -1.554312234475219e-14
 -2.1094237467877974e-15
 -3.3306690738754696e-16
  0.0
  0.0
  0.0

我們可以發現這個數字最小可以到大約 $10^{-16}$, 之後就變成 $0$ 了. 也就是說我們目前用個這個程式語言其捨入誤差大約就是 $10^{-16}$.

machine epsilon

一般我們會用 machine epsilon 這個數字來量化在電腦裡浮點運算的捨入誤差.

machine epsilon 的定義是, 考慮正數$\epsilon$, 使得 $ 1+ \epsilon > 1$ 中最小的那個稱之為 machine epsilon. 當然以數學來看這個 machine epsilon 必須等於零. 不過在電腦裡並不是這樣.

為了方便我們稍微改一下定義, 我們考慮 $\epsilon = 2^{-k}$ 這種形式, 然後 machine epsilon 一樣是使得 $ 1+ \epsilon \ge 1$ 中最小的那個. julia 程式如下:

s=1.0
for k=1:100
    s = s/2.0
    if (1.0+s)<=1.0
        s = 2.0*s
        println("k=", k-1, ",  eps=", s)
        return
    end
end

在我的電腦上我發現 $\epsilon = 2^{-52} = 2.220446049250313e-16$.

k=52,  eps=2.220446049250313e-16

所以的確捨入誤差大約是 $10^{-16}$ 這個等級.

Taylor’s expansion

為什麼當 $x\to 0$ 的時候 $sinc(x)\to 1$ 而不是別的數字? 我們可以用 Taylor’s expansion 來看.

當 $x$ 很小時, $$ \sin(x)= x - \frac{x^3}{6} + O(x^5). $$ 因此, 當 $x$ 很小時, $$ sinc(x)=\frac{\sin(x)}{x} \approx 1 - \frac{x^2}{6}. $$ 所以可以清楚地看到 $sinc(x)\to 1$ for $x\to 0$.

而從電腦的角度來說, 如果 $x\approx 4\times 10^{-8}$, 那

$$ \frac{x^2}{6}=\frac{(4\times 10^{-8})^2}{6} = \frac{16}{6}10^{-16}, $$ 比 machine epsilon 還要大一點點.

不過如果 $x\approx 2\times 10^{-8}$, 那

$$ \frac{x^2}{6}= \frac{(2\times 10^{-8})^2}{6} = \frac{4}{6}10^{-16}\approx 0.66\times 10^{-16}, $$

就比 machine epsilon 還要小. 在雙精度浮點數的運算中, 以 $1$ 減去這數字就沒有任何作用, 導致結果依然精確等於 $1$.

這完美解釋了為什麼在之前的數列計算中, 當 $x$ 足夠小時, 我們會赫然發現 $sinc(x)$ 的輸出變成了連續的三個 $1.0$.

Summary:

稍微總結一下目前我們看到的:

• 我們在微積分裡學過 $sinc(x\to 0)=1$, 所以我們可以定義 sinc 函數為 $$ sinc(x) = \begin{cases} \frac{\sin(x)}{x}, \quad x\ne 0, \\\
  1, \quad x=0. \end{cases} $$ 在這樣的定義之下 $sinc$ 函數是個連續函數.

  • 用程式跑數值發現類似的結果 $sinc(x\to 0)=1$, 不過這只能告訴我們觀察到的趨勢, 要有確定的答案還是需要更嚴謹的數學證明才能下結論.
  • 更多關於 $sinc$ 函數的性質可以參考 wiki 上的介紹.

• 在數值計算上有所謂的捨入誤差, 這是用有限位元來表達無限位數一定會有的差異.

• 可以用 machines epsilon 來量化捨入誤差.