[回目錄]


圖形與說明

 

牛頓法以圖形概念主要是不停在取切線斜率。

考慮一方程式 f(x) ,一開始必須找到初始點 a,然後代入 f(a)。

 

Newton-Raphson01.png

 

在 f(x) 圖形上,對 (a, f(a)) 做切線,與 x 軸相交於 b 點,且該切線斜率為 f'(a),

 

 Newton-Raphson02.png 

 

連接 (a, f(a) ), (b, 0) 二點,又已知斜率為 f'(a),可得一恆等式

f'(a) = ( 0 - f(a) ) / (b-a)

整理移項,得   b = a - f'(a) / f(a)

接著再將 b 對應到 f(b) 上去,得到點 (b, f(b))

Newton-Raphson03.png

 

對 (b, f(b)) 再做切線,切線於 x 軸於 c

 Newton-Raphson04.png  

有了 (b, f(b)), (c, 0) 兩點,斜率 f'(b),可得切線方程式

c = b - f(b) / f'(b)

再重覆上述動作。若將 a 視為 x0, b 視為 x1, c 視為 x2,

兩個直線公式便成為

x1 = x0 - f(x0) / f'(x0)

x2 = x1 - f(x1) / f'(x1)

依此類推

xk+1 = xk -  f(xk) / f'(xk)

一直迭代到 | xk-xk+1 | < eps 或 | f(xk) | < eps 為止。

 

使用條件與特性

 

1.  函數 f(x) 必須可微分。

2.  有幾種圖形沒辦法順利收斂,下圖即為其中一種。

Newton-Raphson05.png 

3. 若收斂點為重根,則收斂速度慢,假設已知為 m 個重根,

可將迭代公式改成 xk+1 = xk -  m * f(xk) / f'(xk)

缺點是,若重根數愈多,則找到的解誤差愈大。

但若未知重根數時,就變得麻煩了,

先設一函式 u(x) = f(x) / f'(x),迭代公式改成

xk+1 = xk - u(xk) / u'(xk)

實作上顯得較為麻煩。

 

 

虛擬碼


Algorithm NewtonRoot

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

    E1 :  x = x0

    E2 :  x0 = x - f(x) / f'(x)

    E3 :  if abs(x-x0) < eps,演算法結束,傳回 x

    E4 :  goto E1

End Algorithm

 

C 語言程式碼

 

/*******************************************************************/
/*                                                                 */
/*     filename : NewtonRoot.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
| E1 :  x  = x0
| E2 :  x0 = x - f(x) / f'(x)
| E3 :  if abs(x-x0) < eps,演算法結束,傳回x0
| E4 :  goto E1

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


#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);
}

// funcd(x) = func'(x)
double funcd(double x)
{
     double x2=x*x;
     return (3*x2-10.96*x-1.4883);
}

// -------------------------------------------------------
double NewtonRoot(double x0,              /*  
初點  */
                  double (*fx)(double),   /* 適應函式*/
                  double (*fd)(double),   /* 微分函式*/
                  double eps,             /* 容許誤差*/
                  int max_itera)          /* 最大迭代*/
{
     double x=x0;
     do{
           x0=x;
           x = x0 - fx(x0) / fd(x0);
     }while(fabs(x-x0)>eps);
     return x;

}

int main()
{
     const double eps=1E-9;
     const int max_iterator=100;
     double x0, x;

     x0 = -2.0;
     x = NewtonRoot(x0, func, funcd, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     x0 = 2.0;
     x = NewtonRoot(x0, func, funcd, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     x0 = 4.0 ;
     x = NewtonRoot(x0, func, funcd, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     x0 = -10.0 ; // test
     x = NewtonRoot(x0, func, funcd, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));
     return 0;
}

 

執行結果

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


改善額外之微分函式 fd(x)

 

上述之牛頓法必須提供一個微分函式才能順利完成,

在使用上並不非常方法,但可有一簡單方式改善此問題,

即利用差分方式取代掉,f'(x) = f(x1) - f(x2) / (x1 - x2)

於是初值便需要二個點,才能計算出第一個差別,程式碼約略如下

 

// -------------------------------------------------------
double NewtonRootDiff(double x0,              /*   初點  */
                      double x1,              /*   初點  */
                      double (*fx)(double),   /* 適應函式*/
                      double eps,             /* 容許誤差*/
                      int max_itera)          /* 最大迭代*/
{
     double y0= fx(x0), y1;
     double diff, x2;
     do{
           y1=fx(x1);
           diff = (y1 - y0) / (x1 - x0);
           x2 = x1 - y1 / diff;
           x0 = x1, y0 = y1, x1 = x2;
           --max_itera;
     }while(fabs(x1-x0)>eps && max_itera);
     return x2;
}

 

執行結果

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

 

 

改善收斂速度慢之問題

 

若遇到牛頓法收斂慢時, 針對普遍性問題大多都可用「下山法」方式完成,

其概念大致上如下

(1) 初始化 rate = 1.0,初始化 min_rate (user define)

(2)  xk+1 = xk -  rate * f(xk) / f'(xk) 

(3) | f(xk+1) | > | f(xk) | 成立時至 (4);不成立時至 (2)

(4) rate > min_rate 時 rate = 0.5 * rate,回到 (2)

 

實作上細節必須細思,但實際上有懶人法,rate 直接用亂數生成,

在合理生成範圍內,往往會有意想不到的效果 (往往比下山法還優)。

 

C 語言程式(改善與測試)


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


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

double func(double x)
{
     double t = x-1.235;
     return pow(t, 5.0);
}


// -------------------------------------------------------
double NewtonRootDiff(double x0,              /*  
初點  */
                      double x1,              /*   初點  */
                      double (*fx)(double),   /* 適應函式*/
                      double eps,             /* 容許誤差*/
                      int max_itera)          /* 最大迭代*/
{
     double y0= fx(x0), y1;
     double diff, x2;
     int i=0;
     do{
           y1=fx(x1);
           diff = (y1 - y0) / (x1 - x0);
           x2 = x1 - y1 / diff;
           x0 = x1, y0 = y1, x1 = x2;
           ++i;
     }while(fabs(x1-x0)>eps && i<max_itera);
     printf("> times: %4d , ", i);
     return x2;
}

// -------------------------------------------------------
double NewtonRootDiffm(double x0,              /*  
初點  */
                       double x1,              /*   初點  */
                       double (*fx)(double),   /* 適應函式*/
                       double eps,             /* 容許誤差*/
                       int max_itera)          /* 最大迭代*/
{
     const int m=5; //
若已知重根數5 個
     double y0= fx(x0), y1;
     double diff, x2;
     int i=0;
     do{
           y1=fx(x1);
           diff = (y1 - y0) / (x1 - x0);
           x2 = x1 - m * y1 / diff;
           x0 = x1, y0 = y1, x1 = x2;
           ++i;
     }while(fabs(x1-x0)>eps && i<max_itera);
     printf("> times: %4d , ", i);
     return x2;
}

// -------------------------------------------------------
double NewtonRootXn(double x0,              /*  
初點  */
                    double x1,              /*   初點  */
                    double (*fx)(double),   /* 適應函式*/
                    double eps,             /* 容許誤差*/
                    int max_itera)          /* 最大迭代*/
{
     double y0= fx(x0), y1;
     double diff, x2, delta;
     double rate=1.0;
     int i=0;

     srand(time(NULL));
     do{
           delta = y1 = fx(x1);
           // if(fabs(y1) < fabs(y0)) rate=1.0;
           // else if(rate > 5E-5) rate*=0.5;
           rate = 5.0 * rand() / RAND_MAX; // more faster

           diff = (y1 - y0) / (x1 - x0);
           x2 = x1 - rate * y1 / diff;         
           x0 = x1, y0 = y1, x1 = x2;
           ++i;
     }while(fabs(delta)>eps && i<max_itera);
     printf("> times: %4d , ", i);
     return x2;
}

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

     low = -10.0, up = low + 1.0;
     x = NewtonRootDiff(low, up, funceps, max_iterator);
     printf("NewtonRootDiff  func(%+.15e) = %+.15e\n", x, func(x));
     x = NewtonRootXn(low, up, funceps, max_iterator);
     printf("NewtonRootXn    func(%+.15e) = %+.15e\n", x, func(x));
     x = NewtonRootDiffm(low, up, funceps, max_iterator);
     printf("NewtonRootDiffm func(%+.15e) = %+.15e\n", x, func(x));
     return 0;
}

執行結果

> times:  138 , NewtonRootDiff  func(+1.234999994620308e+000) = -4.505958212051828e-042
> times:   15 , NewtonRootXn    func(+1.219286789151360e+000) = -9.579099640434134e-010
> times:    4 , NewtonRootDiffm func(+1.220320543356181e+000) = -6.816318936216687e-010

 

雖 NewtonRootDiffm 在此例為最快,但若該方程式非「完全重根」(還有其他解時),

速度並不會比下山法來得快。另下山法其 rate 若用隨機方式產生時,

所找之解並不會每次找到都一樣。

 

[回目錄]

arrow
arrow
    全站熱搜

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