用電腦算極限

這裡我們要介紹如何用電腦算極限, 以及我們來看一下當我們真的這樣做的時候有可能會發生什麼問題. 我們以 $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$ 這個點, 所以沒有問題. 相反的, 如果取奇數個點就一定會取到 $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 *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 \ge 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}$ 這個等級.

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$ 函數的性質可以參考 wiki 上的介紹.

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

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