求根

　
以下程式輸出為：

0.962398
-0.040659
預測值已跨出區間!
0.000000

前兩者在 [0,1] 「附近」找過零點，並不保證不越區。
此函式在那附近正好有兩個： 0.962398 與 -0.040659
Secant 根據給出的初始域，往 -0.040659 靠。
至於 Newton-Raphson 就更不可靠了，因發生越區而傳回 0。
若將初始域放大， Newton-Raphson 便可找到根值。

typedef double DF;              //double float
typedef DF (*FP)(DF);

DF abs (DF f) {if (f<0) f=-f; return f;}

//以精確度 accuracy, 在 loop 次迭代內求 f 在 [x1,x2] 中的根,
//注意! 參數需滿足 x2>x1, 且 [x1,x2] 中至少包含一根

DF false_pos (FP f, DF x1, DF x2, DF accuracy, int loop=64)
{
    DF acc, x,y, y1 = f(x1), y2 = f(x2);

    for (int i=0; i<loop; ++i) {
        y = f (x = x1 + (x2-x1) * y1/(y1-y2));
        if (y < 0) acc=x1-x, x1=x, y1=y;
        else       acc=x2-x, x2=x, y2=y;
        if (abs(acc) < accuracy || 0==y) return x;
    }
    puts ("root not be found!");
    return 0;
}

DF secant (FP f, DF x1, DF x2, DF accuracy, int loop=64)
{
    DF acc, dx, x = x2, y1 = f(x1), y = f(x2);
    if (abs(y1) < abs(y)) x=x1, x1=x2, dx=y1, y1=y, y=dx;

    for (int i=0; i<loop; ++i) {        
        dx = (x1-x) * y/(y-y1);
        x1 = x;
        y1 = y;
        x += dx;
        y = f (x);
        if (abs(dx) < accuracy || 0==y) return x;
    }
    puts ("root not be found!");
    return 0;
}

//df 為 f 的導函數
DF newton (FP f, FP df, DF x1, DF x2, DF accuracy, int loop=64)
{
    DF acc, dx, x = (x1+x2)*.5;   
    for (int i=0; i<loop; ++i) {        
        x -= dx = f(x) / df(x);
        if ((x1-x)*(x-x2) < 0) {puts("\n預測值已跨出區間!"); return 0;}
        if (abs(dx) < accuracy) return x;
    }
    puts ("root not be found!");
    return 0;
}

DF f (DF x) { return (((230*x+18)*x+9)*x-221)*x-9; }
DF df(DF x) { return ((230*x+18)*x+9)*x-221; }

int main()
{
    printf ("\n%f", false_pos (f, 0,1, 1e-6));
    printf ("\n%f", secant (f, 0,1, 1e-6));
    printf ("\n%f", newton (f, df, 0,1, 1e-6));
    return 0;
}