Matlab 教材：用 Matlab 熟悉變數置換

初學電腦程式語言的人，通常遇到的第一個概念難關，就是變數的數值置換。 例如我們先說

x = 1;
然後說
x = 2*x
數學敏感的讀者看到上面那條「等式」，立刻會想：x 勢必等於 0。 其實不是。

因為 = 是當「指派」指令，所以等號的右邊要先做計算， 計算之後的結果是 2，然後再把 2 指派給 x，所以執行之後， x 的值就從 1 變成了 2，這也就是變數的數值置換。 再來，如果現在 x 的數值是 2，執行了

x = x + 2
之後 (數學敏感的讀者，直接反應或許是「無解」或「矛盾」)，x 的值就換成了 4。

等號的左右兩邊不但可以出現同一個變數，在右邊還可以出現一次以上， 例如

x = 4;
x = x * sqrt(x)
所有出現在等號右邊的 x 都取出舊的值，也就是 4。 計算得到結果之後，才指派給 x，也就是 8。

變數的置換技巧，可以拿來做許多有趣的計算。例如想要計算

\begin{displaymath} 1+\frac1{2} +\frac1{2^2} +\frac1{2^3} +\frac1{2^4} +\frac1{2^5} +\frac1{2^6} \end{displaymath}

可以轉換成一個數學迭代的程序：
\begin{displaymath} \left\{\begin{array}{l} x_0 = 1 \\ x_n = x_{n-1} + \frac1{2^n} \end{array}\right. \quad n=1,2,\ldots,6 \end{displaymath}

按照數學寫法，我們需要七個變數：x0, x1, x2, x3, x4, x5, x6。 如果運用電腦程式的變數置換技巧，只需要一個變數 x 就夠了，如下：
x = 1;
x = x + 1/2^1
x = x + 1/2^2
x = x + 1/2^3
x = x + 1/2^4
x = x + 1/2^5
x = x + 1/2^6
我們可以利用向上箭頭，重複呼叫一樣形式的指令， 每次只要修改一個數字即可重複執行上述計算，並且觀察計算的結果。 我得到的結果是 1.9844。

事實上，我們都知道以下無窮等比級數

\begin{displaymath} 1+\frac1{2} +\frac1{2^2} +\frac1{2^3} +\cdots \end{displaymath}

的極限是 2，若持續以上的迭代 (每次升高一個次方)，一直做下去， 就會一步一步看到收斂的過程。做到 15 次方的時候，看到的數值是 2.0000， 然後即使再繼續做，這數值也不會改變。

再舉一個例子，若 f(x) 是一個可微函數，我們知道牛頓求根算法說的是


其中 x0 是一個初始猜想，它必須是一個已經「頗靠近」真解的某個數。 而迭代的過程將會產生新變數 x1, x2, x3, ...； 如果運氣不太壞的話， xn 就會越來越接近 f(x)=0 的一個根， 換句話說，f( xn ) 越來越接近零。

用個實際點的例子，譬如 f(x) = x2 - 2， 那麼我們可以利用牛頓法求他的一個正根，也就是 sqrt(2)。 將 f(x) 和 f'(x) 代入牛頓迭代過程，整理得

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

我們令 x0 是 1，則迭代
x = 1;
x = x/2 + 1/x
x = x/2 + 1/x
x = x/2 + 1/x
x = x/2 + 1/x
就看到答案已經固定在 1.4142 不再改變。

習題

  1. 計算以下算式 (答案寫出五位有效數字)
    \begin{displaymath} 1+\frac1{2^2} +\frac1{3^2} +\frac1{4^2} +\frac1{5^2} +\frac1{6^2} +\frac1{7^2} \end{displaymath}

  2. 以下極限存在，請寫出 x 的前五位有效數字。
    \begin{displaymath} x = 1 + \frac1{2^2} + \frac1{3^2} + \frac1{4^2} + \cdots\end{displaymath}

  3. 令 f(x) = x2 - 2，舉出一個 x0 使得牛頓迭代法收斂到 -sqrt(2) 這個根。
  4. 令 f(x) = x2 + 2，舉出一個 x0 使得牛頓迭代法收斂到 sqrt(2)*i 這個根。
  5. 令 x 是使得
    exp(x) = cos(x)

    的最大負數，寫出 x 的前五位有效數字。
單維彰 (03/03/08) --- [Prev] [Next] [Up]