E19 Root-Finding

Method of Bisection

'Manually', i.e., without a loop

Let's start with a function whose root we would like to find.

import numpy as np
import matplotlib.pyplot as plt
# Define a function whose roots we will find.
def f1(x):
  return np.exp(3-x)*np.cos(2*x) + 4*np.sin(3*x) # I just made this up
# Generate pairs of values so we can plot them.
x_nums = np.linspace(0,6)
y_nums = f1(x_nums)
plt.plot(x_nums,y_nums)
plt.grid()
0.9s

Let's hunt for the first root, which we can see is somewhere near x = 1. Let's start with the interval (0.5, 1.2). Evaluate f(x) at both these values:

x_left = 0.5
x_right = 1.2
print(f"at x = 0.5, f(x) = {f1(x_left):.2f}") # This is the syntax for printing 
print(f"at x = 1.2, f(x) = {f1(x_right):.2f}") # things with desired number of sig. figs.
plt.plot(x_nums,y_nums)
plt.xlim([x_left,x_right])
plt.grid()
0.6s

So, we see that the signs are different, and therefore, a root probably lies between these. Let's find the midpoint and check the sign of f1 there.

x_mid = (1/2)*(x_left+x_right)
print(f"At x = {x_mid}, f(x) = {f1(x_mid):.2f}")
0.3s

This is positive. So the root must be somewhere between x_mid and x_right.

Re-designate x_mid as the new x_left and try again.

x_left = x_mid
print(f"at x = {x_left}, f(x) = {f1(x_left):.2f}") #  
print(f"at x = {x_right}, f(x) = {f1(x_right):.2f}") # 
plt.plot(x_nums,y_nums)
plt.xlim([x_left,x_right])
plt.grid()
0.8s

Now, find the midpoint of the interval and check the sign of f1 there.

x_mid = (1/2)*(x_left+x_right)
print(f"At x = {x_mid}, f(x) = {f1(x_mid):.2f}")
0.3s

So now, we find that at the midpoint, f1 is negative. This means that the interval must lie between x_left and x_mid.

Re-designate x_mid as the new x_right this time.

x_right = x_mid
print(f"at x = {x_left}, f(x) = {f1(x_left):.2f}") #  
print(f"at x = {x_right}, f(x) = {f1(x_right):.2f}") # 
plt.plot(x_nums,y_nums)
plt.xlim([x_left,x_right])
plt.grid()
0.5s

The root lies between 0.85 and 1.025. Let's find the midpoint and check the sign of f1 there.

x_mid = (1/2)*(x_left+x_right)
print(f"At x = {x_mid}, f(x) = {f1(x_mid):.2f}")
0.3s

It's negative! So once again, the root must lie between x_left and x_mid. We can forget about the current x_right, and re-designate x_mid as x_right.

x_right = x_mid
print(f"at x = {x_left}, f(x) = {f1(x_left):.2f}") #  
print(f"at x = {x_right}, f(x) = {f1(x_right):.2f}") # 
plt.plot(x_nums,y_nums)
plt.xlim([x_left,x_right])
plt.grid()
0.5s

So we can tell that the root is somewhere between 0.85 and 0.9375.

0.0s

With a loop

from numpy import sign
x_l = 0.5
x_r = 1.2
ep = 0.05
while (abs(x_l - x_r) > ep):
  fl = f1(x_l) # evaluate f1 at x_l
  fr = f1(x_r) # evaluate f1 at x_r
  print(f"The root lies between {x_l} and {x_r}. The size of the interval is {abs(x_l - x_r):.3f}")
  
  x_m = 0.5*(x_l + x_r) # calculate the mid-point
  fm = f1(x_m) # evaluate f1 at x_m
  
  if sign(fm) == sign(fl):
    # i.e., f1 at the midpoint has the same sign as f1 at the left edge of the interval
    # and so the 'root', i.e., the value of x at which f1 = 0, must lie between
    # x_m and x_r.
    x_l = x_m
  elif sign(fm) == sign(fr):
    # i.e., f1 at the midpoint has the same sign as f1 at the right edge of the interval
    # and so the 'root', i.e., the value of x at which f1 = 0, must lie between
    # x_l and x_m.
    x_r = x_m
# outside the while loop
print("Loop has exited.")
0.3s

Secant method

Let's start with two estimates, this time on the same side: 0.5 and 0.6

import numpy as np
import matplotlib.pyplot as plt
# Define a function whose roots we will find.
def f1(x):
  return np.exp(3-x)*np.cos(2*x) + 4*np.sin(3*x) # I just made this up
# Generate pairs of values so we can plot them.
x_nums = np.linspace(0,6)
y_nums = f1(x_nums)
# Plot the function
plt.plot(x_nums,y_nums)
plt.grid()
# Plot overlay points showing the first two estimates
x1,x2 = 0.5,0.6
y1,y2 = f1(x1),f1(x2)
plt.scatter([x1,x2], [y1,y2], color='red')
0.2s

Next, we will make a straight-line estimate of the function, written here as 'y in terms of x'.

formula not implemented
def f1_straight_line_estimate(x):
  return ((y2-y1)/(x2-x1))*(x-x1) + y1
plt.plot(x_nums,y_nums)
plt.xlim(left=0,right=2)
plt.grid()
plt.scatter([x1,x2], [y1,y2], color='red')
# Now, create a set of points on which to plot the straight-line approximation y(x)
x_set1 = np.linspace(0,2)
y_set1 = f1_straight_line_estimate(x_set1)
plt.plot(x_set1,y_set1,linestyle='--',linewidth=1)
0.5s

Looks like we have an estimate that works really well -- let's solve for the point where this straight line intersects the x-axis.

formula not implemented

x3 = -y1*(x2-x1)/(y2-y1) + x1
print(f'The straight line intersects the x-axis at x = {x3:.4f}')
0.3s

This is already a pretty good estimate. But let's do one more iteration, this time using x2 and x3.

# x3 was already calculated
y3 = f1(x3)
# Plot the original function
plt.plot(x_nums,y_nums)
plt.grid()
# Plot the three points
plt.scatter([x1,x2,x3], [y1,y2,y3], color='red')
0.6s

Next, we will make a straight-line estimate of the function, this time using linear interpolation of x2 and x3.

formula not implemented
def f1_straight_line_estimate_2(x):
  return ((y3-y2)/(x3-x2))*(x-x2) + y2
plt.plot(x_nums,y_nums)
plt.xlim(left=0,right=2)
plt.grid()
plt.scatter([x1,x2,x3], [y1,y2,y3], color='red')
# Now, create a set of points on which to plot the straight-line approximation y(x)
x_set2 = np.linspace(0,2)
y_set2 = f1_straight_line_estimate_2(x_set2)
plt.plot(x_set2,y_set2,linestyle='--',linewidth=1)
0.6s

Now, we want to find out where this straight line intersects the x-axis.

formula not implemented
x4 = -y2*(x3-x2)/(y3-y2) + x2
print(f'The straight line intersects the x-axis at x = {x4:.4f}')
0.4s

Newton-Rapshon Method

In the Newton-Raphson method, we start with a single estimate for the root, and use the derivative of the function to estimate the behavior of the function.

Let's use this method to determine the root of the function inline_formula not implemented

The first step, when using Newton's method, is to determine the derivative of the function whose root you want to find. After some calculus, we find that the derivative of this function is, using the chain rule and product rule,

formula not implemented
import matplotlib.pyplot as plt
from numpy import linspace,cos
f1 = lambda x: (1-x)*(cos(x))**2
xvals = linspace(0,6,200)
yvals = f1(xvals)
plt.plot(xvals,yvals)
plt.grid()
0.6s

Cool! It looks like this has a few different roots. Let's hunt for them, starting with an estimate at 3.8

from numpy import sin,cos
f1 = lambda x: (1-x)*(cos(x))**2
df1= lambda x: (-2*cos(x)*sin(x))*(1-x) + (-1)*(cos(x))**2
# start with an estimate at x = 3.8
x_current = 3.8
for i in range(5):
  # calculate x_new
  print('The new estimate is',x_new)
  # update x_current
0.3s

The library 'scipy.optimize' has a built-in function that can find roots without any help from us. Let's see what it reports is the root for this function near x = 3.8

from scipy.optimize import fsolve
fsolve(lambda x: (1-x)*(cos(x))**2,3.8)
0.0s
Runtimes (1)