Learning Bisection Method

When I learned the Miettinen and Nurminen Test algorithm from the rate-compare article (Unstratified and Stratified Miettinen and Nurminen Test), I found that its CI is given by the roots of an equation. In order to reproduce its algorithm, I'd like to learn more about the Bisection method first.

The bisection method is an approach to finding the root of a continuous function on an interval. Its principle can be shown below (referring from https://rpubs.com/aaronsc32/bisection-method-r).

The method takes advantage of a corollary of the intermediate value theorem called Bolzano's theorem which states that if the values of f(a) and f(b) have opposite signs, the interval must contain at least one root. The iteration steps of the bisection method are relatively straightforward, however; convergence towards a solution is slow compared to other root-finding methods.

Thus the algorithm can be divided into the following steps. - Calculate the minpoint m=(a+b)/2 based on the interval with the lower bound a and higher bound b. - Get the value of function f at midpoint. If f(a) and f(m) have opposite signs, then b point is replaced by the computed midpoint, or if f(b) and f(m) have opposite signs, a is replaced by midpoint. This step ensures the root of function can be kept within the new interval. - Iterate over the above steps 1-2 until the difference between a and b is small enough or the defined iterating number is reached, at which midpoint m is the root.

According to the above steps, we can write a function to implement the bisection method in R. Assume we have the following function f to be solved.

f <- function(x) cos(x)^4 - 4*cos(x)^3 + 8*cos(x)^2 - 5*cos(x) + 1/2

Then we create a function with bisection algorithm and find the root.

bisect_func <- function(f, a, b, tol = 1e-5, n = 100) {
  i <- 1
  while (i <= n) {
    m <- (a + b) / 2
    if (abs(b - a) < tol | f(m) == 0) {
      break
    }
    
    if (f(a) * f(m) <0) {
      b <- m
    } else if (f(b) * f(m) < 0) {
      a <- m
    } else {
      stop("f(a) and f(b) must have opposite sign.")
    }
    
    i <- i + 1
    if (i > n) {
      warning("The max iterations reached, maybe the lower and upper bound should be changed and try again.")
    }
  }
  root <- (a + b) / 2
  return(root)
}

bisect_func(f, 0, 2, tol = 0.001)
[1] 0.6293945

And the same outcome can be obtained from the mature function cmna::bisection().

cmna::bisection(f, 0, 2)
[1] 0.6293945

Actually this f function should have two roots, but the bisect_func() or cmna::bisection() functions trend to have the root on the lower bound side. In this case, we must either run the bisection function multiple times for different intervals, or use another function able to finding multiple roots.

# two roots in (0,1) and (1,2)
bisect_func(f, 0, 1, tol = 0.001)
[1] 0.6293945
bisect_func(f, 1, 2, tol = 0.001)
[1] 1.447754

# using uniroot.all() function
rootSolve::uniroot.all(f, c(0, 2), n = 100)
[1] 0.6293342 1.4480258

---

Above is my brief summary and hope it can be useful. And in the following post I will learn how to use bisection method or uniroot function to solve in the Miettinen-Nurminen (Score) Confidence Limits.

Reference

Bisection Method
Finding multiple roots of univariate functions in R