[C語言數值分析] 非線性方程式求解 - 假位法 False Position－Edison.X. Blog

[回目錄]

圖形與說明

在 [C語言數值分析] 非線性方程式求解 - 割線法 Secant 這篇文章我曾提到，

割線法和假位法方法幾乎一模一樣，他們使用的迭代公式都是

x2 = x1 - y1 * (x1 - x0) / (y1 - y0)

再看一次割線法虛擬碼

Algorithm SecantRoot

E0 : 初始化最小誤差 EPS，初始點 xo, x1

E1 : 算 y0 = f(x0), y1 = f(x1)

E2 : x2 = x1 - y1*(x1-x0) / (y1-y0)， delta = x2-x1

E3 : x0 = x1, x1 = x2, y0=y1, y1=f(x1)

E4 : if abs(delta) < eps 或 abs(y1-y0) < eps，演算法結束，傳回 x2

E5 : goto E2

End Algorithm

關鍵只在於 E3 取代部份不同，假位法乃採用了割線法之迭代方式，

而取代方式乃採用了二分法，即假設 y2 = f(x2)。

(1) 處理意外方式

若 y2 * y0 < 0 成立，則 x1 = x2，y1 = y2；

若 y1 * y0 < 0 成立，則 x0 = x2，y0 = y2；

若上述都不成立，則於 [x0, x1] 間沒有解。

(2) 不處理意外方式

若 y2 * y0 < 0 成立，則 x1 = x2，y1 = y2；

反之，則 x0 = x2，y0 = y2；

此處不採取意外處理。

這種方式的特點在於，保證搜尋空間必落在 [x0, x1] 之間，

然而割線法只要 "接近" 解時，便有機會找到其解，

同時假位法之收斂性與解品質，並沒割線法來得好，

故現之書籍，大多只提割線法，避開提起假位法。

虛擬碼

Algorithm SecantRoot

E0 : 初始化最小誤差 EPS，初始點 xo, x1

E1 : 算 y0 = f(x0), y1 = f(x1)

E2 : x2 = x1 - y1*(x1-x0) / (y1-y0)，y2 = f(x2)

E3 : if (y2 * y0) < 0 成立，delta = x2-x1, x1 = x2, y1=y2, goto E5

E4 : delta = x2 - x0, x0 = x2, y0 = y2

E5 : if abs(delta) < eps 或 abs(y1-y0) < eps，演算法結束，傳回 x2

E6 : goto E2

End Algorithm

C 語言程式碼

/*******************************************************************/
/*                                                                 */
/*     filename : FalsePosition.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.             */
/*                                                                 */
/*******************************************************************/

/*----------------------------------------------------------------*\
|
| E0 : 初始化最小誤差EPS，初始點xo, x1
| E1 : 算y0 = f(x0), y1 = f(x1)
| E2 : x2 = x1 - y1*(x1-x0) / (y1-y0)，y2 = f(x2)
| E3 : if (y2 * y0) < 0 成立，delta = x2-x1, x1 = x2, y1=y2；goto E5
| E4 : delta = x2 - x0, x0 = x2, y0 = y2
| E5 : if abs(delta) < eps 或abs(y1-y0) < eps，演算法結束，傳回x2
| E6 : goto E2
|
\*----------------------------------------------------------------*/

#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 FalsePosition(double x0,              /* 初點 */
                     double x1,              /* 初點 */
                     double (*fx)(double),   /* 適應函式*/
                     double eps,             /* 容許誤差*/
                     int max_itera)          /* 最大迭代*/
{
     double y0=func(x0), y1=func(x1);
     double x2, y2;
     double delta=0.0;
     do{
           x2 = x1 - y1 * (x1-x0) / (y1-y0);
           y2 = func(x2);
           if(y2 * y0 < 0 ) delta = y2-y1, y1 = y2, x1=x2;
           else delta = y2-y0, y0 = y2, x0 = x2;
           --max_itera;
     }while(fabs(delta) > eps && max_itera);
     return x2;
}

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

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

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

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

     low = -10.0 ; // test
     x = FalsePosition(low, up, func, eps, max_itera);
     printf("\n >func(%+.15e) = %+.15e", x, func(x));
     return 0;
}

執行結果

>func(-1.781503813803826e+000) = +2.577493773969763e-011
>func(+2.313837835590115e+000) = -1.649169689699193e-011
>func(+4.947665978192316e+000) = -4.228546401918720e-010
>func(+4.947665978435647e+000) = +3.889812916213487e-009

[回目錄]