[回目錄]

 

圖形與說明

 

[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) 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

 

[回目錄]

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


留言列表 (2)

發表留言
  • 訪客
  • 處理意外方式有寫錯嗎?
    若y1*y0<0成立,則x0=x2,y0=y2;
    或y2*y1<0成立
  • 沒錯啊 @@

    edisonx 於 2014/02/03 11:41 回覆

  • 訪客
  • 抱歉問了蠢問題,後來想通了
    若解不在[x0,x2],而[x0,x1]間有解,則解在[x2,x1]

您尚未登入,將以訪客身份留言。亦可以上方服務帳號登入留言

請輸入暱稱 ( 最多顯示 6 個中文字元 )

請輸入標題 ( 最多顯示 9 個中文字元 )

請輸入內容 ( 最多 140 個中文字元 )

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼