R/newtonJ.R

In rootoned: Roots of one-dimensional functions in R-only code

##
##  n e w t o n . R  Newton Root finding
##


newtonRaphson <- function(fun, x0, dfun = NULL, ...,
	               maxiter = 100, tol = .Machine$double.eps^0.5) {
	# Newton method for finding function zeros
	if  (is.null(dfun)) {
		dfun <- function(x, ...) { h <- tol^(2/3)
			(fun(x+h, ...) - fun(x-h, ...)) / (2*h)
		}
	}
	x   <- x0
	fx  <- fun(x, ...)
	dfx <- dfun(x, ...)
	niter <- 0
	diff  <- tol + 1
	while (diff >= tol && niter <= maxiter) {
		niter <- niter + 1
        if (dfx == 0) {
            warning("Slope is zero: no further improvement possible.")
            break
        }
		diff  <- - fx/dfx
		x <- x + diff
		diff <- abs(diff)
		fx  <- fun(x, ...)
		dfx <- dfun(x, ...)
	}
	if (niter > maxiter) {
		warning("Maximum number of iterations 'maxiter' was reached.")
	}
	return(list(root=x, f.root=fx, niter=niter, estim.prec=diff))
}


newton <- newtonRaphson


halley <- function(fun, x0,
                   maxiter = 100, tol = .Machine$double.eps^0.5, ...) {
    f0 <- fun(x0,...)
    if (abs(f0) < tol^(3/2))
        return(list(root = x0, f.root = f0, maxiter = 0, estim.prec = 0))

    f1 <- fderiv(fun, x0, 1, ...)
    f2 <- fderiv(fun, x0, 2, ...)
    x1 <- x0 - 2*f0*f1 / (2*f1^2 - f0*f2)

    niter = 1
    while (abs(x1 - x0) > tol && niter < maxiter) {
        x0 <- x1
        f0 <- fun(x0, ...)
        f1 <- fderiv(fun, x0, 1, ...)
        f2 <- fderiv(fun, x0, 2, ...)
        x1 <- x0 - 2*f0*f1 / (2*f1^2 - f0*f2)
        niter <- niter + 1
    }
    return(list(root = x1, f.root = fun(x1, ...),
                iter = niter, estim.prec = abs(x1 - x0)))
}