Root finding

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

  1. Nice post :)

    Here's a crude approach using R which I've found to be quite handy when precision and speed aren't a big concern.

    Reply

    ygc China Unknow Browser Unknow Os Reply:

    Thanks for sharing.
    It's much straightforward to do it by:

    > fun = function(x) x/(1+x) -(5-x)/5
    > x = seq(0, 5, length=5000)
    > which.min(abs(fun(x)))
    [1] 1792
    > fun(x[1792])
    [1] 2.312344e-05

    You need not to split the function to f=x/(1+x) and g=(5-x)/5. Many functions, such as cos(x), cannot be splitted.

    Reply

    Paul United States Unknow Browser Unknow Os Reply:

    Very true, using which.min() would be more efficient. With some tinkering, you can always deal with multiple roots by do something like...

    > x=seq(0,15,length=5000)
    > curve(cos(x),0,15);
    > points(x[which((abs(cos(x))<2e-3))],
    cos(x[which((abs(cos(x))

    :smile:

    Reply

    ygc China Unknow Browser Unknow Os Reply:

    Absolutely right.
    Typo error occur in last line.
    should be:
    points(x[which((abs(cos(x))<2e-3))],cos(x[which((abs(cos(x))<2e-3))]))

    Reply

Leave a Comment


NOTE - You can use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>