使用函式做數學計算的範例

牛頓求根迭代法

若 f(x) 是一次可微且導函數連續的函數， 則牛頓求根法的迭代公式是

x_{n+1} = x_n - {f(x_n)\over f'(x_n)},\quad n=0,1,2,\ldots
其中 x0 是初始猜想值 (initial guess)。 牛頓求根法產生一個數列 x0, x1, x2, ... 這個數列不一定收斂，即使收斂也不一定收斂到 f(x) 的一個根。 但是通常如果 x0 很靠近所求的根，則 n 不必超過 10 就足夠使得 f(xn) 幾乎是 0。

以下這個程式，將 $f(x) = \sqrt x - \cos x$ 和它的導函數 $f'(x) = {1\over2\sqrt x} + \sin x$ 定義成兩個 C 語言函式 f() 和 df()， 然後用牛頓法求一個近似根。 程式中的所有浮點數都用雙精度，取 x0 = 0.5， 讓電腦迭代 10 次之後停下來，輸出每次的 xn。 這個程式需要數學函式。 

#include <stdio.h>
#include <math.h>
 
/* 示範牛頓求根迭代法  (newton-iteration.c) */
#define STEP 10
  
double f(double);
double df(double);
 
main() {
    int n;
    double x;
    x = 0.5;
    for (n=0; n < STEP; ++n) {
        x = x - f(x)/df(x);
        printf("x_%-2d = %17.14f\n", n+1, x);
    }
}
 
double f(double x) {
    return sqrt(x) - cos(x);
}
double df(double x) {
    return 1/(2*sqrt(x)) + sin(x);
}
我必須提醒讀者，以上程式在技巧上和數學問題的處理上，都還有改進之處。 例如我們沒有考慮發生零分母的可能。 它只是一個範例而已。 注意在 f() 和 df() 中，使用 return 的技巧。 執行這個程式，可以明顯地看到數字收斂的效果，其實在 x4 之後，迭代數值的前 14 位有效數字就已經不變了， 因此可以說前 14 位有效數字「已經收斂」了。

梯形數值積分法

若 f(x) 在 [a,b] 中連續，則可積。 梯形數值積分法 (trapezoidal rule) 是求定積分之近似值的一種方法， 其公式如下

\int_a^b f(x)\,dx \approx h \cdot \Bigl[\frac{f(x_0)}2 + \frac{f(x_n)}2 + \sum_{k=1}^{n-1} f(x_k)\Bigr]
其中 n 表示將 [a,b] 等分成 n 段，令
h = \frac{b-a}n
而且 x0 = a、 xn = b，一般而言 xk = a + k h。

以下這個程式，示範以梯形數值積分法求

\int_0^1 \sqrt{1+x^3}\,dx
的定積分近似值。 此程式將被積分函數定義成 C 函式 f()， 所有浮點數用雙精度計算， 並依序用 n=2, 4, 8, 16,..., 216 算梯形積分值， 輸出對應每個 n 的定積分近似值。 這個程式需要呼叫數學函式。 
#include <stdio.h>
#include <math.h>
  
/* 示範梯形數值積分法  (trapeziodal.c) */
#define STEP 16
 
double f(double);
 
main() {
    int n, k, i;
    double h, T, a, b, x;
    T = 0;
    a = 0;
    b = 1;
    for (n=1, i=0; i < STEP; ++i) {
        n = 2*n;
        h = (b-a)/n;
        x = a;
        T = f(x)/2;
        for (k=1; k<n ; ++k) {
           x = x + h;
           T = T + f(x);
        }
        T = T + f(b)/2;
        T = h*T;
        printf("T_%-5d = %17.14f.\n", n, T);
    }
}
 
double f(double x) {
    return sqrt(1+x*x*x);
}
我必須提醒讀者，以上程式在技巧上和數學問題的處理上，都還有改進之處。 此處只是示範而已。 0鶡瘜o個程式，可以清楚地看到梯形積分數值的收斂情況， 也就是說，隨著 n 的增加，Tn 就有越來越多的有效數字是不變的。而這些不再隨 n 而改變的有效數字， 就是「已經收斂」的數字，它們是可靠的答案。 例如，我們應該可以相信 T65536 已經準確到小數點下第 9 位。

Collatz 猜想

任取一個正整數 p，令 x0 = p 且定義 Collatz 數列如下：

x_{n+1} = \cases{{x_n\over 2}} & if $x_n$...
所謂的 Collatz 猜想，就是說數列 xn 必定會出現 1。根據 Collatz 數列的定義，一旦出現 1， 接下去就是 1,4,2,1,4,2,... 這樣的循環。 到寫教材的時候為止，Collatz 猜想還是個未解決的問題： 亦即，尚未證明它對，也沒找到它錯的反例。

以下這個程式，為 Collatz 猜想做一些實驗。 我們依序取 p = 1,2,3,...,1000，對每一個 p， 定義 t(p) 是使得 xt(p) = 1 *熙怳p整數。例如若 p=1 則因為 x0=1 所以 t(1) = 0。 若 p=2 則因為 x0=2 然後 x1=1 所以 t(2) = 1。 若 p=3 則因為

x0=3 x1=10 x2=5 x3=16 x4=8 x5=4 x6=2 x7=1
所以 t(3) = 7。依此類推。 以下這個程式，定義一個函式 collatz()，輸入正整數 p， 輸出整數 t(p)，輸出每一對的 p 和 t(p)，並輸出 t(p) 最大的情況。 
#include <stdio.h>
   
/* 示範 Collatz 猜想的實驗 (collatz.c) */
#define SIZE 1000

int collatz(int);

main () {
    int p, tp, maxp, maxtp;
    maxp = maxtp = 0;
    for (p=1; p<=SIZE; ++p) {
        tp = collatz(p);
        printf("%4d: %4d\n", p, tp);
        if (tp > maxtp) {
            maxp = p;
            maxtp = tp;
        }
    }
    printf("for p=0..%d, max t(%d) = %d\n", SIZE, maxp, maxtp);
}

int collatz(int p) {
    int x=p, n=0;
    while (x != 1) {
        if (x%2 == 1)
            x = 3*x + 1;
        else
            x = x/2;
        ++n;
    }
    return n;
}
以上程式仍然是在技巧上還有改進之處。 此處的主要目的，是示範如何利用計算機從事數學實驗。 從這一千對的 p 和 t(p) 看來，Collatz 猜想都是正確的， 事實上讀者可以做更多的測試，到目前為止沒有找到 Collatz 猜想的反例。

習題

  1. 用牛頓求根法迭代 10 次， 計算 e-x = sin x 最靠近 0 的近似根。 說明這個近似根有幾位準確的有效數字。
  2. 用 216 個等分點的梯形數值積分法，估計以下定積分
    \int_0^1 \sqrt{1+9x^4}\,dx
    說明這個近似值有幾位準確的有效數字。
  3. 令 p=5,6,7,...,15, 列出當 x0 = p 所產生的 Collatz 數列，直到第一次出現 1 為止。

[ 前一節 ]‧[ 後一節 ]‧[ 回目錄 ]


注意：此處所有文件均為原著，個別的版權宣告日後會一一公布， 整體版面設計亦尚未完成。但仍請勿抄襲文字與圖片，以免觸犯著作權法。

Created: Feb 24, 2001
Last Revised: Feb 25, 2001
© Copyright 2001 Wei-Chang Shann 單維彰