迴圈應用: 二分勘根法

1. 靈活地運用函數來組織程式

2. 迴圈應用設計

二分勘根法

有的時候要你找一個方程式 f(x) = 0 的解時, 你沒有辦法寫出它的 closed-form 解, 例如 f(x) = x² - 2 - cos(2 x) = 0. 這個時候我們可以用數值方法很快地去找到接近的 x 使得 f(x) 非常接近 0, 底下描述一個簡單的二分逼近法:

1. 由一個可能的根的區間 [x₁, x₂] 開始, 也就是 x₁, x₂ 滿足 f(x₁) * f(x₂) < 0
2. 令 x_left = x₁ 且 x_right = x₂
3. 然後測試 x_left 和 x_right 的中間點 x₃ = (x_left + x_right) / 2, 如果 f(x_left) * f(x₃) < 0, 則令 x_right = x₃ (也就是搜尋的區間變成 [x₁, x₃]), 否則令 x_left = x₃ (也就是搜尋的區間變成 [x₃, x₂])
4. 接下去我們還是測試 x_left 和 x_right 的中間點 x₄ = (x_left + x_right) / 2, 如果 f(x_left) * f(x₄) < 0, 令 x_right = x₄, 否則令 x_left = x₄
5. 如此不斷地重複
   ...
6. x_n = (x_left + x_right) / 2, 如果 f(x_left) * f(x_n) < 0, x_right = x_n, 否則令 x_left = x_n ...
7. 直到 f(x_n) is 0 或是 | x_right - x_left | < epsilon, epsilon 是一個指定的誤差值

上面說明中 x_left 和 x_right 可以作為程式中的變數, x_i 純粹是說明用的

以下圖片應該可以給你更多的直覺:

函數 x² - 2 和 cos(2 x)

函數 f(x) = x² - 2 - cos(2 x)

上圖中顯然根的範圍落在區間 [0, 2].

如何設計一個程式 (同時設計類似上面的演算法) 來得到接近的數值解?

第一步: 請你手動一步一步去做, 步驟越詳細越好, 不要遺漏任何步驟, 例如:

若 f(x) = x² - 2 - cos(2 x)

1. 由 x_left = x₁ = 0 開始, f(x_left) = -3, x_right = x₂ = 2, f(x_right) = 2.6536, f(x_left) * f(x_right) < 0
   
2. x₃ = (x_left + x_right) / 2 = (x₁ + x₂) / 2 = 1, f(x₃) = -0.5839, f(x₃) * f(x_right) < 0, x_left = x₃
   
3. x₄ = (x_left + x_right) / 2 = (x₃ + x₂) / 2 = 1.5 , f(x₄) = 1.24, f(x_left) * f(x₄) < 0, x_right = x₄
   
4. x₅ = (x_left + x_right) / 2 = (x₃ + x₄) / 2 = 1.75 , f(x₅) =1.999, f(x_left) * f(x₅) < 0, x_right = x₅
   
5. x₆ = (x_left + x_right) / 2 = (x₃ + x₅) / 2 = 1.375 , f(x₆) =0.8149, f(x_left) * f(x₆) < 0, x_right = x₆
   
6. x₇ = (x_left + x_right) / 2 = (x₃ + x₆) / 2 = 1.1875, f(x₇) =0.1304, f(x_left) * f(x₇) < 0, x_right = x₇
   
7. x₈ = (x_left + x_right) / 2 = (x₃ + x₇) / 2 = 1.0938, f(x₈) =-0.2252, f(x₈) * f(x_right) < 0, x_left = x₈
   
8. x₉ = (x_left + x_right) / 2 = (x₈ + x₇) / 2 = 1.1406, f(x₉) = -0.0469, f(x₉) * f(x_right) < 0, x_left = x₉
   
9. x₁₀ = (x_left + x_right) / 2 = (x₉ + x₇) / 2 = ...
   

在上面的步驟裡, 你可以注意到兩個步驟間的變動其實很少, 常常這就是可以用 迴圈 來實作的地方, 在轉換成程式之前你需要把這些動作規律化, 你希望兩個步驟間重複的部份越多越好

第二步: 由上面手動計算過程中擷取 兩個步驟重複的共通部份 來得到 迴圈的主體

x_i = (x_left + x_right) / 2 = 1.75 , f(x_i) =1.999,
如果 f(x_left) * f(x_i) < 0 則 x_right= x_i否則如果 f(x_i) * f(x_right) < 0 則 x_left= x_i

在運作過程中這些 x_i 其實不需要分別是變數, 在計算 x_i+1 的時候, x_i已經記錄在 x_left 或是 x_right變數中了, 所以可以用單一的變數 x_mid來存放 x_i 的數值, 直接用 x_i+1 覆蓋 x_i 的數值就好了

第三步: 迴圈初始化

令 x_left = x₁= 0 及 x_right = x₂= 2

第四步: 迴圈控制變數更新

if f(x_left) * f(x_i) < 0 則 x_left 不變且 x_right = x_i

f(x_i) * f(x_right) < 0 則 x_left = x_i 且 x_right 不變

第五步: 迴圈結束條件

1. | f(x_i) | <= 10^-13 (在有限的有效位數下測試 f(x_i) == 0, 如果用 float 來設計的話就只能測試 | f(x_i) | <= 10^-7)

2. 或是 | x_left - x_right | < epsilon (請注意, 因為是二分逼近法, 所以迴圈執行 n 次以後, | x_left - x_right | 應該是 | x₁ - x₂ | * 2^-n, 是一個可以預測的數值, 所以也可以由 | x₁ - x₂ | 以及可以接受的誤差 epsilon 先算出迴圈需要執行的次數 n)

請注意上面這兩個測試對於一個平滑的曲線來說是重複的

範例執行程式 (尋找 f(x) = x² - 2 - cos(2 x) 在 [0,2] 區間中的根)

基於上面的分析, 撰寫下列 bisect 函數

double bisect(double x_left,  /* input  - endpoints of interval in */ 
               double x_right, /*          which to look for a root */ 
               double epsilon) /* input  - error tolerance          */

函數回傳 x 值使得 f(x) 的數值在 [-epsilon, epsilon] 區間中.

迴圈初始化:

x_left 是左邊的邊界, x_right 是右邊的邊界

x_mid = (x_left + x_right) / 2;

f_mid = f(x_mid);

f_left = f(x_left);

f_right = f(x_right);

迴圈結束條件設立:

while (fabs(f_mid) > 1e-13) 或是 (x_right - x_left > epsilon) 或是先計算出迴圈需要執行的次數 n, 用 for (i=0; i<n; i++)

迴圈主體:

x_mid = (x_left + x_right) / 2;

f_mid = f(x_mid);

迴圈控制條件修改:

if (f_mid * f_left < 0)
{
    x_right = x_mid;
    f_right = f_mid;
}
else if (f_mid * f_right < 0) 
{
    x_left = x_mid;
    f_left = f_mid;
}

請完成下列函式:

double bisect(double x_left,  /* input  - endpoints of interval in */ 
              double x_right, /*          which to look for a root */ 
              double epsilon) /* input  - error tolerance          */ 

{
   ...
}

進一步的程式修改

我們可以把每次迴圈執行過程的 x_left, x_right, x_mid, 和 f_mid 都列印在螢幕上, 如此可以快速地看到程式是否表現正常

我們可以把條件測試 (x_right - x_left < epsilon) 放進迴圈裡面, 並且設定一個表示狀態的變數 root_found

如果條件滿足, 我們列印出成功的訊息, 否則列印出下一個區間 [x_left, x_right]

我們可以移除不需要的條件測試 (f_mid * f_right < 0) 甚至移除不需要的變數 f_right

我們可以把想要求根的函數 f(x) 寫成一個 C 的函數:

double f(double x)
{
    return x * x - 2 - cos(2 * x);
} 

測試 bisect() 的 main() 函數如下:

#include <stdio.h>
#include <stdlib.h>


int main(void) 
{
    double epsilon,   /* error tolerance */ 
           root; 
    
    /*  Get error tolerance from user    */ 
    printf("\nEnter tolerance> "); 
    scanf("%lf", &epsilon); 
    
    root = bisect(0, 2, epsilon); 
    printf("\n   (%.7f)^2 + 2 = cos(2 * %.7f)\n", root, root); 

    system("pause");
    return(0); 
} 

線上繳交版本請參考下列說明:

1. 程式輸入四個倍精準浮點數 a, b, c, d 以指定多項式 f(x) = a x² + b x + c - cos(d x)
   
2. 程式再輸入一個倍精準浮點數 epsilon 以指定根 r 的函數值 f(r) 距離 0 的可以容忍的誤差範圍 [-epsilon, epsilon]
   
3. 程式再輸入兩個倍精準浮點數 x1 及 x2 代表根的範圍 (假設範圍內一定有一解)
   
4. 程式印出 r 的數值 (小數點後請印 5 位數)

範例執行程式