A root finding example

For our first programming "project", you'll try to find the roots of a fun function using SciPy's built in brentq as well as your own Newton's method code. Here's a (somewhat easier) sample problem, where we'll play with $$f(x) = x^4-5x-1.$$ As we often do, we'll start by importing some functionality.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import brentq

Let's begin our analysis by looking at a plot of the function. Since we're looking for roots, let's also plot the line $y=0$ so we can see the roots clearly.

In [2]:
def f(x): return x**4-5*x-1
xs = np.linspace(-2,2,100)
plt.plot(xs,f(xs))
plt.plot([-2,2],[0,0])
Out[2]:
[<matplotlib.lines.Line2D at 0x108e31a90>]

Let's look for the positive root; there's clearly one between $1$ and $2$ and Brent's method can find it.

In [3]:
x0 = brentq(f,1,2)
x0
Out[3]:
1.7720290019648661

Here's a geometric double check.

In [4]:
plt.plot(xs,f(xs))
plt.plot([-2,2],[0,0])
plt.plot(x0,0,'o')
Out[4]:
[<matplotlib.lines.Line2D at 0x108ff4048>]

Cool! Here's how to find it with Newton iteration.

In [5]:
def n(x): return x-(x**4-5*x-1)/(4*x**3-5)
xi = 1.5
cnt = 0
while abs(f(xi))>10**(-15) and cnt < 100:
    xi = n(xi)
    cnt = cnt + 1
xi = n(n(xi))
xi
Out[5]:
1.772029001964867