Newton's Method

In this notebook, we develop Newton's Method for locating a root of an equation.

Newton's Method is an algorithm for approximating solutions to an equation of the form $f(x) = 0$. The idea is that, since the tangent line approximates the function near a point, then the zeros of the tangent line "should" approximate the zeros of the function. The algorithm is:

  1. Make an initial guess, $x = x_0$;

  2. For $i \geq 0$, compute $x_{i + 1} = x_i - \dfrac{f(x_i)}{f'(x_i)}$;

  3. Continue until $\big\vert x_{i+1} - x_i\big\vert$ is less than some prescribed tolerance, $\varepsilon$.

We define the algorithm as a python function that takes in the function $f$, the initial guess $x_0$, and the tolerance $\varepsilon$. It returns the approximate root and the number of iterations done to find it.

In [1]:
def Newton(func,guess,tolerance):
    F = func;
    G = guess;
    tol = tolerance;
    xx = {}
    xx[0]=G.n();
    xx[1] = (xx[0] - F(x = xx[0])/F.diff(x)(x = xx[0])).n();
    for i in range(500): # sets maximum number of iterations at 500. If this is too low, we can change it.
        if (xx[i+1] - xx[i]).abs() >= tol:
            xx[i+2] = (xx[i+1] - F(x = xx[i+1])/F.diff(x)(x = xx[i+1])).n();
        else:
            return xx[i+1], i+1

Let's try it on a familiar example: $f(x) = \cos(x) - x$. First, we look at the graph to help make an inital guess.

In [2]:
var('x');
f1(x) = cos(x);
f2(x) = x;
show(plot(f1,(x,-pi,pi),color='red') + plot(f2,(x,-pi,pi),color='blue'))

Based on the graph, the solution should be between $\tfrac{1}{2}$ and $1$. We can choose either value as our initial guess. Let's choose $x_0 = 1$ and start by setting our tolerance at $\varepsilon = 0.01$.

In [3]:
Newton(cos(x) - x,1,0.01)
Out[3]:
(0.739085133385284, 3)

Let's check this answer against Sage's built-in version of Newton's Method, $\small\texttt{.find root()}$:

In [4]:
(cos(x) - x).find_root(0,1)
Out[4]:
0.7390851332151559
In [5]:
Newton(cos(x)-x,1,0.01)[0] - (cos(x) - x).find_root(0,1)
Out[5]:
1.70128133802905e-10

Our answer is actually much better than the tolerance, in only 3 steps. However, if our initial guess is not such a good one (for example, if we weren't able to look at the graph before guessing), then the method could require many more steps to reach the same tolerance.

In [6]:
Newton(cos(x)-x,5,0.01)
Out[6]:
(0.739085140008201, 109)

If our guess is really, really horrible, then the method does not return any output at all.

In [7]:
Newton(cos(x)-x,500,0.01)
In [8]:
# To do---Maybe we can figure out how to get an error message: 'Make a better initial guess'

There is another important reason to try to make the best initial guess possible. The equation we are trying to solve could have multiple roots. If our goal is to identify a specific root, then a bad initial guess could result in the algorithm finding a different root than that which is desired.

For example, consider the function $f(x) = x^3-x$. We know that the zeros of this function are $0, \pm 1$, but for the sake of example, let's use Newton's Method to try to approximate the middle one.

Our initial guess needs to be between $-1$ and $1$, but if it's too close to either of those, the algorithm could pick out one of them. Again, let's look at the graph to get an idea of what we should guess.

In [9]:
f3(x) = x^3 - x;
graph3 = plot(f3,(x,-2,2));
show(graph3)

The graph appears to have two critical points between $-1$ and $1$. We can look at the derivative of $f$ to get an idea of where they occur.

In [10]:
f4 = f3.diff(x);
graph4 = plot(f4,(x,-2,2),color='red');
show(graph3 + graph4)

The red graph is the derivative of the blue one, so the zeros of the red graph are critical numbers of the blue graph. To ensure that Newton's method produces the desired root, our initial guess should be sufficiently away from those numbers. It looks like $\pm0.4$ should work.

In [11]:
Newton(f3(x),0.4,0.01)
Out[11]:
(1.84212966603849e-12, 4)

This number is 0.0000000000018421..., so it is very close to 0.

We can improve the tolerance to get closer to 0.

In [12]:
Newton(f3(x),0.4,0.00000000000001)
Out[12]:
(0.000000000000000, 6)

Let's see what would happen if our initial guess was $0.6$ instead of $0.4$.

In [13]:
Newton(f3(x),0.6,0.01)
Out[13]:
(1.00011819089665, 8)

As expected, Newton's method found the wrong root!

Notice that while the function $\small\texttt{Newton}$ compiles the entire list $\small\texttt{xx=\{ \}}$ as it's running, it does not waste memory storing it after it's finished.

In [14]:
whos

Variable   Type          Data/Info
----------------------------------
Newton     function      <function Newton at 0x7f28d2220f50>
f1         Expression    x |--> cos(x)
f2         Expression    x |--> x
f3         Expression    x |--> x^3 - x
f4         Expression    x |--> 3*x^2 - 1
graph3     Graphics      Graphics object consistin<...>g of 1 graphics primitive
graph4     Graphics      Graphics object consistin<...>g of 1 graphics primitive

We end these notes by illustrating Newton's Method with a progression of graphical steps. We'll use the function $f(x) = \tfrac{1}{x}-2$ with initial guess $x_0 = 0.8$ and tolerance $\varepsilon = 0.0001$.

In [15]:
f(x) = 1/x - 2;
a = 0.8;
e = 0.0001;

def Newton_Points(function,guess,tolerance):
    F = function;
    G = guess;
    tol = tolerance;
    xx = {}
    xx[0]=G.n();
    xx[1] = (xx[0] - F(x = xx[0])/F.diff(x)(x = xx[0])).n();
    for i in range(500): # sets maximum number of iterations at 500.
        if (xx[i+1] - xx[i]).abs() >= tol:
            xx[i+2] = (xx[i+1] - F(x = xx[i+1])/F.diff(x)(x = xx[i+1])).n();
        else:
            return xx

points = Newton_Points(f(x),a,e);
points
Out[15]:
{0: 0.800000000000000,
 1: 0.320000000000000,
 2: 0.435200000000000,
 3: 0.491601920000000,
 4: 0.499858944504627,
 5: 0.499999960206695,
 6: 0.499999999999997}
In [16]:
len(points)
Out[16]:
7
In [36]:
graph = plot(f,(x,0,2),ymin=-2,ymax=2,color='red',aspect_ratio=1/3);
graph
Out[36]:
In [33]:
p={}
for i in range(7):
    p[i] = plot(point((points[i],f(x)(x=points[i])),color='black',size=15));
    
pp={}
for i in range(7):
    pp[i] = plot(point((points[i],0),color='black',size=15));
    
show(graph + p[0],aspect_ratio=1/3)
In [19]:
tanlines = {}
for i in range(7):
    tanlines[i] = plot(f(x)(x = points[i]) + f.diff(x)(x = points[i])*(x - points[i]),(x,0,2),ymin=-2,ymax=2,thickness=0.5,color='blue');
In [34]:
show(graph + p[0] + tanlines[0],aspect_ratio=1/3)
In [35]:
show(graph + p[0] + tanlines[0] + pp[1],aspect_ratio=1/3)
In [22]:
var('t');
vlines = {}
for i in range(7):
    vlines[i] = parametric_plot([points[i],t*f(points[i])],(t,0,1),xmin=0,xmax=2,ymin=-2,ymax=2,linestyle="--",color='black');
In [28]:
show(graph + p[0] + tanlines[0] + pp[1] + vlines[1] + p[1],aspect_ratio=1/3)
In [29]:
show(graph + p[0] + tanlines[0] + pp[1] + vlines[1] + p[1] + tanlines[1] + pp[2],aspect_ratio=1/3)
In [30]:
show(graph + pp[1] +p[1] + tanlines[1] + pp[2] + vlines[2] + p[2],aspect_ratio=1/3)
In [31]:
show(graph + pp[1] +p[1] + tanlines[1] + pp[2] + vlines[2] + p[2] + tanlines[2] + pp[3],aspect_ratio=1/3
)
In [32]:
show(graph + pp[1] + pp[2] + vlines[2] + p[2] + tanlines[2] + pp[3],aspect_ratio=1/3)

As this process continues the intersection points of the tangent lines with the $x$-axis approach the zero of the function.