圖形與說明
在給定函式 f(x) 之情況下，給定了下界 a 與上界 b，同時確定了 [a, b] 含有一解
(可使用勘根定理判別)，以二分法方式求出 f(x) 於 [a,b] 區間之解。
一開始，以 a, b 值代入 f(x)，得到 f(a), f(b) ，對應之圖形如下
接下來，選擇 a, b 之中點 c，即 c = (a+b) / 2，再求 f(c) 之值。
注意到上圖，這裡又是一次勘根定理的應用，算出 f(c) 後，
若 f(a) * f(c) < 0 ，則代表根位在 [a,c] 之間，所以把下界 b 值調整為 c，
若 f(b) * f(c) > 0，則代表根位在 [ c, b] 之間，此時該把上界 a 值調成 c。
調完之後把原本的 b 擦掉，如下圖
接下來繼續找 c = (a+b) / 2，並算出 f(c) 值
和上面一樣，又是一次勘根定理的應用，算出 f(c) 後，
若 f(a) * f(c) < 0 ，則代表根位在 [a,c] 之間，所以把下界 b 值調整為 c，
若 f(b) * f(c) > 0，則代表根位在 [ c, b] 之間，此時該把上界 a 值調成 c。
這次若在 [c, b] 間，於是將 a 換成 c。
重覆以上動作，一直到 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,
up, func, eps, max_itera);
printf("\n> func(%+e) = %+e", x, func(x));
low
= +2.0 , up = +3.0;
x
= BiSector(low,
up, func, eps, max_itera);
printf("\n> func(%+e) = %+e", x, func(x));
low
= +4.0 , up = +5.0;
x
= BiSector(low,
up, func, 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,
up, func, 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));
這部份太過細節鎖碎，且可爭議探討部份實在太多，
因非主題原故，此處不再深入研究探討。
留言列表
開發手札 (2)
{{ article.title }}