root finding - 4种演算法求函数的根
Numerical root finding methods use iteration, producing a sequence of numbers that hopefully converge towards a limits which is a root. This post only focuses four basic algorithms on root finding, and covers bisection method, fixed point method, Newton-Raphson method, and secant method.
The simplest root finding algorithms is the bisection method. It works when f is a continuous function and it requires previous knowledge of two initial gueeses, u and v, such that f(u) and f(v) have opposite signs. This method is reliable, but converges slowly. For detail, see 《用中值定理求根》.
Root finding can be reduced to the problem of finding fixed points of the function g(x) = c*f(x) +x, where c is a non-zero constant. It is clearly that f(a) = 0 if and only if g(a) = a. This is the so called fixed point algorithm.
fixedpoint <- function(fun, x0, tol=1e-07, niter=500){
## fixed-point algorithm to find x such that fun(x) == x
## assume that fun is a function of a single variable
## x0 is the initial guess at the fixed point
xold <- x0
xnew <- fun(xold)
for (i in 1:niter) {
xold <- xnew
xnew <- fun(xold)
if ( abs((xnew-xold)) < tol )
return(xnew)
}
stop("exceeded allowed number of iterations")
}
> f <- function(x) log(x) - exp(-x)
> gfun <- function(x) x - log(x) + exp(-x)
> fixedpoint(gfun, 2)
[1] 1.309800
> x=fixedpoint(gfun, 2)
> f(x)
[1] 3.260597e-09
The fixed point algorithm is not reliable, since it cannot guaranteed to
converge. Another disavantage of this method is that the speed is relatively slow.
Newtom-Raphson method converge more quickly than bisection method and fixed point method. It assumes the function f to have a continuous derivative. For detail, see 《Newton-Raphson Method估算函数的根》.
The secant method does not require the computation of a derivative, it only requires that the function f is continuous. The secant method is based on a linear approximation to the function f. The convergence properties of the secant method are similar to those of the Newton-Raphson method.
secant <- function(fun, x0, x1, tol=1e-07, niter=500){
for ( i in 1:niter ) {
x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
if (abs(fun(x2)) < tol)
return(x2)
x0 <- x1
x1 <- x2
}
stop("exceeded allowed number of iteractions")
}
> f <- function(x) log(x) - exp(-x)
> secant(f, x0=1, x1=2)
[1] 1.309800
往期精彩