圖形與說明
牛頓法以圖形概念主要是不停在取切線斜率。
考慮一方程式 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
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, func, eps, max_iterator);
printf("NewtonRootDiff func(%+.15e) = %+.15e\n", x, func(x));
x = NewtonRootXn(low, up, func, eps, max_iterator);
printf("NewtonRootXn func(%+.15e) = %+.15e\n", x, func(x));
x = NewtonRootDiffm(low, up, func, eps, 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 若用隨機方式產生時，
所找之解並不會每次找到都一樣。
留言列表
開發手札 (2)
{{ article.title }}