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

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'{(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