Edison.X. Blog

[回目錄]

說明

筆者目前見過求一解最好的方法便是穆勒法 ，

它是在 f(x) 一次取三個點，(a, f(a) )，(b, f(b) )，(c , f(c)) ，

再以二次函式方式去逼近。

穆勒法目前在數值分析的書上也較為罕見，

且過程甚為攏長，光是初始就要進行一大段敘述。

雖方法非常不容易，一次迭代就要調用三次 function，

但收斂速度快，且初值可設的範圍非常廣，解品質甚佳，搜尋性強。

此處不進行圖解說明，只附程式碼作為參考。

C 語言程式碼

/*******************************************************************/
/*                                                                 */
/*     filename : Muller.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 <float.h>
double func(double x)
{
     double x2=x*x, x3=x2*x;
     return (x3 - 5.48*x2 -  1.4883*x + 20.394828);
}

double Muller(double x1,               /*   初點  */
              double x2,               /*   初點  */
              double x3,               /*   初點  */
              double (*fx)(double),    /* 適應函式*/
              double eps,              /* 容許誤差*/
              int max_itera)           /* 最大迭代*/
{
     double x4_1,x4_2,h1,h2,h3_1,h3_2,h4,D,d1,d2,a1,a2,a0;
     double fx1,fx2,fx3, tmp, y1, y2, y3;

     fx1 = fx(x1), fx2=fx(x2), a0=fx3=fx(x3);
     h1 = x1-x3,   h2 = x2-x3;
     d1 = fx1-fx3, d2 = fx2-fx3;
     D = h1*h2*(h1-h2);

     a1 = (d2*h1*h1 - d1*h2*h2)/D;
     a2 = (d1*h2 - d2*h1)/D;

     tmp = sqrt(fabs(a1*a1 - (4*a2*a0)));
     h3_1 = -((2*a0)/(a1 + tmp));
     h3_2 = -((2*a0)/(a1 - tmp));

     if(a1+tmp >a1-tmp) h4 = h3_1;
     else h4 = h3_2;

     x4_1 = x3 + h4;
     x1=x2, x2=x3, x3=x4_1;

     do{
           fx1 = fx(x1), fx2 = fx(x2);
           fx3 = a0 = fx(x3);

           h1 = x1-x3, h2 = x2-x3;

           d1 = fx1-fx3, d2 = fx2-fx3;
           D = h1*h2*(h1-h2);

           a1 = (d2*h1*h1 - d1*h2*h2)/D;
           a2 = (d1*h2 - d2*h1)/D;

           tmp = sqrt(fabs(a1*a1 - (4*a2*a0)));
           h3_1 = -((2*a0)/(a1 + tmp));
           h3_2 = -((2*a0)/(a1 - tmp));

           if(a1+tmp >a1-tmp) h4 = h3_1;
           else h4 = h3_2;
           x4_2 = x3 + h4;

           if(fabs(x4_1 - x4_2) < eps) max_itera=0;
           else
           {
                x4_1=x4_2;
                x1=x2, x2=x3, x3=x4_1;
                --max_itera;
           }
          
     }while(max_itera);
     return x4_2;
}

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

     a=-2.0, b=-1.5, c=-1.0;
     x = Muller(a, b, c , func, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     a= 0.0, b= 2.5, c= 2.0;
     x = Muller(a, b, c , func, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     a= 4.0, b= 4.5, c=5.0;
     x = Muller(a, b, c , func, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));

     a=-2E10, b=0.0, c=2E10;
     x = Muller(a, b, c , func, eps, max_iterator);
     printf("\n> func(%+.15e) = %+.15e", x, func(x));
     return 0;
}

[回目錄]

執行結果

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

結果有二點值得注意

1. 即使初始解給定了 [a, b]，最後找到的解未必落於該區間。

2. 可搜尋範圍甚廣，即使區間定為 [-2E10, 2E10]，仍可找到解，幾乎是 int 之解空間範圍。