Newton-Raphson Method
曲线f(x)有根c,取曲线上一点(x1,f(x1)),过此点的切线交x轴x2,过曲线上(x2,f(x2))的切线交x轴x3，如此反复得到一个序列x1,x2,...,xn逼近c值,过(xn,f(xn))的切线方程为y-f(xn) = f'(xn)(x-xn),假设此方程与x轴的交点为xn+1,即有0 - f(xn) = f'(xn)(xn+1 - xn)，即xn+1 = xn - f(xn)/f'(xn) <1>
下面利用此法来求一个数的开方。
f(x) = x^2 - a$有根\sqrt{a}$,
由f'(xn) = 2xn, 代入式<1>可得xn+1 = (xn + a/xn)/2;
当i -> inf 时, xi -> sqrt(a);
下面用C语言来实现这一算法。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | #include<stdio.h> 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; } |
2010-01-11
用R来实现一下，不单是求开方，估算函数的根。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | 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
- pingback on November 25, 2008 at 12:54 pm
Scala, 3行代码，5秒搞定：
scala> def lenOfPa...
3 Comments.