藍影

在了解 如何判定解區間 後
接下來要找的是
在這個解區間裡要如何找出正確解
可以用最笨的暴力法

do{
   x := x + step;   
}while(f(x) > error_rate);

不過暴力法花的時間太長了,
如果 error_rate 設小一點, 真的會跳很久
在此先介紹二分法 bi-sector method.

二分法原理是一直取 low 和 up 之中間值 m,
m = (low + up) / 2
然後不停的帶入, 分成以下情況討論

(1) f(m)*f(low) <0:
    有一解於 m 與 low 之間, 故將 up 設為 m
    再將 m 設為 (up+low)/2

(2) f(m)*f(up) <0:
    有一解於 m 與 up 之間, 故將 low 設為 m
    再將 m 設為 (up+low)/2

(3) 以上二種情形都不符合的時候, 代表
    其實這個區間裡面並沒有解,
    在第一個迭代裡面就可以發現出來

基於以上說明, 通常我們都會再設一個 "精準度"
如果根是無理數, 那麼我們只可以要求去 "逼近"
沒辦法完整的用小數去表示,
於是我們設一個最小誤差的變數為10的-4次方:

#define SMALL (double)(10E-4)

當f(m) 的值小於 SMALL 時,
便代表其誤差是我們可以接受的,
這時便可跳出回圈..

先以 f(x) = x^3 + 2x^2 -5x +1 ,
在 x = (1,2) 中找根為例,
在此只說明幾項, 其餘可自行推演

==========================
0 times:
L   :  1.0000,H   :  2.0000,X    :  1.5000
f(L): -1.0000,f(H):  7.0000,f(X) :  1.3750

-> f(X)*f(L) < 0,
H = X = 1.5
X = (1.5+1)/2 = 1.25
==========================
1 times:
L   :  1.0000,H   :  1.5000,X    :  1.2500
f(L): -1.0000,f(H):  1.3750,f(X) : -0.1719

-> f(X)*f(H) < 0,
L = X = 1.25
X = (1.25+5)/2 = 1.375
==========================
2 times:
L   :  1.2500,H   :  1.5000,X    :  1.3750
f(L): -0.1719,f(H):  1.3750,f(X) :  0.5059
....
....
==========================
7 times:
L   :  1.2813,H   :  1.2891,X    :  1.2852
f(L): -0.0197,f(H):  0.0201,f(X) :  0.0001
==========================

最後算出來, x = 1.2852

以下程式碼請參閱, 若有問題歡迎回覆討論

// ==================================
// filename: BiSector.cpp
// use bisector method to find solution
// author  : Edison.Shih.
// Date    : 2010.3.6
// ** all rights resever **
// ==================================

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

#define SMALL (double)(10E-4)
#define NO_FIND_NUM      (double)(-9999.0)

// =======================================
// f(x) = x^3 + 2x^2 -5x +1,
// have sol at (-4, -3), (0,1), (1, 2)
double func(double x){
      return (x*x*x + 2*x*x - 5*x + 1);
}

// =======================================
double bisector(double (*fptr)(double),
                        double low,
                        double high,
                        double small)
{// small: error rate
      double x = (low+high) / 2.0;
      double xvalue = 0.0; // f(x) = f( (h+low)/2 )
      double hvalue = 0.0; // f(h)
      double lvalue = 0.0; // f(low)
     
      int times = 0;
      //printf("==========================\n");
      do{
            xvalue = (*fptr)(x);
            hvalue = (*fptr)(high);
            lvalue = (*fptr)(low);
            /*
            printf("==========================\n");
            printf("%d times:\n", times++);
     
            printf("L   :%8.4lf,H   :%8.4lf,X    :%8.4lf\n",
                  low, high, x);
            printf("f(L):%8.4lf,f(H):%8.4lf,f(X) :%8.4lf\n",
                  lvalue, hvalue, xvalue);           
            */

            if(xvalue>0 && lvalue>0 && hvalue<0) low = x;
            else if(xvalue>0 && lvalue<0 && hvalue>0) high = x;
            else if(xvalue<0 && lvalue>0 && hvalue<0) high = x;
            else if(xvalue<0 && lvalue<0 && hvalue>0) low = x;
            else return (NO_FIND_NUM); // no find ans
            x = (high + low) / 2;           
      }while(fabs(xvalue) > small);
      // printf("==========================\n");

      return x;
}

// =======================================
int main(int argc, char **argv)
{
      double sol = 0.0;
      for(double x=-5.0; x<=5.0; x=x+1.0)
      {
            sol = bisector(func, x, x+1.0, SMALL);
            printf("sol at (%.0lf,%.0lf):",x,x+1);
            if(sol==NO_FIND_NUM) {
                  printf("no find sol\n");
            }
            else {
                  printf("%.4lf\n",sol);
            }
      }
      return 0;
}

// =======================================
執行結果

sol at (-5,-4):no find sol
sol at (-4,-3):-3.5070
sol at (-3,-2):no find sol
sol at (-2,-1):no find sol
sol at (-1,0):no find sol
sol at (0,1):0.2222
sol at (1,2):1.2832
sol at (2,3):no find sol
sol at (3,4):no find sol
sol at (4,5):no find sol
sol at (5,6):no find sol