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

The Interval is divided into half.

\[c=\frac{(a+b)}{2}\]

• The function is evaluated at \(c\)

• Location of the root is determined as lying within the subinterval for the next iteration (\(c\) is replaced either by \(a\) or \(b\)).

• This process is repeated until the desired precision.

1.2 The Bisection Algorithm

For a given function f(x),the Bisection Method algorithm works as follows:

1. Two values a and b are chosen for which \(f(a)>0\) and \(f(b)<0\) (or the other way around)
2. Interval halving: A midpoint \(c\) is calculated as the arithmetic mean between \(a\) and \(b\),\(c=(a+b)/2\)
3. The function \(f\) is evaluated for the value of \(c\)
4. If \(f(c)=0\) means that we found the root of the function, which is \(c\)
5. If \(f(c)\neq 0\) we check the sign of \(f(c)\):
   • if \(f(c)\) has the same sign as \(f(a)\), we replace \(a\) with \(c\) and we keep the same value for \(b\).
   • if \(f(c)\) has the same sign as \(f(b)\), we replace \(b\) with \(c\) and we keep the same value for \(a\)
6. We go back to step 2. and recalculate \(c\) with the new value either it is value of \(a\) or \(b\)

The algorithm ends when the values of \(f(c)\) is less than a defined tolerance (e.g. 0.001). In this case we say that \(c\) is close enough to be the root of the function for which \(f(c) ~= 0\).

In order to avoid too many iterations, we can set a maximum number of iterations (e.g. 1000) and even if we are above the defined tolerance, we keep the last value of c as the root of our function.

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

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

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

library(NLRoot)              # loading the package 
BFfzero(func, 2, 3)

## [1] 1
## [1] 2.462045
## [1] 8.100567e-06
## [1] "finding root is successful"

Imelda_Bisection <- function(f, a, b, n = 1000, 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')
}

Imelda_Bisection(func, 2, 3)

## [1] 2.462045

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

If \(x_n\) is an approximation a solution of \(f(x)=0\) and if \(f'(x_n)\neq 0\) the next approximation is given by,

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}\]

2.2 The Newton’s Algorithm

1. Find points \(a\) and \(b\) such that \(a < b\) and \(f(a).f(b) < 0\).
2. Take the interval \([a,b]\) and find next value \(x_0=\frac{(a+b)}{2}\)
3. Find \(f(x_0)\) and \(f'(x_0)\), \(x_1=x_0-\frac{f(x0)}{f'(x_0)}\)
4. If \(f(x_1)=0\) then \(x_1\) is an extract root, else \(x_0=x_1\).
5. Repeat steps 2 to 4 until \(f(x_i)=0\) or \(|f(x_i)|\leq Accuracy\).

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)
# 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

1. 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 ($)

2. 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:

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).

3.1 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))

3.2 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)=+\).

3.3 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

3.4 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.

3.5 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")

3.6 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.5 lbs 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")