[回目錄]


圖形與說明


 

在給定函式 f(x) 之情況下,給定了下界 a 與上界 b,同時確定了 [a, b] 含有一解

(可使用勘根定理判別),以二分法方式求出 f(x) 於 [a,b] 區間之解。

 

一開始,以 a, b 值代入 f(x),得到 f(a), f(b) ,對應之圖形如下

BiSect01.png

 

接下來,選擇 a, b 之中點 c,即 c = (a+b) / 2,再求 f(c) 之值。

 

BiSect02.png

 

注意到上圖,這裡又是一次勘根定理的應用,算出 f(c) 後,

若 f(a) * f(c) < 0 ,則代表根位在 [a,c] 之間,所以把下界 b 值調整為 c,

若 f(b) * f(c) > 0,則代表根位在 [ c, b] 之間,此時該把上界 a 值調成 c。

調完之後把原本的 b 擦掉,如下圖

 

BiSect03.png

 

接下來繼續找 c = (a+b) / 2,並算出 f(c) 值


  BiSect04.png  

 

和上面一樣,又是一次勘根定理的應用,算出 f(c) 後,

若 f(a) * f(c) < 0 ,則代表根位在 [a,c] 之間,所以把下界 b 值調整為 c,

若 f(b) * f(c) > 0,則代表根位在 [ c, b] 之間,此時該把上界 a 值調成 c。

這次若在 [c, b] 間,於是將 a 換成 c。

 

BiSect05.png

 

重覆以上動作,一直到 f(c) 小於某個可接受誤差 eps,

就將 c 點視為此方程式之一組解,此算法即可停止。

 

要提醒的是,二分法在執行前,一定要選確定在 [a, b] 區間有根,

做下去才有意義。

 

終止條件

 

上述是將 f(c) 小於某個可接受誤差 eps 視為終止條件,

這只是終止條件之其中一種方式。而上述的小於某個誤差值,

這種終止方式叫「收斂」,因其解答已在可接收的誤差範圍內。

筆者見過的終止條件有下述幾種

 

(1) f(x) <= eps 視為收斂

如上所述,這種方式是較為在意求到之解代入函數後,誤差必須極小,但不一定可以收斂到 eps 這麼小的數。

等號非常重要,在大多情況下都會加,因大多情況下會將 eps 設非常小,可能為 1E-6 ~ 1E-12,甚至

可能有心人士誤設為 0 ,為避免這種情況,筆者建議要加上等號。

 

(2) 絕對值 (x2-x1) <= eps 視為收斂

將上一次求得之 x1 與這次求得之 x2 相減,若小於等於 eps 視為收斂,

當然也有可能會一直都大於 eps ,但結果已算可接受程式。

 

(3) 絕對值 (x2-x1)/x2 <= eps 視為收斂

若解之變動率小於一數值,則視為收斂。

 

(4) 絕對值 f(x2)-f(x1) <= eps 視為收斂

將上一次求得之 f(x1) 與 這次求得之 f(x2) 相減,若小於等於 eps 視為收斂,

若該函式在求解區間之斜率非常大的話,幾乎不可能收斂。

 

(5) 強制終止方式

上述四種方式,有時會混用,如挑 1, 2 混用 或 1, 4 混用等等,

只要有一條件符合,即視為收斂情形。但在某些函數情況下,

求解過程可能無法順利收斂,3 種情況可能都無法達成,

於是會再用一個強制終止之條件,最常見到的方式是,限制最大迭代次數。

假設最大迭代次數為 100,進行 100 次還無法收斂的話即強制停止求值動作。

 

虛擬碼


Algorithm BiSect

    E0 : 初始化最小誤差 EPS,初始化上界 a,下界 b。

    E1 :  c = (a+b)  / 2

    E2 :  if abs ( f(c)  ) <= EPS ,演算法結束,傳回 c 值。

    E3 :  if ( f(c) * f(a) <= 0)   b = c;
            else  a = c;

    E4 :  goto E1

End Algorithm

 

C語言初階程式碼

 

底下程式碼是依上述虛擬碼直譯而言。

 

double BiSector(double low,              /*   下界  */
                 double up,              /*   上界  */
                 double (*fx)(double),   /* 適應函式*/
                 double eps,             /* 容許誤差*/
                 int max_itera)          /* 最大迭代*/
{
     double mid;
     do{
           mid = (low+up)*0.5;
           if( fx(mid) * fx(up) <= 0.0) low = mid;
           else up = mid;
           --max_itera;
     }while(fabs( fx(mid)) > eps && max_itera!=0);
     return mid;
}


大致上該有的變數都有,結果也正常,但並不是很好的程式碼。

較大問題在於調用 fx 次數太多,在實際問題上,

有可能非線性方程式代入一次值 ( 上述 fx(mid)、fx(up) 都是在代值)

之時間可能不小,故在設計時,盡可能使用暫存變數,

以取代調用 fx 所花費之時間成本。完整程式碼如下

 

C 語言程式碼

 

/*******************************************************************/
/*                                                                 */
/*     filename : BiSector.c                                       */
/*     author   : edison.shih/edisonx                              */
/*     compiler : Visual C++ 2008                                  */
/*     date     : 2011.03.07                                       */
/*                                                                 */
/*         A.L.L.      R.I.G.H.T.S.     R.E.S.E.R.V.E.             */
/*                                                                 */
/*******************************************************************/

/*----------------------------------------------------------------*\
|
|  Assume there is a root at [low, up] in f(x)
|
|  (1) yup = f(up)
|
|  (2) mid <- (low + up) / 2 , y <- f(mid)
|
|  (3) if fabs(y) < eps ---> algorithm terminate, mid is answer
|
|  (4) if y1 * f(up) < 0  ( f(mid) * f(up) < 0 ) ---> low = mid
|      else ---> up = mid , yup = y
|
|  (5) goto step (2)

\*----------------------------------------------------------------*/

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

// [ -2.00 , -1.00 ] , [ 2.00 , 3.00 ] , [ +4.00 , +5.00 ]
double func(double x)
{
     double x2=x*x, x3=x2*x;
     return (x3 - 5.48*x2 -  1.4883*x + 20.394828);
}

// -------------------------------------------------------
double BiSector(double low,             /*  
下界  */
                double up,              /*   上界  */
                double (*fx)(double),   /* 適應函式*/
                double eps,             /* 容許誤差*/
                int max_itera)          /* 最大迭代*/
{
     double mid, y;
     double yup = fx(up);

     if( yup *fx(low) > 0.0) {
           printf("\n>   has no root at [%lf, %lf]", low, up);
           return low;
     }
     do{
           mid = (low+up)*0.5;
           y   = fx(mid);
           if(y * yup <= 0.0) low = mid;
           else up = mid, yup=y;
           --max_itera;
     }while(max_itera && fabs(y) > eps);
     return mid;
}

int main()
{
     double low, up, x;
     const int max_itera=100;
     const double eps = 1E-9;

     low = -2.0 , up = -1.0;
     x = BiSector(low, upfunc, eps, max_itera);
     printf("\n> func(%+e) = %+e", x, func(x));

     low = +2.0 , up = +3.0;
     x = BiSector(low, upfunc, eps, max_itera);
     printf("\n> func(%+e) = %+e", x, func(x));

     low = +4.0 , up = +5.0;
     x = BiSector(low, upfunc, eps, max_itera);
     printf("\n> func(%+e) = %+e", x, func(x));

     low = -10.0 , up = -20.0; // for test, can't find root.
     x = BiSector(low, upfunc, eps, max_itera);
     printf("\n> func(%+e) = %+e", x, func(x));
     return 0;
}

 

執行結果

 
> func(-1.781504e+000) = -5.953638e-011
> func(+2.313838e+000) = +2.521503e-010
> func(+4.947666e+000) = +1.224798e-010
>   has no root at [-10.000000, -20.000000]
> func(-1.000000e+001) = -1.512722e+003

 

其他補充說明

 

這份程式碼還不是最好的設計模式,原因如下述

(a) 呼叫副函式結束時,只傳回 x 值,沒記下 y 值,回到 main 裡後要看函數值還要再調用一次,浪費時間成本。

(b) 傳回值無法判別是收斂結束,還是迭代次數過長強制結束。

 

此系列文章只做初步示範,必會留下一些思考與改進空間,往後此系列文亦如此,

有興趣之讀者可實際再對於函式做適度之修改。

 

至於常見到的 

(max_itera && fabs(y) > eps);

是不是要寫成

(max_itera && (y>eps || y<-eps));

這部份太過細節鎖碎,且可爭議探討部份實在太多,

因非主題原故,此處不再深入研究探討。

 

[回目錄]

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()