計算數學入門系列 - 數值表示方式（三）

 

另外一個比較特別的，是可以看到在電腦裏面可以表示的數字並不是平均分布。上面我們可以看到在0.0 旁邊的數字，要比在1.0旁邊的密 （距離1.0 最近的數字所差距只有大約\(10^{-7}\)）。 而這個在1.0附近的數字密度就大大影響着所有在電腦裏面計算的準確度。 想像一下，如果我們要計算x+y = x (1+y/x)，如果y/x 是少於 \(\epsilon\)，電腦將無法把1+y/x 準確計算出來。所以x+y就等於x。

由於上面這兩個關於數值分布的特性，在我們編寫程式時有些時候也要特別留意。其中一個例子是同學在編寫偏微分方程的數值解的時候，有時會出現下面的情況。 （下面用到的，是一些MATLAB 的程式，在其他編程語言裏面問題也差不多。）

finalT=…;

N=…; % number of time iterations

dt=finalT/N;

time=0;

while (time~=finalT)

      time=time+dt;

      …

end

上面這個結構，很多時候我們在解決跟時間有關的偏微分方程時都會用到。我們希望找到在finalT 時候的數值解。由於一些穩定性的要求，我們必須把從開始的時間到最後時間分隔成N分，然後也要把中間的所有數值解也找出來，而時間與時間之間的相隔用dt 來表示。 這個程式的問題是在於time 這個variable 會不斷改變，而這個不斷改變的time，我們會期待在N次運算以後剛好等於finalT，然後中間的while-loop 就會完成。這個疊代 （Iteration）過程完結以後，我們得到的數值解就是我們希望找到的答案。可是由於數字很可能無法被正確表示，有可能time 和finalT 是沒有辦法相等的。然後，這個程式就loop 死了。 那有什麼辦法解決呢？方法有很多，既然已經知道需要做N次疊代，為什麼不用它作為疊代 的index 呢？

finalT=…;

N=…; % number of time iterations

dt=finalT/N;

time=0;

for i=1:N

      time=time+dt;

      …

end

 

另外還有些情況，穩定性的要求有機會需要隨着時間而有所改變，這樣子我們就沒有辦法在程式開始時知道必須要做多少次的疊代，上面這個問題就需要有其他的解決方法。

finalT=…;

time=0;

flag=0;

while (flag==1)

      dt=…;

      if (time+dt>=finalT)

            dt=finalT-time;

            flag=1;

      end

      time=time+dt;

      …

end

需要解決的問題有兩個，第一我們如何知道加進來的時間dt 會否令到time 大於我們所需要的finalT，第二如何令到程式停止在需要的時間。要解決這個問題我們必須要有方法調整dt，令到他不會使time 太大。所以我們就可以在程式裏面多加一個if 的condition，使到程式做完最後這一次的疊代，就會停在我們所需要的時間。

另外一個關於machine epsilon 的影響，是在於我們用一些疊代方法去找一些數學問題的近似值。 比如說上面提到要找f(x)=0 的答案，我們會設計一些方法去產生一個序列，先從我們給予的一個猜想x(0)，找出x(1)，然後希望這個序列會越來越接近f(x) 的 根。當然我們不可能做無限次疊代，所以這個過程我們在程式裏面需要有一個開關，當我們猜想找出來的答案已經很接近實質數值時，程式就需要停止。這個開關有很多 不同的方法，同學可以隨便選用。其中一個是說接待過程對於找出來的答案已經不會有太大改變，數學上來說，我們覺得|x(n+1)-x(n)| 少於某個tol 的時候，我們就可以把疊代過程停止。那問題是，如何挑選這個tol 呢？很多同學都希望答案越準確越好， 可能把這個tol 設定成\(10^{-20}\)。可是這根本不需要，因為我們根本沒有辦法保證中間運算過程裏面的加減乘除誤差不會少於這個數值，所以只要答案跟machine epsilon 差不多其實我們已經可以把程式停止。