The Modified Newton-Raphson Method

First, we define the function we want to use.

fn <- function(x){x^4-8.6*x^3-35.51*x^2+464.4*x-998.46} # other function: exp(-x)-x
fn(4)
## [1] -3.42

Second, get the derivative.

fne<- expression(x^4-8.6*x^3-35.51*x^2+464.4*x-998.46)
fn1<- D(fne, 'x')
# this will print the derivative 
fn1 #3 * x^2 - 2 * (2 * x) + 1
## 4 * x^3 - 8.6 * (3 * x^2) - 35.51 * (2 * x) + 464.4

#Third, define new function of the first derivative

fn_1<- function(x){4 * x^3 - 8.6 * (3 * x^2) - 35.51 * (2 * x) + 464.4}
fn_1(4)
## [1] 23.52

And define new function of the second derivative

fne2<- expression(4 * x^3 - 8.6 * (3 * x^2) - 35.51 * (2 * x) + 464.4)
fn2<- D(fne2, 'x')
fn2
## 4 * (3 * x^2) - 8.6 * (3 * (2 * x)) - 35.51 * 2
fn_2<- function(x){4 * (3 * x^2) - 8.6 * (3 * (2 * x)) - 35.51 * 2}
fn_2(4)
## [1] -85.42

To graph the function we have to see the initial point

curve(fn, from=2, to=10, xlab="x", ylab="y")
abline(h=0,col="red",lwd=1.5)

NOTE: you can change the range of x (from, to) _____________________

Eventually we define the function of Modified Newton-Raphson

library(rlist)
## Warning: package 'rlist' was built under R version 4.1.3
ModNewRaph<- function(xi,es=0.05,j=10){ 
  i=0 #intial value for i
  ea<-1.1*es
  x0=xi
  newlist <- list()
  
  while (ea>es & i<j) { 
    x1<- x0-(fn_1(x0)*fn(x0)) / (fn_1(x0)^2-fn(x0)*fn_2(x0))
    
    ea<-abs((x1-x0)) #/x1 -> divide by x1 if you want relative error
      looplist <- list(i,x0,fn(x0),fn_1(x0),fn_2(x0),x1,ea)
      newlist <- list.append(newlist, looplist)
    
    x0=x1
    i=i+1
  }
  df <-  as.data.frame(do.call(rbind, newlist))
  list_name <- list("iteration",
                    "x(i)",
                    "f(xi)",
                    "f'(xi)",
                    "f''(xi)",
                    "x(i+1)",
                    "Error")
  names(df) <- c(list_name)
#  View(df)
  print(df)
    cat("x=",x0 ,"\n") #
}
ModNewRaph(4,0.005,5)
##   iteration     x(i)         f(xi)       f'(xi)   f''(xi)   x(i+1)        Error
## 1         0        4         -3.42        23.52    -85.42 4.308129    0.3081294
## 2         1 4.308129  -0.002342102   -0.5756394 -70.59973 4.300008  0.008121288
## 3         2 4.300008 -2.309548e-09 -0.000572742 -71.01958      4.3 8.065196e-06
## x= 4.3
## 
##  the following is to calculat *m*: 
## 
  fne <-
    expression(x ^ 4 - 8.6 * x ^ 3 - 35.51 * x ^ 2 + 464.4 * x - 998.46)
get_m <-function(x){
  fn1 <- fne
  a <- 0
  m<-0
  while (a == 0) {
    fn1 <- D(fn1, 'x')
    print(fn1)
    a <- eval(fn1)
    cat("the derivative= ",a,"\n")
    m<- m+1
    print(m)
  }
  
  #cat("x=",x0 ,"\n") #
  cat("m=",m ,"\n") #
}
get_m(4.3)
## 4 * x^3 - 8.6 * (3 * x^2) - 35.51 * (2 * x) + 464.4
## the derivative=  0 
## [1] 1
## 4 * (3 * x^2) - 8.6 * (3 * (2 * x)) - 35.51 * 2
## the derivative=  -71.02 
## [1] 2
## m= 2

let’s graph the function again

 x_1 = 4.3
 curve(
   fn,
   from = 3,
   to = 8,
   xlab = "x",
   ylab = "y"
 )
 abline(h = 0, col = "red", lwd = 1.5)
 abline(v = x_1, col = "red" , lwd = 1.5)
 text(x = 5, y = 15, labels = "(4.3,0)")

rmarkdown::render(“Modified Newten.R”)