Matlab 教材：for-if-break 技術配合牛頓法

我們再回顧所謂的「牛頓法」。這是以數值計算估計 f(x) = 0 之一個根的方法。 假設 f(x) 可微，且在某個數 x0 的「附近」有一個根， 則根據以下迭代公式

\begin{displaymath} x_n = x_{n-1} - \frac{f(x_{n-1})}{f'(x_{n-1})},\quad n=1,2,3,\ldots \end{displaymath}

可以計算出來一個數列
x1, x1, x2, ...

理論上這是一個無窮數列。如果 x0 挑得「夠好」 (詳細的數學要在其他領域中學習)，那麼此數列將會收斂：
\begin{displaymath} \lim_{n\to\infty} x_n = p \end{displaymath}

收斂的值 p 將會是一個根，亦即 f(p)=0。

問題是，用電腦做計算總是無法真的計算無窮多項，因此不能真的算出極限。 我們再度以 f(x) = x2 - 2 為例，用牛頓法計算正根 $\sqrt2$。 則迭代公式為

\begin{displaymath} x_n = \frac{x_{n-1}}2 + \frac1{x_{n-1}},\quad n=1,2,3,\ldots \end{displaymath}

(這個公式和 [用 Matlab 熟悉變數置換] 裡面寫的其實一樣， 只是足標的符號不同而已。) 如果我們要計算 x3，可以按向上箭頭來重複執行， 也可以用以下 for 迴圈辦到： 
1 x = 1;
2 for n = 1:3,
3     x = x/2 + 1/x;
4 end
以上程式可以直接寫在 Matlab 介面內， 但是最好寫成一個腳本程式。 在以上程式內，紅色的數目代表列號，稍後再說，它們並不是程式的一部分。 而我們假設 x0 = 1，並且採用「變數置換」的技巧， 不必使用三個變數 x1, x2, x3， 而是重複使用一個變數 x。這應該是我們起碼已經熟悉的技巧。

當 Matlab 執行上述四列指令時，變數 n 和 x 的變化如下：
列號n  x  
1 未使用x0
2 1x0
3 1x1
4 迴圈結束，返回 2
2 2x1
3 2x2
4 迴圈結束，返回 2
2 3x2
3 3x3
4 迴圈結束，返回 2，
但是發現沒有下一個數要
代入 n，所以完成迴圈
5 輸入下一步指令

從上述表格中觀察，得知每當迴圈結束，x 的值就是 xn， 而 n 就是當時代入的值。所以 

x = 1;
for n = 1:10,
    x = x/2 + 1/x;
end
完成之後，就計算出來 x10。讀者注意不要忘記， 每次重新開始迭代之前，要記得把 x 指派成初始值。

以上的說明，告訴我們如何計算 xn，但是沒有說明如何讓 Matlab 自動檢查 xn 是否已經足夠靠近 p。 一般而言，我們並不知道 p，所以不可能計算 |xn - p|。 因此在計算的技巧上，我們通常要預設一個「容忍誤差」err，如果

\begin{displaymath} \vert x_n - x_{n-1}\vert \leq \hbox{err} \end{displaymath}

我們就認為 xn 在可計算的範圍內已經「達到」極限， 因此可以停止迭代，而取 xn 作為 p。 為了要計算 |xn - xn-1|， 我們必須用另一個變數來儲存「前一次迭代」的結果。以下程式是一個範例； 同樣地，紅色數目字代表指令的編號，不是指令的一部分。 
1 x = 1;
2 for n = 1:10
3     xold = x;
4     x = x/2 + 1/x;
5     if ( abs(xold - x) <= err )
6         break;
7     end
8 end

我們再一步步地審視 Matlab 的執行狀況：
列號n  x  xold
1 未使用x0未使用
2 1x0未使用
3 1x0x0
4 1x1x0
5 變數都沒被更動，只是檢查條件。假如 False，則 6 沒執行
7 if 結束，變數都沒被更動
8 迴圈結束，返回 2
2 2x1x0
3 2x1x1
4 2x2x1
5 變數都沒被更動，只是檢查條件。假如 False，則 6 沒執行
7 if 結束，變數都沒被更動
8 迴圈結束，返回 2
...
2 6x5x4
3 6x5x5
4 6x6x5
5 變數都沒被更動，只是檢查條件。假如 True
6 提早完成迴圈，跳到 8 的下一個指令
9 輸入下一步指令

可見得，當上述迭代結束的時候，要不是 n 達到了 10， 就是 |xn - xn-1| <= err。 我們可以檢查 n 的值。如果它是 10，則有可能 |xn - xn-1| 沒有達到容忍誤差，否則就代表 |xn - xn-1| 達到了容忍誤差， 而我們可以將 xn 視為極限 p。

假設 err = 5e-15 也就是 5 * 10-15，則上述程式結束時， n 的值是 6 而 abs(x - sqrt(2)) 是 2.2204e-16。

採用另一個變數來紀錄「前一次迭代結果」的技巧非常基本且重要，要能夠流利使用。

習題

  1. 令 x0 = -1，以牛頓法求
    \begin{displaymath} f(x) = e^x - \cos x \end{displaymath}

    的根，預設容忍誤差是 5e-15。輸出第一個達到此要求的 n 和 xn 以及 f(xn)。
  2. 令 x0 = 2，以牛頓法求
    \begin{displaymath} f(x) = x^3 - 3x -1 \end{displaymath}

    的根，預設容忍誤差是 5e-15。輸出第一個達到此要求的 n 和 xn 以及 f(xn)。
  3. 令 x0 = 0，以牛頓法求
    \begin{displaymath} f(x) = \tan x - \cos x \end{displaymath}

    的根，預設容忍誤差是 5e-15。輸出第一個達到此要求的 n 和 xn 以及 f(xn)。
  4. 令 x0 = i (單位虛數)，以牛頓法求
    \begin{displaymath} f(x) = x^2 + 2 \end{displaymath}

    的根，預設容忍誤差是 5e-15。輸出第一個達到此要求的 n 和 xn 以及 f(xn)。
單維彰 (2003/04/23) --- [Prev] [Next] [Up]