Edison.X. Blog

[回目錄]

大概只有 BiSector 方式可以硬找出來，同時也有點依附於 compiler 特性。

這問題是來自於某個學生，因初值那裡不好找，不知道該如何找下去，

目前所有方法裡面，就只有二分法對任何函數都必收斂，

於是採二分法做 Force。

關鍵在於初值時，low bound 與 up bound 設定，

double low bound = -FLT_MAX

double up  bound =   FLT_MAX

接著再進行二分法，下面就是關鍵了。

在二分法過程中，代入 f(x) 時，很可能會造成 overflow (-1.#IND)，

但好笑的是，這種 over flow 情況若不管它，

繼續代下去，結果還是收斂了 (oh~我真不懂 compiler 做了什麼事)，

再把 low / up bound 換更誇張的 ， -DBL_MAX / -DBL_MIN，

結果也收斂了。

-FLT_MAX / FLT_MAX ，最大計算次數應為 log2(2 * FLT_MAX)，

大概 128 次左右；而 -FLT_MAX / FLT_MAX，大概 1026 次左右。

這種方式不確定是不是每個含有實數解之函數都適用，

但這方式實在認為很妙，特別是莫名奇妙還弄不懂為什麼 OV 後，

答案還是會慢慢浮出..

以下程式碼僅供參考。

/*******************************************************************/
/*                                                                 */
/*     filename : FindRoot.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 <float.h>
#include <math.h>

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

double BiSect(double eps,
              double (*fx)(double))
{
     double low = -DBL_MAX;
     double up  =  DBL_MAX;
     double mid = 0.0;
     double ylow=fx(low), yup=fx(up), ymid=fx(mid);
     unsigned long long times=0;
     while(fabs(ymid) > eps){
           if(ymid * yup < 0.0) ylow = ymid, low=mid;
           else yup = ymid, up=mid;
           mid  = 0.5 * (up+low);
           ymid = fx(mid);
           ++times;
     }
     printf("\ntimes : %u\n", times);
     return mid;
}

int main()
{
     const double eps=0.0;
     double root, val;

     root = BiSect(eps, func), val  = func(root);
     printf("root : %+.15e\n", root);
     printf("func : %+.15e\n", val);
     return 0;
}

[回目錄]