Root Finding: The Bisection Method

Anthony Jimenez
17 August 2021
http://www.fa-jimenez.com/

Import necessary libraries

import numpy as np
import pandas as pd
import plotly.express as px
from plotly.subplots import make_subplots

Define complex function we which to find the root of

$f(x) = x^3 + 4x^2 - 10$

Pseudocode for bisection algorithm

INPUT: x_left, x_right, tolerance, max_iterations

set i == 0:

• Compute y_left

while i < max_iterations:

• Compute x_mid

• Compute y_mid

  if y_left * y_mid > 0:

  • x_left = x_mid
  • y_left = y_mid

  else:

  • x_right = x_mid
  • y_right = y_mid
    
    
• Check convergence
• Index i

1) Defining the left and right starting limits

x_left = 1
x_right = 2
tolerance = 1e-4
max_iteration = 20

2) Define the smooth, continuous function

def root_function(x):
    return x**3 + 4 * x**2 - 10

3) Code the main bisection algorithm

def bisection_algorithm(x_left, x_right, tolerance, max_iteration):
    # Initial computations and creating lists to be appended to
    x_left = [x_left]
    x_right = [x_right]
    y_left = [root_function(x_left[0])]
    y_right = [root_function(x_right[0])]
    x_mid = []
    y_mid = []
       
    # Create index counter for while loop
    i = 0
    while i < max_iteration:        
        # Compute the midpoint of the interval and the y-value
        x_mid.append(x_left[i] + (x_right[i] - x_left[i]) / 2)
        y_mid.append(root_function(x_mid[i]))
        
        # If we converge, assign the appropriate values and exit
        if (y_mid[i] == 0) or ((x_right[i] - x_left[i]) / 2 < tolerance):
            break
        
        # If we do not converge, index the counter and update list values
        i += 1
        # If greater than 0, that means x_left and x_mid are both negative or both positive (same side of y-axis)
        if y_left[i-1] * y_mid[i-1] > 0:
            # Update the x_left to be the previously computed mid point
            x_left.append(x_mid[i-1])
            y_left.append(y_mid[i-1])
            
            # The right side x- and y-values remain the same
            x_right.append(x_right[i-1])
            y_right.append(y_right[i-1])
        
        # If x_left and x_mid are on opposite sides, we update the x_right values instead
        else:
            # Update x_right
            x_right.append(x_mid[i-1])
            y_right.append(y_mid[i-1])
            
            # Maintain x_left
            x_left.append(x_left[i-1])
            y_left.append(y_left[i-1])
    
    # Create pandas dataframe that will serve as a summary report table
    iter_index = np.arange(1, len(x_left)+1)
    cols = ['Iteration', 'x_left', 'x_right', 'x_midpoint', 'y_midpoint', 'y_left', 'y_right']
    df = pd.DataFrame(data=[iter_index, x_left, x_right, x_mid, y_mid, y_left, y_right]).T
    df.columns = cols
    return df

4) Call the main function

df = bisection_algorithm(x_left, x_right, tolerance, max_iteration)
df

	Iteration	x_left	x_right	x_midpoint	y_midpoint	y_left	y_right
0	1.0	1.000000	2.000000	1.500000	2.375000	-5.000000	14.000000
1	2.0	1.000000	1.500000	1.250000	-1.796875	-5.000000	2.375000
2	3.0	1.250000	1.500000	1.375000	0.162109	-1.796875	2.375000
3	4.0	1.250000	1.375000	1.312500	-0.848389	-1.796875	0.162109
4	5.0	1.312500	1.375000	1.343750	-0.350983	-0.848389	0.162109
5	6.0	1.343750	1.375000	1.359375	-0.096409	-0.350983	0.162109
6	7.0	1.359375	1.375000	1.367188	0.032356	-0.096409	0.162109
7	8.0	1.359375	1.367188	1.363281	-0.032150	-0.096409	0.032356
8	9.0	1.363281	1.367188	1.365234	0.000072	-0.032150	0.032356
9	10.0	1.363281	1.365234	1.364258	-0.016047	-0.032150	0.000072
10	11.0	1.364258	1.365234	1.364746	-0.007989	-0.016047	0.000072
11	12.0	1.364746	1.365234	1.364990	-0.003959	-0.007989	0.000072
12	13.0	1.364990	1.365234	1.365112	-0.001944	-0.003959	0.000072
13	14.0	1.365112	1.365234	1.365173	-0.000936	-0.001944	0.000072

5) Visualizations

First show the overall function behavior

x_plt = np.linspace(-10, 10, 250)
y_plt = root_function(x_plt)
fig = px.scatter(x=x_plt, y=y_plt).update_traces(mode='lines')
fig.update_layout(title='General overview of the defined function (note the cubic nature)')
fig.show()

Create animation zoomed in on interval [1,2] to see the root finding

fig = px.scatter(df, x='x_midpoint', y='y_midpoint', animation_frame='Iteration',
                range_x=[1, 1.6], range_y=[-2, 3])
fig.add_scatter(x=x_plt, y=y_plt, name=r'$f(x) = x^3 + 4x^2 - 10$')
fig.update_traces(marker={'size': 12})
fig.update_layout(title='Animation figure for root finding (bisection algorithm)', xaxis_title='x', yaxis_title='y')
fig.update_layout(legend={
                     'orientation': 'h',
                     'yanchor': 'bottom',
                     'xanchor': 'right',
                     'x': 1,
                     'y': 1.02
                 })
fig.show()

6) Potential Downfalls

Long convergence time and unbiased approach to newly-guessed root value

The bisection algorithm is very dependent on a good first guess. Additionally, a function that may have multiple zeros can provide difficulties to the algorithm.

Further, because the algorithm always divides the space between x_left and x_right by 2, there can be a long time for convergence when oscillating around the zero.

x_left = -5
x_right = 3
tolerance = 1e-4
max_iteration = 20

df2= bisection_algorithm(x_left, x_right, tolerance, max_iteration)
df2

	Iteration	x_left	x_right	x_midpoint	y_midpoint	y_left	y_right
0	1.0	-5.000000	3.000000	-1.000000	-7.000000	-35.000000	53.000000
1	2.0	-1.000000	3.000000	1.000000	-5.000000	-7.000000	53.000000
2	3.0	1.000000	3.000000	2.000000	14.000000	-5.000000	53.000000
3	4.0	1.000000	2.000000	1.500000	2.375000	-5.000000	14.000000
4	5.0	1.000000	1.500000	1.250000	-1.796875	-5.000000	2.375000
5	6.0	1.250000	1.500000	1.375000	0.162109	-1.796875	2.375000
6	7.0	1.250000	1.375000	1.312500	-0.848389	-1.796875	0.162109
7	8.0	1.312500	1.375000	1.343750	-0.350983	-0.848389	0.162109
8	9.0	1.343750	1.375000	1.359375	-0.096409	-0.350983	0.162109
9	10.0	1.359375	1.375000	1.367188	0.032356	-0.096409	0.162109
10	11.0	1.359375	1.367188	1.363281	-0.032150	-0.096409	0.032356
11	12.0	1.363281	1.367188	1.365234	0.000072	-0.032150	0.032356
12	13.0	1.363281	1.365234	1.364258	-0.016047	-0.032150	0.000072
13	14.0	1.364258	1.365234	1.364746	-0.007989	-0.016047	0.000072
14	15.0	1.364746	1.365234	1.364990	-0.003959	-0.007989	0.000072
15	16.0	1.364990	1.365234	1.365112	-0.001944	-0.003959	0.000072
16	17.0	1.365112	1.365234	1.365173	-0.000936	-0.001944	0.000072

fig = px.scatter(df2, x='x_midpoint', y='y_midpoint', animation_frame='Iteration'
                )
fig.add_scatter(x=x_plt, y=y_plt, name=r'$f(x) = x^3 + 4x^2 - 10$')
fig.update_traces(marker={'size': 12})
fig.update_layout(title='Example of how convergence can take a long time', xaxis_title='x', yaxis_title='y')
fig.update_layout(legend={
                     'orientation': 'h',
                     'yanchor': 'bottom',
                     'xanchor': 'right',
                     'x': 1,
                     'y': 1.02
                 })
fig.show()