Lab4: Solving Nonlinear Equations

1 Bisection Method

The bisection method is another approach to finding the root of a continuous function \(f(x)\) on an interval \([a,b]\). 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.

1.1 Matematical Formula

first we calculating the midpoint, \[m=a/2+b/2\], of the interval. The function is then evaluated at that point \(f(m)\). Then , recalling Bolzano’s theorem, if \(f(a)\) and \(f(m)\) have opposite signs, the method replaces b with the computed midpoint, or if \(f(b)\) and \(f(m)\) have opposite signs, a is replaced by the midpoint. This step ensures there is still a root contained in the interval. Then the procedure continues to the next iteration. The solution is said to be found once the function equals to 0 at \(f(m)\) or is small enough to be sufficient.

1.2 The Bisection Algorithm

The bisection algorithm is a simple method for finding the roots of one-dimensional functions. The goal is to find a root
x0∈[a,b] such that \(f(x0)=0\). The algorithm starts with a large interval, known to contain x0, and then successively reduces the size of the interval until it brackets the root. The theoretical underpinning of the algorithm is the intermediate value theorem which states that if a continuous function f takes values \(f(a)\) and \(f(b)\) at the end points of the interval \([a,b]\), then f must take all values between \(f(a)\) and \(f(b)\) somewhere in the interval. So if \(f(a)<γ<f(b)\), then there exists \(a c ∈[a,b]\) such that \(f(c)=γ\).

Using this information, we can present the bisection algorithm. First we must check that sign(f(a))≠sign(f(b)). Otherwise, the interval does not contain the root and might need to be widened. Then we can proceed:

1. Let \(c=a/2+b/2\).
2. If \(f(c)=0\), stop and return \(c.\)
3. If sign\((f(a))≠sign(f(c))\), then set \(b←c\). Else if \(sign(f(b))≠sign(f(c)), then set a←c.\)
4. Goto the beginning and repeat until convergence (see below).

After n iterations, the size of the interval bracketing the root will be \(2^−n(b−a)\).

The bisection algorithm is useful, conceptually simple, and is easy to implement. In particular, you do not need any special information about the function

except the ability to evaluate it at various points in the interval. The downsides are that it is only useful in one dimension and its convergence is linear, which is the slowest rate of convergence for algorithms we will discuss (more on that later).

1.3 Bisection Algorithm in R

Suppose we want to find the root of the function \(x^3−2x−5\). It is often valuable to plot the function to find an interval containing the root.

1.3.1 Visualize the function

func <- function(x) {
  x^3 - 2 * x - 5
}

curve(func, xlim=c(-3,3), col='blue', lwd=1.5, lty=2)
abline(h=0)
abline(v=0)

As we can see from the visualization above, we know that the root appears to lie close in 2. Next, we will use the interval [2,3] in Bisection Method.

1.3.2 Use NLRoot packages

On of the most usefull package for the Bisection Method in R is the NLRoot package, here you will apply BFfzero() to find the solution, as you can see below:

library(NLRoot)              # loading the package 
BFfzero(func, 2, 3)          # to find the solution (just as simple like this)

## [1] 1
## [1] 2.094553
## [1] 1.262094e-05
## [1] "finding root is successful"

The output of the function shows the root is located at \(x=2.094553\), which matches the results found in previous examples on other root finding methods such as the Newton-Raphson method and Secant method.

1.3.3 Write Function

We can also write a function that implements the bisection method using the iteration steps detailed at the beginning of the post.

Bakti_Bisection <- function(f, a, b, n = 1000, tol = 1e-7) {
  # If the signs of the function at the evaluated points, a and b, stop the function and return message.
  if (!(f(a) < 0) && (f(b) > 0)) {
    stop('signs of f(a) and f(b) differ')
  } else if ((f(a) > 0) && (f(b) < 0)) {
    stop('signs of f(a) and f(b) differ')
  }
  
  for (i in 1:n) {
    c <- (a + b) / 2 # Calculate midpoint
    
    # If the function equals 0 at the midpoint or the midpoint is below the desired tolerance, stop the 
    # function and return the root.
    if ((f(c) == 0) || ((b - a) / 2) < tol) {
      return(c)
    }
    
    # If another iteration is required, 
    # check the signs of the function at the points c and a and reassign
    # a or b accordingly as the midpoint to be used in the next iteration.
    ifelse(sign(f(c)) == sign(f(a)), 
           a <- c,
           b <- c)
  }
  # If the max number of iterations is reached and no root has been found, 
  # return message and end function.
  print('Too many iterations')
}

Now, we can apply our function find the root of the function \(x^3−2x−5\) as we have done above by using NLRoot packages.

Bakti_Bisection(func, 2, 3)

## [1] 2.094552

1.3.4 Your Own Function

funct <- function(x) {
  x^3 - 2 * x - 6
}

curve(funct, xlim=c(-3,3), col='red', lwd=1.5, lty=2)
abline(h=0)
abline(v=0)

library(NLRoot)              
BFfzero(funct, 2, 3) 

## [1] 1
## [1] 2.179984
## [1] 3.000934e-05
## [1] "finding root is successful"

x_Bisection <- function(f, a, b, n = 10000, tol = 1e-9) {
  if (!(f(a) < 0) && (f(b) > 0)) {
    stop('signs of f(a) and f(b) differ')
  } else if ((f(a) > 0) && (f(b) < 0)) {
    stop('signs of f(a) and f(b) differ')
  }
  
  for (i in 1:n) {
    c <- (a + b) / 2 
    
    if ((f(c) == 0) || ((b - a) / 2) < tol) {
      return(c)
    }
    
  ifelse(sign(f(c)) == sign(f(a)), 
           a <- c,
           b <- c)
  }
  print('Too many iterations')
}
x_Bisection(funct, 2,3)

## [1] 2.179981

2 Newton’s Method

The Newton-Raphson method is a method for approximating the roots of polynomial equations of any order. In fact the method works for any equation, polynomial or not, as long as the function is differentiable in a desired interval. ’s often become increasingly better approximations of the function’s root. Newton’s method usually use for solving equations is another numerical method for solving an equation \(f(x)=0\). It is based on the geometry of a curve, using the tangent lines to a curve. As such, it requires calculus, in particular differentiation. Roughly, the idea of Newton’s method is as follows.

2.1 Matematical Formula

Newton’s method build a sequence of values \({x_n}\) via functional iteration that converges to the root of a function \(f\). Let that root be called \(x_∞\) and let \(x_n\) be the current estimate. By the mean value theorem, we know there exists some \(z\) such that \(f(x_n)=f′(z)(x_n−x_∞)\), where z is somewhere between xn and x∞. Rearranging terms, we can write \(x_∞=x_n− f(x_n)/f′(z)\) Obviously, we do not know x_∞ or z, so we can replace them with our next iterate xn+1 and our current iterate xn, giving us the Newton update formula, \(x_n+1=x_n−f(x_n)f′(x_n).\)

2.2 The Newton’s Algorithm

Newton’s method (also known as the Newton-Raphson method or the Newton-Fourier method) is an efficient algorithm for finding approximations to the zeros (or roots) of a real-valued function f(x).

The iteration goes on in this way:

\(x_{k + 1} = x_{k} - \frac{FUN(x_{k})}{FUN'(x_{k})}\)

From the starting value \(x_0\), vertical lines and points are plotted to show the location of the sequence of iteration values \(x_1, x_2, \ldots\); tangent lines are drawn to illustrate the relationship between successive iterations; the iteration values are in the right margin of the plot.

2.3 Newton’s Method in R

Suppose we want to find the root of the function \(f(x)=x^2−9 \text{ for, } x>0\). We assume that \(f\) is smooth (differentiable). Our goal is to find a root of \(f\) i.e a \(x\) value for which \(f(x)=0\). We illustrate the method with a simple example for which you know the root.

2.3.1 Visualize the function

# first generate the sequence
x.seq <- seq(from = 0.1, to = 10, length.out = 101)
f <- sapply(x.seq, function(x) x^2 - 9)

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.6.3

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

##       x         f
## 1 0.100 -8.990000
## 2 0.199 -8.960399
## 3 0.298 -8.911196
## 4 0.397 -8.842391
## 5 0.496 -8.753984
## 6 0.595 -8.645975

# create the base layer of ggplot with geom_Line
plot_f <- ggplot(plot_data, aes(x = x, y = f)) + 
  geom_line(color = "red", size = 1.0) + 
  labs(x = "x", y = "f(x) = x^2 - 9")
# beautify : Increase font size in axes 
plot_f + 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"))

By considering the visualization above, it’s looks like the function is 0 around 3. How do we know that exactly correct, let’s try zoom the visualization.

plot_f + 
  geom_hline(yintercept = 0.00, color = "green", size = 1) +
  coord_cartesian(xlim = c(2,4), ylim = c(-1, 3))

Next, let us visualize this to find the root of our function \(f(x)=x^2−9.\) First we need to create a function that returns the tangent line to the function f at some point a.

# Let us create our tangent line function at some point a
f_tangent <- function(x, f, df, a){
  # We need the function and its derivative df
  # a is the point where we create the tangent to the function
  # returns the tangent line 
  df(a)*(x - a) + f(a)
}

Then, plot the tangent line to the initial value \(f\) at \(x(0)=4\) which will be our first guess of the root.

# add the tangent line at 4 to our data frame
plot_data$f.tan <- sapply(plot_data$x, f_tangent, f = function(x) x^2 - 9, df = function(x) 2*x, a = 4)

# create the base layer of ggplot with geom_Line
plot_tangent_line <- ggplot(plot_data, aes(x = x)) + 
  geom_line(aes(y = f), color = "red", size = 0.95) +
  geom_line(aes(y = f.tan), color = "blue", size = 0.95) +
  labs(x = "x", y = "f(x) and tangent line at 4") 
plot_tangent_line

Next, let us zoom in on \(2≤x≤4\) and try to see where the tangent line is equal to 0.

# zoom in between 3 and 4
plot_tangent_line + 
  # add a horizontal line at 0
  geom_hline(yintercept = 0.00, color = "green", size = 1) +
  coord_cartesian(xlim = c(2,4), ylim = c(-5, 10))

By considering the visualization above, looks like the tangent line is 0 around 3.1. This will be our updated guess at the root. Understand the rationale behind this. If f is approximately equal to its tangent line at 4, then the root of \(f\) is approximately equal to the root of the tangent line. Can you see this in the plot above?

The exact value when the tangent line is 0 is given by

\[ x^{(1)}=4−{f(4) \over f′(4)} =3.125\]

Then, add the tangent line at 3.125 to our data frame

plot_data$f.tan <- sapply(plot_data$x, f_tangent, f = function(x) x^2 - 9, df = function(x) 2*x, a = 3.125)

# create the base layer of ggplot with geom_Line
plot_tangent_line <- ggplot(plot_data, aes(x = x)) + 
  geom_line(aes(y = f), color = "red", size = 0.95) +
  geom_line(aes(y = f.tan), color = "blue", size = 0.95) +
  labs(x = "x", y = "f(x) and tangent line at 4") 
plot_tangent_line

Now let us zoom in on around 3 and try to see where the tangent line is equal to 0.

# zoom in between 2.9 and 3.1
plot_tangent_line + 
  # add a horizontal line at 0
  geom_hline(yintercept = 0.00, color = "green", size = 1) +
  coord_cartesian(xlim = c(2.9,3.1), ylim = c(-2.5, 2.5))

Looks like the tangent line is 0 around 3. In just two steps we are very close. This will be our updated guess at the root. The exact value when the tangent line is 0 is given by

\[ x^{(2)}=x^{(1)}−{f(x^{(1)}) \over f′(x^{(1)})}\] \[ x^{(2)}=3.125−{f(3.125) \over f′(3.125)} = 3.0025\]

We can continue until our guesses keep getting smaller and smaller.

2.3.2 Use NLRoot packages

The same packages as we use to the Bisection Method, this packages also available for Newton’s method.

library(NLRoot)
# first argument is the function whose root we need to find
# second argument is the the derivative of the function whose 
# root we need to find.
# third argument is an initial guess
NIMfzero(function(x) x^2-9, function(x) 2*x, x0 = 4)

## [1] 3
## [1] 1.081801e-12
## [1] "finding root is successful"

3 Case Study

In this lesson we will construct two functions that you will need for this course and begin to explore the R programming language. Run through the Chapter 1 of the book which talks about model building for a one-variable optimization problem. Example 1.1 (Meerschaert) A pig weighing 200 lbs gains 5 lbs per day and costs 45 cents per day to keep (feed and board). The market price for pigs is 65 cents per pound, but is falling 1 cent per day. When should the pig be sold?

Modeling the Problem:

Variables: t = time (days) w = weight of the pig(lbs) p = price for pigs (\(/lb) C = cost of keeping pig t days (\)) R = revenue obtained by selling pig (\() P = profit from selling pig (\))

Assumptions: \(w=200+5t\) \(p=0.65−0.01t\) \(C=0.45t\) \(R=p⋅w\) \(P=R−C\) \(t \geq 0\)

Objective: Maximize profits from the sell of the pig Let \(x\) be the time to sell the pig and develop the profit function using x as the independent variable. In the case of an optimization, the variable on which you will make a decision on is generally called a decision Profit Function = Objective function Given the description and mathematical explanation of the variables, lets write out the profit function: profit \((x)\) = (0.65-0.01x)(200+5x)-0.45x Let’s use R to “see” what the profit function looks like over the next 20 days.

profit = function (x){
  return((0.65-0.01*x)*(200+5*x)-.45*x)
}
x = seq(0,20,1)
plot(x,profit(x),type="o")

print(profit(8))

## [1] 133.2

As you can see, there is a good and bad day to sell the pig if our assumptions remain the same over the next 20 days. Although we could use this graph in order to select the “winner”, we want to be able to use this tool (R programming language) to assist us in ensuring that we have all of the information needed to make a decision. That means that we need to get our hands dirty and figure out how to use our tools (in this case, R).

##Build the functions that we need to solve the problem We will create two functions that will be able to solve for an optimal point in R using mathematics / programming that we already know how to do.

Using this simple formula, let’s build a derivative of the profit function that we defined above. Notice how our function fprime(f,x) will work for ANY function that we assign to it.

x=seq(0,20)
fprime = function (f,a,h=0.0001){(f(a+h)-f(a-h))/(2*h)}
plot(x,fprime(profit,x))

##Root-Finding Solver Although we could use the graph above to find the solution for this particular problem, it is more appropriate to use math to make our solution automatic. There are many root-finding algorithms out there and feel free to use another, but we will use a simple algorithm known as the bisection method. This method only requires that you select two points that have opposite signs: if \(f(a)=− then f(b)=+\).

##Bisection Algorithm for finding the root of an equation. The simple idea of the bisection method is that you start with an interval [a,b] that has a positive and negative value for the function \(f(a)\) and \(f(b)\) - it does not matter which one is positive. You cut the interval in half and proceed with the interval that keeps one side positive and the other negative. In Figure 1, we start on the interval [2,9] which works because \(f(2)\) is positive and \(f(9)\) is negative. The first step cuts the interval in half and the decision must be made to continue with [2,5.5] or [5.5,9]. It should be pretty obvous from the picture that [2,5.5] is the right answer (we want to keep a root in the interval), but from an algorithm perspective, we simply test f(5.5) and find out that it is negative which causes us to replace 9 (\(f[9]\) was negative) with 5.5. We would continue this algorithm until we are satisfied that the width of the interval is small enough to make a conclusion about the value of the root.

bisection = function(f,a,b,tol=0.0001){
  if (f(a)*f(b) > 0){
    return ("Boundary Conditions Not Met")
  }
  else{
    middle = a
    while (abs(f(middle))>tol){
      middle = (a+b)/2
      if (f(middle)*f(a)>0) (a = middle)
      else (b = middle)
      x=middle
      y=f(middle)
      ## if you want to "see" what happens at every step, take off the # of the next line ##
      #cat(sprintf("x-Val: %.4f ; f(x-val): %.4f\n",x,y))
    }
    return (middle)
  }
}
#Example of a finding a root of a function
f = function (x){x^3-5} #sinlge line function definition
zero=bisection (f,-10,20)
x=seq(-1,2,0.01)
plot(x,f(x),"l")
abline(h=0,col="red")
abline(v=zero,col="purple")

print(zero)

## [1] 1.709967

Use our tool to solve the pig problem Use the tools at your disposal to solve the problem. Be able to make several alternatives to the proposed model and investigate the relationship between multiple variables. When do we sell the pig?

dProfit = function(x){fprime(profit,x)} #derivative of profit function
bisection(dProfit,0,20)

## [1] 8.000488

As you can see, we got approximately the same answer as the book. The difference is that our solution uses an estimation of the derivative instead of the closed-form solution as well as a solver for finding the day x where the profit is maximized (profit function derivative = 0). In the case of most problems, you will use numerical methods (like we did) because real world problems do not ususally have a easy answer. To become a good problem solver, you are going to need to invest in your coding abilities because you will usually need to develop some original code to help you solve a particular problem.

Sensitivity Analysis: Price of the Pig Let’s go back and examine more closely one of our original assumptions that may or may not hold to see how our solution is afftected by a different assumption. Does the optimal day to sell change if we assume some different price of the pig. Here we let the pig price drop to recreate the numbers in the book, but we are able to do so much more which will be demonstrated in the next analysis.

#Price falling per day = p 
p = seq(0.008,0.012,0.001)
ans = array(0,length(p))
for (i in 1:length(p)){
    profit = function (x){
      return((0.65-p[i]*x)*(200+5*x)-.45*x)
    }
    dProfit = function(x){fprime(profit,x,)}
    ans[i] = bisection(dProfit,0,20,0.0001)
}
print(ans)

## [1] 15.000000 11.110840  8.000488  5.454102  3.332520

plot(p,ans,"o",xlab="p($/day)",ylab="x(Days to Sell)")
title("Sensitivity of Falling Price of Pig")

Sensitivity Analysis: Growth Rate of the Pig Since we are also concerned about the growth rate of the pig, let’s conduct see how the days to sell are affected by the growth rate. The weight (w) is calculated using the assumption that the pig will grow 5lbs per day so that: \[w=200+gt\] where \[g=5\]. Let’s see what happens if g is between 3 and 7. Remember the bisection method needs two boundary points that have the opposite values and if the growth rate is too small, think about what that does to the profit. Experiment with large and small growth rates and the profit function to figure it out for yourself if you don’t see it. We needed to make the starting point of the bisection method to the left (negative number) in order for it to find the zero for 3lbs and 3.5lbs per day.

#Growth rate of the pig = g
g = seq(1,12,0.5)
ans = array(0,length(g))
for (i in 1:length(g)){
    profit = function (x){
      return((0.65-0.01*x)*(200+g[i]*x)-.45*x)
    }
    dProfit = function(x){fprime(profit,x,)}
    ans[i] = bisection(dProfit,-100,50,0.0001)
}
print(ans)

##  [1] -90.002441 -49.169922 -28.749084 -16.499329  -8.332825  -2.500916
##  [7]   1.875305   5.278778   8.000183  10.227203  12.083435  13.653564
## [13]  14.999390  16.166687  17.187500  18.088150  18.889236  19.605637
## [19]  20.249939  20.833588  21.363449  21.847534  22.291565

plot(g,ans,"o",xlab="g(growth/day)",ylab="x(Days to Sell)")
title("Sensitivity of Growth Rate of Pig")