Edison.X. Blog

[回目錄]

圖形與說明

這解法的名稱多到讓人混亂，其實指的是一樣的東西。

(1) Secant，中文稱 割線法，另也有人稱正割法。

(2) 概念主要使用 拉格南奇 插值法(Lagrange's Interpolation)，故有些書稱此求根法為 插值法 。

(3) 若有書中另外再提及 假位法 / 錯位法 (False Position) ，方法與此法幾乎一模一樣，

     仔細比對迭代式，發現根本就是 Secant 割線法，不同的只有在換值部份而已。

圖形上不斷取兩點連線，與牛頓法相似，

比牛頓法好的地方在不用設微分函式；比牛頓法差之地方在初值要設二個。

x 軸上任意給二點 a, b ，對應至 f(x) 函式上，得二點為 (a, f(a) ), (b, f(b) )

將對應之二點做連線動作，得一直線 L，

根據 拉格南奇 二點插值 公式，割線 L 可得下列恆等式

y = f(b) * (x - b) / (a - b) + f(a) * (x - a) / (b - a)

再令此割線交 x 軸於 (x1, 0)

令 y =0 代入割線方程式，

y = 0 = f(b) * (x - b) / (a - b) + f(a) * (x - a) / (b - a)

整理可得

x1 = f(a) *b / f(a) - f(b) + f(b) *a / (fb) - f(a)

把 原本 a 點擦掉, 只看 b, x1

再將 (b, f(b)) , (x1, f(x1)) 相連

重覆上述步驟，可得 

x2 = f(b) *x1 / f(b) - f(x1) + f(x1) *b / f(x1) - f(b)

若將上述之 a, b, x1, x2 對照成 x0, x1, x2, x3，

可得一迭代公式

xi+2 = f(xi) * xi+1 / (f(xi) - f(xi+1)) + f(xi+1) *xi / ( f(xi+1) - f(xi) )

另也有一迭代式

xi+2 = xi+1  - ( f(xi+1) * (xi+1 - xi) ) / ( f(xi+1) - f(xi) )

一直到 (xi+2 - xi+1) < eps ，或 ( f(xi+1) - f(xi) ) < eps

雖迭代式二個化開來長得一樣，不過第二式所基本運算較少，筆者採用。

再清楚一點的迭代式： x2 = x1 -  y1 * (x1 - x0) / (y1 - y0)

注意，由於分母為 (y1 - y0) ，故這裡絕不能為 0 ，必設為收斂條件

虛擬碼

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

C 語言程式碼

/*******************************************************************/
/*                                                                 */
/*     filename : Secant.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，初始點x0, 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) or abs(y1-y0), 演算法結束, 傳回x2
| E5 :  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 Secant(double x0,              /*  初點0  */
              double x1,              /*  初點2  */
              double (*fx)(double),   /* 適應函式*/
              double eps,             /* 容許誤差*/
              int max_itera)          /* 最大迭代*/
{
     double y0=func(x0), y1=func(x1), x2;
     double delta=0.0;
     do{
           x2 = x1 - y1 * (x1-x0) / (y1-y0);
           delta = x2-x1;
           y0 = y1, x0 = x1, x1 = x2;
           y1 = fx(x2);
           --max_itera;
     }while(fabs(delta) > eps && fabs(y1-y0) && max_itera);
     return x2;

}

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

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

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

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

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

執行結果

 >func(-1.781503813804761e+000) = +3.552713678800501e-015
 >func(+2.313837835588586e+000) = +0.000000000000000e+000
 >func(+4.947665978216175e+000) = +3.552713678800501e-015
 >func(+4.947665978216176e+000) = +3.197442310920451e-014

[回目錄]