In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Add some physical examples to motivate root finding. Just using algebraic examples up to this point.
Solving a linear equation (e.g. $y=mx+b$) is one of the first topic one encounters in any sort of introductory algebra course. This of course can be done by hand quite easily. One of course then moves onto more sophisticated equations involving terms like $x^2,x^3,\ldots$. For example, the equation $$x^2+10x+9=0$$ can be solved using the familiar quadratic formula $x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$ where $a=1$, $b=10$, and $c=9$. Applying the quadratic formula, we see that $x=-1,-9$. When we plot the function $f(x)=x^2+10x+9$, we can see that $f=0$ at $x=-1$ and $x=-9$.
In [20]:
x=np.linspace(-15,5,1000)
x_zeros=[-1,-9]
fig = plt.figure(figsize=(10,10))
ax = fig.gca()
ax.plot(x,0*x,'--k',linewidth=2.0)
ax.plot(x,x**2+10*x+9,linewidth=2.0)
ax.plot(x_zeros,[0,0],'ro',markersize=10)
ax.set_xlabel(r'$x$',fontsize=22)
ax.set_ylabel(r'$f(x)$',fontsize=22)
ax.set_title(r'$f(x)=x^2+10x+9$',fontsize=22)
ax.set_ylim([-20,20])
plt.show()
These types of equations, often referred to as polynomials, can also be solved through factorization. For example, the equation $$x^3+6x^2+11x+6=0$$ can be rewritten as $$(x+1)(x+2)(x+3)=0$$ such that $x=-1,-2,-3$ are all roots of the equation. Again, let's plot the function and the roots to check our answer.
In [21]:
x=np.linspace(-4,0,1000)
x_zeros=[-1,-2,-3]
fig = plt.figure(figsize=(10,10))
ax = fig.gca()
ax.plot(x,0*x,'--k',linewidth=2.0)
ax.plot(x,x**3+6*x**2+11*x+6,linewidth=2.0)
ax.plot(x_zeros,[0,0,0],'ro',markersize=10)
ax.set_xlabel(r'$x$',fontsize=22)
ax.set_ylabel(r'$f(x)$',fontsize=22)
ax.set_title(r'$f(x)=x^3+6x^2+11x+6$',fontsize=22)
ax.set_ylim([-2,2])
plt.show()
Note also that finding the value of $x$ where $f(x)=0$ is the just as easy as finding the value of $x$ where $f(x)=a$, where $a$,a constant, is some value of $f$. This is because $f(x)=a$ can be rewritten as $f(x)-a=0$. Thus, we can rewrite $f(x)=a$ as a root-finding problem $g(x)=0$ where $g(x)=f(x)-a$.
So far, we've only considered the case of finding roots of polynomials. In these cases, we could always find the roots using either the quadratic formula (when our leading term was $x^2$) or a factorization method when possible. But what if our function is more complicated, say $f(x)=2x^2\sin(x)+\ln(x)$, or better yet, what if our function is unknown altogether? Often, we resort to root finding algorithms where we make a first guess about the root of our function and then repeatedly refine our guess until we have a reasonably accurate estimate of the root.
One of the most popular root-finding methods is known as the bisection method. In order to apply the bisection method, there are two important requirements:
Let's consider what both of these conditions actually mean. Consider a function $f$ that is monotonically increasing (or decreasing) from $x=a$ to $x=b$. This is just a fancy way of saying that as we move from $a$ to $b$, $f$ only increases (or only decreases). For example, the above function ($f(x)=x^3+6x^2+11x+6$) is neither monotonically increasing or decreasing on the interval shown, $x\in[-4,0]$. However, from $x=-4$ to $x=-2-\sqrt{3}/3$, the function increases, from $x=-2-\sqrt{3}/3$ to $x=-2+\sqrt{3}/3$ the function decreases, and then from $x=-2+\sqrt{3}/3$ to $x=0$, the function increases. Thus, when applying our bisection method, we must use it on each of these intervals individually. Lets consider the middle interval, $x\in[-2-\sqrt{3}/3,-2+\sqrt{3}/3]$.
In [22]:
xzoom = np.linspace(-2.0-np.sqrt(3.0)/3.0,-2.0+np.sqrt(3.0)/3.0,1000)
def cubic_func(x):
return x**3+6*x**2+11*x+6
fig = plt.figure(figsize=(10,10))
ax = fig.gca()
ax.plot(xzoom,0*xzoom,'--k',linewidth=2.0)
ax.plot(xzoom,cubic_func(xzoom),linewidth=2.0)
ax.set_xlabel(r'$x$',fontsize=22)
ax.set_ylabel(r'$f(x)$',fontsize=22)
ax.set_title(r'$f(x)=x^3+6x^2+11x+6$',fontsize=22)
ax.set_ylim([-.5,.5])
plt.show()
Let's call $a=-2.3$ and $b=-1.5$. Looking at the above figure, we can see that between $f(a)$ and $f(b)$, $f=0$ (i.e. the blue curve intersects the black dotted line). Thus, we have also satisfied the second condition for applying the bisection method to find the roots of our function. Now that we've shown that conditions for applying the bisection method are satisfied, how exactly does the method work?
Calculate $f(a)\times f(c)$,
Evalute $f(c)$.
Now let's write a function that performs the bisection method which accepts the bounds on our interval $a,b$ and our function $f$. Note that we can write our bisection method function in such a way that it is not specific to the kind of function to which we apply it. This will allow us to reuse this piece of code for several different functions.
In [36]:
def bisection(a,b,f,options):
"""Bisection method for interval (a,b) as applied to a monotonically increasing (or decreasing) function f on (a,b)"""
#First guess for root value
c = a + (b - a)/2
#Evaluate function at first guess
f_guess = f(c)
#Initialize the counter
count = 0
#Initialize lists of endpoints and roots
a_list = [a]
b_list = [b]
c_list = [c]
#Begin iteration
while(count<options['max_count'] and np.fabs(f_guess)>options['tol']):
#Output iteration and endpoints
print "Iteration %d, Endpoints (%f,%f)"%(count,a,b)
#Check relative sign
if f_guess*f(a)>0:
#advance left endpoint
a=c
elif f_guess*f(a)<0:
#advance right endpoint
b=c
else:
#case where f(c)==0
print "f(c) == 0 exactly. Breaking out of loop."
break
#update root value guess
c = a + (b-a)/2
#evaluate function at current guess
f_guess = f(c)
#save a,b,c values
a_list.append(a)
b_list.append(b)
c_list.append(c)
#update the counter
count = count + 1
#print final zero approximation
print "Value of the function at approximate zero is f(c=%f)=%f"%(c,f_guess)
#return lists of a,b,c
return a_list,b_list,c_list
Now that we have an algorithm for our bisection method, let's apply it to the above function for our chosen endpoints $a=-2.3$ and $b=-1.5$.
In [37]:
opt = {'max_count':100,'tol':1.0e-10}
alist,blist,clist = bisection(-2.5,-1.9,cubic_func,opt)
Thus, we see that our bisection method performs very well for our cubic function on the given interval. This provides a nice check of our method as we can analytically verify (by hand) that our root on this interval is $x=-2$. Now let's apply this method over the two remaining intervals to find our other zeros. For the leftmost interval, we'll choose $a=-5,b=-2-\sqrt{3}/3$ and for the rightmost interval, we'll choose $a=-2+\sqrt{3}/3,b=4$.
In [38]:
opt['tol'] = 1.0e-6
aleft,bleft,cleft=bisection(-5.0,-2.0-np.sqrt(3.0)/3.0,cubic_func,opt)
In [39]:
aright,bright,cright=bisection(-2.0+np.sqrt(3.0)/3.0,4.0,cubic_func,opt)
Our bisection method is able to calculate the zeros of our cubic function in a reasonable number of iterations. Now that we've shown that our bisection algorithm works, let's apply it to a more useful case: an unknown function where we can't use a simple method like the quadratic formula or factorization to find the zeros.
In [68]:
%run make_unknown_function
In [69]:
unknowns = np.loadtxt('example_data/root_finding_unknown_function.txt')
x_unknown = unknowns[0,:]
y_unknown = unknowns[1,:]
Now let's plot our unknown function to determine where it is monotonically increasing/decreasing. We'll make smart choices about our intervals by "bracketing" our function at $(0.5,2.9)$ and $(3.5,6)$. Note that we are free to make any choice when it comes to our intervals provided that on $(a,b)$, $f$ is either monotonically increasing or decreasing.
In [70]:
a1 = 0.5
b1 = 2.9
a2 = 3.5
b2 = 6.0
In [71]:
fig = plt.figure(figsize=(10,10))
ax = fig.gca()
ax.plot(x_unknown,y_unknown,linewidth=2.0)
ax.plot(x_unknown,0*x_unknown,'--k',linewidth=2.0)
ax.plot([a1,a2],[y_unknown[np.where(x_unknown<a1)[0][-1]],y_unknown[np.where(x_unknown<a2)[0][-1]]],'ko',marker=r'$[$',markersize=25)
ax.plot([b1,b2],[y_unknown[np.where(x_unknown<b1)[0][-1]],y_unknown[np.where(x_unknown<b2)[0][-1]]],'ko',marker=r'$]$',markersize=25)
ax.set_xlabel(r'$x$',fontsize=22)
ax.set_ylabel(r'$f(x)$',fontsize=22)
ax.set_title(r'Unknown Function',fontsize=22)
ax.set_xlim([x_unknown[0],x_unknown[-1]])
plt.show()
Now, we can apply our bisection method to this new unknown function to try to find the roots of this equation. Note that here we cannot use our quadratic method or factorization; the form of the function has been defined when we ran the script make_unknown_function.py. It is only important to know the form of the function so that we can feed it to the bisection method; suffice it to say that is complicated.
In [73]:
a1list,b1list,c1list = bisection(a1,b1,sine_integral_func,opt)
a2list,b2list,c2list = bisection(a2,b2,sine_integral_func,opt)
Now, let's plot the zeros on top of our plot of the unknown function to get a sense of how good our guess is.
In [75]:
fig = plt.figure(figsize=(10,10))
ax = fig.gca()
ax.plot(x_unknown,y_unknown,linewidth=2.0)
ax.plot(x_unknown,0*x_unknown,'--k',linewidth=2.0)
ax.plot([a1,a2],[y_unknown[np.where(x_unknown<a1)[0][-1]],y_unknown[np.where(x_unknown<a2)[0][-1]]],'ko',marker=r'$[$',markersize=25)
ax.plot([b1,b2],[y_unknown[np.where(x_unknown<b1)[0][-1]],y_unknown[np.where(x_unknown<b2)[0][-1]]],'ko',marker=r'$]$',markersize=25)
ax.plot([c1list[-1],c2list[-1]],[sine_integral_func(c1list[-1]),sine_integral_func(c2list[-1])],'ro',markersize=10)
ax.set_xlabel(r'$x$',fontsize=22)
ax.set_ylabel(r'$f(x)$',fontsize=22)
ax.set_title(r'Unknown Function',fontsize=22)
ax.set_xlim([x_unknown[0],x_unknown[-1]])
plt.show()
In [ ]: