Root finding
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Numerical root finding methods use iteration, producing a sequence of numbers that hopefully converge towards a limits which is a root. In this post, only focus four basic algorithm 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 http://ygc.cwsurf.de/2008/11/25/bisect-to-solve-equation/ .
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 fixed point method 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 http://ygc.cwsurf.de/2007/06/02/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
Related Posts
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.