[C語言數值分析] 非線性方程式求解 - 定點法 Fix Point－Edison.X. Blog

[回目錄]

圖形與說明

假設一函式為 f(x) = 3x^3 + 2x^2 + 5x - 1，

要求其根，則令 f(x) = 0，即 3x^3 + 2x^2 + 5x - 1=0 進行求解。

而定點迴路法較特別，它是先將 f(x) = 0 ，化作 g(x) = x 之形式，

以此例而言，共有三種化法

(1) 移 x 項

3x^3 + 2x^2 + 5x - 1 = 0 ,

5x = -3x^3 - 2x^2 +1

x = (-3x^3 - 2x^2 +1) / 5

g(x) = (-3x^3 - 2x^2 +1) / 5

(2) 移 x^2 項

3x^3 + 2x^2 + 5x - 1 = 0

2x^2 = -3x^3 - 5x +1

x = sqrt [ (-3x^3 - 5x +1) /2 ]

g(x) = sqrt [ 0.5 * (-3x^3 - 5x +1) ]

(3) 移 x^3 項

3x^3 + 2x^2 + 5x - 1 = 0

3x^3 = -2x^2 - 5x +1

x = [( -2x^2 - 5x +1)/3] ^ (1/3)

g(x) = [( -2x^2 - 5x +1)/3] ^ (1/3)

將 f(x)=0 化作 x = g(x) 之後，只要求得 y=x 與 y=g(x) 之交點 x0，

便使得 f(x0) = 0 成立，即為 f(x) 其根。

從 f(x) = 0 化成 g(x) = x 後，必須先給一初始點 x0，假設圖型如下

接著計算 x1 = g(x0)

由於 (y=x) ，x1 代入 (y=x)，得到 x2，即 x2 = x1

再將 x2 代入 g(x), 得 x3，即 x3 = g(x2)

x3 代入 (y=x) 得 x4，即 x4 = x3，依此類推下去，

直到 abs(xn - xn-1) < eps 時終止。

再看一下剛剛迭代路徑

於是可得到一迭代式：

x1 = g(x0)

x0 = x1

這次程式碼例子，將以

f(x) = x^3 - 5.48*x^2 - 1.4883*x + 20.394828

為例，而 g(x) 則針對 x^2 移項做處理，即

g(x) = sqrt( (x*x*x - 1.4883*x+20.394828) /5.48);

虛擬碼

Algorithm FixPoint

E0 : 初始化最小誤差 EPS，初始點 x0

E1 : 計算 x = g(x0)

E2 : if abs(x0-x) < eps 成立，演算法終止，傳回 x0

E3 : x0 = x

E4 : goto E1

End Algorithm

C 語言程式碼

/*******************************************************************/
/*                                                                 */
/*     filename : FixPoint1.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 = g(x0)
| E2 : if abs(x0-x) < eps 成立，演算法終止，傳回x0
| E3 : x0 = x
| 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);
}

double gfunc(double x)
{
     double x2=x*x, x3=x2*x;
     return sqrt( (x*x*x - 1.4883*x+20.394828) /5.48);
}

double FixPoint1(double x0,              /* 初點 */
                 double (*gx)(double),   /* 適應函式*/
                 double eps,             /* 容許誤差*/
                 int max_itera)          /* 最大迭代*/
{
     double x = gx(x0);
     while(fabs(x-x0)>eps && max_itera){
           x0 = x;
           x = gx(x0);
           --max_itera;
     }
     return x;
}

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

     x0 = -2.0 ;
     x = FixPoint1(x0, gfunc, eps, max_itera);
     printf("\n >func(%+.15e) = %+.15e", x, func(x));

     x0 = +2.0 ;
     x = FixPoint1(x0, gfunc, eps, max_itera);
     printf("\n >func(%+.15e) = %+.15e", x, func(x));

     x0 = +4.0 ;
     x = FixPoint1(x0, gfunc, eps, max_itera);
     printf("\n >func(%+.15e) = %+.15e", x, func(x));

     x0 = 0.0 ; // test
     x = FixPoint1(x0, gfunc, eps, max_itera);
     printf("\n >func(%+.15e) = %+.15e", x, func(x));
     return 0;
}

執行結果

>func(+2.313837834730220e+000) = +9.258705802039913e-009
>func(+2.313837834658870e+000) = +1.002830884999639e-008
>func(+2.313837836693511e+000) = -1.191818910228903e-008
>func(+2.313837834508028e+000) = +1.165536644975873e-008

定點法收斂性

由程式結果發現，當 g(x) 以 x^2 項化開時，x0 分別給 -2，2，4，0，

最後都是收斂到同一解，原因乃在於，化開之 y= g(x) 與 y=x 兩條曲線，

在這區間中，圖型上只有一個交點，故所指向之解為同一個。

針對此題，若對 x^3 項、x 項化成 g(x) 時，也未必找得到解。

定點法要找到解有幾個事項要注意。

(1) 收斂充份條件

由 f(x) = 0 化作 g(x) = x ，於 [a, b] 尋找，必須確保 g(x) 可微分外，

「最好」再確保 | g'(x) | < 1，若 | g'(x) | < 1 成立，則必定收斂；

但若 | g'(x) | > 1 時，則「有可能收斂，有可能不收斂」，故稱充份條件。

以下圖形為例，最後只會發散，找不到 y=x 與 y=g(x) 之交點。