範例EX2-3 牛頓法求 x^3 +4x^2 -10 =0 的根

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

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

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

  

連接 (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))

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

   

有了 (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.  有幾種圖形沒辦法順利收斂，下圖即為其中一種。

 

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

程式碼 

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

| 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 + 4*x2 -10);

}

// funcd(x) = func'(x)

double funcd(double x)

{

     double x2=x*x;

     return (3*x2 + 8*x);

}

// -------------------------------------------------------

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(-nan) = -nan                                                                                                                                           

> func(+1.365230013414097e+00) = +0.000000000000000e+00                                                                                                       

> func(+1.365230013414097e+00) = +0.000000000000000e+00                                                                                                       

> func(+1.365230013414097e+00) = +0.000000000000000e+00                                                                                                       

                                                                                                                                                              

...Program finished with exit code 0                                                                                                                          

Press ENTER to exit console. 

源自於

http://edisonx.pixnet.net/blog/post/35763002-%5Bc%E8%AA%9E%E8%A8%80%E6%95%B8%E5%80%BC%E5%88%86%E6%9E%90%5D-%E9%9D%9E%E7%B7%9A%E6%80%A7%E6%96%B9%E7%A8%8B%E5%BC%8F%E6%B1%82%E8%A7%A3---%E7%89%9B%E9%A0%93%E6%B3%95-n