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

牛頓求根迭代法

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 = axn = 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),輸出每一對的 pt(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;
}

以上程式仍然是在技巧上還有改進之處。 此處的主要目的,是示範如何利用計算機從事數學實驗。 從這一千對的 pt(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 單維彰