Newton Raphson iterative method based on R language (for binary differentiable function)

Newton-Rapson Method

Newton Raphson method is an iterative method based on the initial value guess of the root. The function used in this method is the original function and the derivative of the original function. If it succeeds, it usually converges quickly, but it may fail like other root seeking methods, which should be noted. (because Newton's method does not always converge, its convergence theory acts on local convergence.)

Algorithmic thinking

This method tries to solve f (x) = 0 by iterative method
xn+1=xn−f(xn)f′(xn) x_{n+1} = x_{n} - \frac{f{(x_n)}}{f&#x27;{(x_n)}} xn+1​=xn​−f′(xn​)f(xn​)​

Geometric interpretation:

R code

You can use deriv()

#Write the function newton.rephson() , x0 is the initial x value,
#fun is the f(x), precision is the estimated accuracy
newton.raphson <- function(x1, x2, fun, precision, max.loop)
{
#Compute derivatives of simple expressions, symbolically and algorithmically
dy <- deriv(fun, c("x1", "x2"), hessian = TRUE)
x1 <- x1
x2 <- x2

#Build a vector to storage each iteration x value
estimate_x1 <- vector(mode = "numeric", length = 0)
estimate_x1[1] <- x1
estimate_x2 <- vector(mode = "numeric", length = 0)
estimate_x2[1] <- x2

#Evaluate an R expression in a specified environment
d <- eval(dy)
i <- 1
#Establish a loop until the estimated accuracy reaches the specified level
while((abs(d) > precision) && (i < max.loop))
{
i <- i+1
d <- eval(dy)
hessian <- matrix(attributes(d)$hessian, byrow = F, 2, 2) x1 <- x1 - solve(hessian, t(attributes(d)$gradient))[1,]
x2 <- x2 - solve(hessian, t(attributes(d)\$gradient))[2,]
estimate_x1[i] <- x1
estimate_x2[i] <- x2
}
#Return each iteration x value, the last one is the number closest to the zero point of the f(x)
return(list(Iteration_x1 = estimate_x1, Iteration_x2 = estimate_x2))
}



You can also use D()

newton.raphson <- function(f, init.value, df = NULL, tol = 1e-5, maxiter = 1000) {
if (is.null(df)) {
df <- D(f, 'x')
}
niter <- 0
diff <- tol + 1
x <- init.value
while (diff >= tol && niter <= maxiter) {
niter <- niter + 1
fx <- eval(f)
dfx <- eval(df)
if (dfx == 0) {
warning("Slope is zero: no further improvement possible.")
break
}
diff <- -fx/dfx
x <- x + diff
diff <- abs(diff)
}


Note that the selection of starting value is very important. The wrong starting value may lead to the failure of iteration method, and try to make your initial value close to your estimated root

Tags: R Language

Posted on Sun, 01 Dec 2019 15:29:02 -0800 by samoht