[回目錄]

 

大概只有 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;
}

 

[回目錄]

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()