Root finding

December 4, 2010
By

(This article was first published on YGC » R, and kindly contributed to R-bloggers)

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.name/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.name/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

To leave a comment for the author, please follow the link and comment on his blog: YGC » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , ,

Comments are closed.