# 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