Matlab 教材：for-if-break 技術配合數值積分

採用另一個變數來紀錄「前一次迭代結果」的技巧非常基本且重要， 現在我們用常見的「梯形積分法」再示範一遍。

假設 f(x) 是一個在 [a,b] 閉區間內的連續函數， 則 f(x) 在 [a,b] 閉區間內的積分 (quadrature)， 可以用解析方法求得，也可以用數值方法估計。 所謂解析方法，就是我們在微積分課程中學習的： 先找到 f(x) 的反導函數 F(x)，然後計算 F(b) - F(a)。 但是許多 f(x) 是找不到 (基本函數型式的) F(x)， 因此解析方法就使不上力。 這時候只能用數值方法來估計積分值，這些方法通稱為 數值積分法 (numerical quadrature rules)。

現在我們要介紹一種常見的數值積分法：梯形積分法 (Trapezoidal quadrature rule)。 請參照下圖：先把 [a,b] 區間等分成 n 段 (n 是一個正整數)， 每一小段的兩端點形成一組節點 (nodal points)， 將這些節點編號、命名為 x0、x1、x2、...、 xn，很明顯地 x0 = a 而 xn = b。 第一段是 [x0, x1]、 第二段是 [x1, x2]、 第三段是 [x2, x3]、...、 第 n 段是 [xn-1, xn]。 每一段的寬度都一樣，就是 (b-a)/n，記做 $\Delta x$。

左圖 右圖

先假設 f(x) 是正值函數，如上圖，則
\begin{displaymath}\int_a^b f(x) dx\end{displaymath}

就「大約」是每一小段上形成的梯形面積和。這個概念，就是梯形積分法。

參照上面右側的圖，在第一段上有一個「高」為 $\Delta x$ 而兩「底」之寬度分別為 f(x0) 和 f(x1) 的梯形， 它的面積是 $\frac{\Delta x}2 (f(x_0)+f(x_1))$ 同理，在第二段上的梯形面積就是 $\frac{\Delta x}2 (f(x_1)+f(x_2))$ 在第三段上的梯形面積就是 $\frac{\Delta x}2 (f(x_2)+f(x_3))$ ... 在第 n 段上的梯形面積就是 $\frac{\Delta x}2 (f(x_{n-1})+f(x_n))$ 把這些小段上的面積加在一起，就是

\begin{displaymath}\int_a^b f(x)\,dx \approx \Delta x \left(\frac{f(x_0)}2+f(x_1)+f(x_2)+\cdots+f(x_{n-1})+\frac{f(x_n)}2\end{displaymath}

上式之右側，就是 n 段梯形積分估計值。 當 n 越來越大，也就是當 $\Delta x$ 越來越小，n 段梯形積分估計值就越來越接近解析的積分值。

即使 f(x) 不是正值的函數，我們還是用一樣的公式來計算 n 段梯形積分估計值。

舉例來說， $\sqrt{1+x^3}$ 就找不到基本函數型式的反導函數， 所以必須用數值積分法來估計它在 [0,1] 閉區間內的積分值：

\begin{displaymath}\int_0^1 \sqrt{1+x^3} dx\end{displaymath}

以下是用 Matlab 計算 100 段梯形積分估計值的指令：
x = linspace(0,1,101);
y = sqrt(1+x.^3);
h = 1/100;
T = h * (y(1)/2 + sum(y(2:100)) + y(101)/2)
以下解釋上述四句指令的用處：
  1. linspace(0,1,101) 將 [0,1] 區間等分成 100 段，取得 101 個節點，依序是
  2. y 就是 f(x) 在節點的函數值，依序是
  3. h 就是 $\Delta x$
  4. T 就是梯形積分估計值，其實，稍微更方便一點的寫法是
    T = h * (sum(y) - y(1)/2 - y(101)/2)

上面示範了如何計算特定段數的梯形積分估計值。 但是特定的段數並無意義，因為我們不知道那個估計值究竟有多準確？ 為了敘述的方便，我們把 n 段梯形積分估計值寫成以下記號：

\begin{displaymath}\int_a^b f(x) dx \approx T_n(f,a,b)\end{displaymath}

在計算的技巧上，我們通常要預設一個「容忍誤差」err， 並計算 T1(f,a,b), T2(f,a,b), T4(f,a,b), T8(f,a,b), T16(f,a,b), ..., 一直到
\begin{displaymath}\vert T_{2^{k-1}}(f,a,b) - T_{2^k}(f,a,b)\vert \leq \hbox{err}\end{displaymath}

為止。到那時，我們就把 $T_{2^k}(f,a,b)$ 當作是 f(x) 在 [a,b] 區間內的積分估計值。

以下就是一批實現上述想法的 Matlab 指令，我們假設最多分成 220 段。 

T = (1+sqrt(2))/2;
for k = 1:20
    Told = T;
    n = 2^k;
    x = linspace(0,1,n+1);
    y = sqrt(1+x.^3);
    h = 1/n;
    T = h * (sum(y) - y(1)/2 - y(n+1)/2);
    if ( abs(Told - T) <= err )
        break;
    end
end
以上程式可以直接寫在 Matlab 介面內，但是最好寫成一個腳本程式。 其中第一條指令就是計算 T1(f,0,1)，請讀者自己驗證。

可見得，當上述迭代結束的時候，要不是 k 達到了 20， 就是相鄰兩次的梯形積分估計值的差異，已經落在容忍誤差之內。 假設 err = 1e-8 也就是 10-8，則上述程式結束時， k 的值是 13 而 T 是 1.11144797184967， 我們應該可以相信這個 T 和積分的解析值相差在 10-7 以下， 所以 1.1114480 應該是一個正確的估計值。

習題

  1. 用梯形積分法求以下積分之估計值，確保小數點下前四位數是正確的估計值：
    \begin{displaymath}\int_0^2 \sqrt{1+x^3}\,dx\end{displaymath}

  2. 用梯形積分法求以下積分之估計值，確保小數點下前四位數是正確的估計值：
    \begin{displaymath}\int_{-1}^1 \sqrt{1+x^3}\,dx\end{displaymath}

  3. 用梯形積分法求以下積分之估計值，確保小數點下前四位數是正確的估計值：
    \begin{displaymath}\int_0^1 e^{-x^2}\,dx\end{displaymath}

  4. 用梯形積分法求以下積分之估計值，確保小數點下前四位數是正確的估計值：
    \begin{displaymath}\int_{-\frac1{\sqrt2}}^{\frac1{\sqrt2}} e^{-x^2}\,dx \end{displaymath}

單維彰 (2003/06/13) --- [Prev] [Next] [Up]