<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Add some physical examples to motivate root finding. Just using algebraic examples up to this point.
Solving Algebraic Equations
Step1:  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.
Step2:  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.
Bisection Method
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
Step4:  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?
Input interval endpoints $a,b$ for a function $f$ which satisfies the above conditions.
Calculate first guess of root $c=a+(b-a)/2$ and evaluate $f(c)$.
Calculate $f(a)\times f(c)$,
If $f(a)\times f(c)>0$ (i.e. $c$ is on the same side of the root as $a$), set $a=c$ and update $c$.
If $f(a)\times f(c)<0$ (i.e. $c$ is on the same side of the root as $b$), set $b=c$ and update $c$.
Evalute $f(c)$.
Repeat steps 3 and 4 for a finite number of steps or until $f(c)<\epsilon$, where $\epsilon$ is some specified error.
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.
Step5:  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$.
Step6:  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$.
Step7:  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
Step8:  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.
Step9:  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.
Step10:  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.
<ASSISTANT_TASK:>
Python Code:
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()
Explanation: Add some physical examples to motivate root finding. Just using algebraic examples up to this point.
Solving Algebraic Equations: Root-finding Methods
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$.
End of explanation
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()
Explanation: 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.
End of explanation
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()
Explanation: 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.
Bisection Method
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:
the function $f$ must be monotonically increasing (or decreasing) between $x=a$ and $x=b$, 
between $f(x=a)$ and $f(x=b)$, $f(x)=0$ for some value of $x=x_0$ such that $a<x_0<b$.
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]$.
End of explanation
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
Explanation: 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?
Input interval endpoints $a,b$ for a function $f$ which satisfies the above conditions.
Calculate first guess of root $c=a+(b-a)/2$ and evaluate $f(c)$.
Calculate $f(a)\times f(c)$,
If $f(a)\times f(c)>0$ (i.e. $c$ is on the same side of the root as $a$), set $a=c$ and update $c$.
If $f(a)\times f(c)<0$ (i.e. $c$ is on the same side of the root as $b$), set $b=c$ and update $c$.
Evalute $f(c)$.
Repeat steps 3 and 4 for a finite number of steps or until $f(c)<\epsilon$, where $\epsilon$ is some specified error.
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.
End of explanation
opt = {'max_count':100,'tol':1.0e-10}
alist,blist,clist = bisection(-2.5,-1.9,cubic_func,opt)
Explanation: 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$.
End of explanation
opt['tol'] = 1.0e-6
aleft,bleft,cleft=bisection(-5.0,-2.0-np.sqrt(3.0)/3.0,cubic_func,opt)
aright,bright,cright=bisection(-2.0+np.sqrt(3.0)/3.0,4.0,cubic_func,opt)
Explanation: 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$.
End of explanation
%run make_unknown_function
unknowns = np.loadtxt('example_data/root_finding_unknown_function.txt')
x_unknown = unknowns[0,:]
y_unknown = unknowns[1,:]
Explanation: 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.
End of explanation
a1 = 0.5
b1 = 2.9
a2 = 3.5
b2 = 6.0
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()
Explanation: 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.
End of explanation
a1list,b1list,c1list = bisection(a1,b1,sine_integral_func,opt)
a2list,b2list,c2list = bisection(a2,b2,sine_integral_func,opt)
Explanation: 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.
End of explanation
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()
Explanation: 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.
End of explanation
<END_TASK>