Introduction

Most numerical root-finding methods use iteration, producing a sequence of numbers that hopefully converge towards the root as a limit.

In this study, I am going to implement three following root finding methods. The goal of finding root is to find the local minimum or maximum points of the function. Setting zero of the first derivative of the continuous function returns the local/global maximum or minimum point(s).

  1. Bisection Method: In mathematics, the bisection method is a root-finding method that applies to any continuous functions for which one knows two values with opposite signs. The method consists of repeatedly bisecting the interval defined by these values and then selecting the subinterval in which the function changes sign, and therefore must contain a root.

  2. The Newton-Raphson method (also known as Newton’s method) is a way to quickly find a good approximation for the root of a real-valued function f(x) = 0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

  3. In numerical analysis, the secant method is a root-finding algorithm that uses a succession of roots of secant lines to better approximate a root of a function f. The secant method can be thought of as a finite-difference approximation of Newton’s method.

Objective

Use R function optimize to find the maximum of g(x) = log(x)/(x+1), Find the maximum of g over the interval [1; 5] with an accuracy of 0:001 (i.e., " = 0.0001).

Use the following optimization techniques and compare the results.
1. Interval Bisection Method
2. Newton-Raphson Method
3. Secant Method

Note that the following functions are the root finding method. In order to find their maximum point, I am going to use their derivative function.

Evaluation

First,let’s evaluate the function.

We can graph the function and see where might be the local maximum point.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
x.seq <- seq(from = 0.1, to = 10, length.out = 101)
g <- sapply(x.seq, function(x) log(x)/(1+x))


# put the x and y values in a data frame
plot_data <- data.frame("x" = x.seq, "g" = g)

head(plot_data)
##       x          g
## 1 0.100 -2.0932592
## 2 0.199 -1.3464975
## 3 0.298 -0.9327132
## 4 0.397 -0.6612878
## 5 0.496 -0.4687028
## 6 0.595 -0.3255134
# create the base layer of ggplot with geom_Line
plot_g <- ggplot(plot_data, aes(x = x, y = g)) +
          geom_line(color = "red", size = 1.5) +
          labs(x = "x", y = "g(x)")

# Increase font size in axes
plot_g + theme(axis.text.x =element_text(face = "bold",size = 12),
               axis.text.y =element_text(face = "bold", size = 12),
               axis.line = element_line(colour = "black",
                                        size = 1, linetype = "solid"))

And now plot g`(x) function.

#First Create g`(x) function
g = expression(log(x)/(x+1)) # g function
dg.x = D(g,'x') # derivative of g
dg.dx = function(x){eval(dg.x)}

#Add a new column
plot_data$dg.dx <- sapply(plot_data$x, dg.dx)

head(plot_data)
##       x          g     dg.dx
## 1 0.100 -2.0932592 10.993872
## 2 0.199 -1.3464975  5.314114
## 3 0.298 -0.9327132  3.303866
## 4 0.397 -0.6612878  2.276435
## 5 0.496 -0.4687028  1.660984
## 6 0.595 -0.3255134  1.257797
# create the base layer of ggplot with geom_Line
plot_dg.dx <- ggplot(plot_data, aes(x = x, y = dg.dx)) +
              geom_line(color = "blue", size = 0.95) +
              labs(x = "x", y = "g'(x)")
plot_dg.dx

Now let us zoom in on 3 <= x<= 4 and try to see where g`(x) is equal to 0.

# Zoom in between 3 and 4
plot_dg.dx + 
  # add a horizontal line at 0
  geom_hline(yintercept = 0.00, color = "green", size = 1) +
  coord_cartesian(xlim = c(3,4), ylim = c(-0.025, 0.025))

So, as it seen on the graph that our root is around 3.59.

Now, let’s implement our root finding methods.

Using Optimize Function

f <- function(x) log(x)/(x+1)
xmax <- optimize(f, c(1,5),tol = 0.0001, maximum = T)
xmax
## $maximum
## [1] 3.591119
## 
## $objective
## [1] 0.2784645

We can see that maximum of g(x) is 3.591119 at 0.278.



I. Interval Bisection Method

#Let's create a function
IBM <- function (a,b){
n=0; value=1; tolerance=0.0001
while (abs(value) > tolerance){
  n=n+1 
  x = (a+b)/2
  g = expression(log(x)/(x+1)) # g function
  dg.dx  = D(g,'x') # derivative of g
  value <- eval(dg.dx ) # evaluating of derivative
  if (value>0) {
    a=(a+b)/2
    b=b
    }
  else {
    a=a
    b=(a+b)/2
  }
}
list( "Maxvalue" =x, "iterations"=n)
}
IBM(1,5)
## $Maxvalue
## [1] 3.59375
## 
## $iterations
## [1] 7



II. Newton-Raphson Method

NR <- function (x){
  n <- 1
  g = expression(log(x)/(x+1)) # g function
  dg.dx  =D(g,'x') # derivative of g
  ddg.dx  =D(D(g,'x'),'x')
  f0 <- eval(dg.dx )
  f1 <- eval(ddg.dx)
  y <- x - (f0/f1)
  while (abs(f0) > 0.0001) {
    x<- y
    f0 <- eval(dg.dx )
    f1 <- eval(ddg.dx)
    y <- x - (f0/f1)
    n= n+1
  }
  list( "Maxvalue" =x, "iterations"=n)
}
NR(5)
## $Maxvalue
## [1] 3.58957
## 
## $iterations
## [1] 8



III. Secant Method

ga = expression(log(a)/(a+1))
gd_a =D(ga,'a') # derivative of g

gb = expression(log(b)/(b+1))
gd_b =D(gb,'b') # derivative of g


SEC <- function (a,b){
  n <- 1
  ga = expression(log(a)/(a+1))
  gd_a =D(ga,'a') # derivative of g
  fa <- eval(gd_a) # evaluating of derivative
  ####
  fb <- eval(gd_b) # evaluating of derivative
  c <- b - fb/((fb-fa)/(b-a))
  while (abs(c-b) > 0.0001) {
    a<- b
    b<- c
    fa <- eval(gd_a) 
    fb <- eval(gd_b) 
    c <- b - fb/((fb-fa)/(b-a))
    n= n+1
  }
  list( "Maxvalue" =c, "iterations"=n)
}

SEC(1,5)
## $Maxvalue
## [1] 3.591121
## 
## $iterations
## [1] 13
SEC(1,5)$Maxvalue
## [1] 3.591121



Method Comparisons

df_new <- data.frame(
  id = 1:4,
  Model = c('Optimize function', 'IBM', 'Newton-Raphson', 'Secant'),
  Value =  c(xmax$maximum, IBM(1,5)$Maxvalue,  NR(5)$Maxvalue,  SEC(1,5)$Maxvalue), 
  Iteration =  c(NA,IBM(1,5)$iterations,  NR(5)$iterations,  SEC(1,5)$iterations))
library(formattable)
## Warning: package 'formattable' was built under R version 4.0.5
first_formatter <- formatter("span", 
                                 style = ~ style(color = "grey",
                                                 font.weight = "bold"))

formattable(df_new,
            list(Model = first_formatter,
                 Value = first_formatter,
                 Iteration = first_formatter))
id Model Value Iteration
1 Optimize function 3.591119 NA
2 IBM 3.593750 7
3 Newton-Raphson 3.589570 8
4 Secant 3.591121 13



References

  1. Colorado State MAS Program Computational Methods Class Notes
  2. https://en.wikipedia.org/wiki/Mathematical_optimization
  3. https://en.wikipedia.org/wiki/Secant_method
  4. https://brilliant.org/wiki/newton-raphson-method/#




**************