Newton-Raphson Method

曲线f(x)有根c,取曲线上一点$(x_1,f(x_1))$, 过此点的切线交x轴$x_2$,过曲线上$(x_2,f(x_2))$的切线交x轴$x_3$，如此反复得到一个序列 $x_1,x_2,\cdot \cdot \cdot,x_n$ 逼近c值.

过$(x_n,f(x_n))$的切线方程为 $y-f(x_n) = f’(x_n)\,(x-xn)$,假设此方程与x轴的交点为$x{n+1}$, 即有: $0 - f(x_n) = f’(x_n)(x_n+1 - xn)$， 即$x{n+1} = x_n - \frac{f(x_n)}{f’(x_n)}$ <Eq. 1>.

下面利用此法来求一个数的开方。 $f(x) = x^2 - a$ 有根$\sqrt{a}$, 由$f’(x_n) = 2xn$, 代入式<Eq. 1>可得$x{n+1} = (x_n + a/x_n)/2$; 当i -> INF 时, $x_i$ -> $\sqrt{a}$;

C implementation

#include

int main(void){
    int a,error;
    double x0,x1 = 1;
    do {
        printf("Input a positive integer: ");
        scanf("%d",&a);
        if (error = (a <=0))
            printf("\nERROR: Do it again!\n\n");
        }
    while (error);

 while (x0 != x1) {
        x0 = x1; /* save the current value of x1 */
        x1 = 0.5 * (x1 + a / x1); /* compute a new value of x1 */
    }
    printf("%lf\n",x1);
    return 0;
}

R implemantation

2010-01-11 用R来实现一下，不单是求开方，估算函数的根。

newton <- function(fun, x0, tol=1e-7, niter=100) { 
    fun.list = as.list(fun)
    var <- names(fun.list[1])
    fun.exp = fun.list[[2]] 
    diff.fun = D(fun.exp, var) 
    df = list(x=0, diff.fun) 
    df = as.function(df) 
    for (i in 1:niter) { 
        x = x0 - fun(x0)/df(x0) 
        if (abs(fun(x)) < tol) 
            return(x) 
        x0 = x      
    } 
    stop("exceeded allowed number of iterations") 
}

> f = function(x) x^2 – 5
> newton(f, 4)
[1] 2.236068
> g = function(x) x^3 – 5
> newton(g, 4)
[1] 1.709976