Shared09 - Newton's Method Assignment / Newton's Method Notes.sagewsOpen in CoCalc
This material was developed by Aaron Tresham at the University of Hawaii at Hilo and is Creative Commons License
licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Prerequisites:

  • Intro to Sage
  • Graphing and Solving Equations
  • Tangent Lines
  • Differentiation

Newton's Method

Newton's Method is a tool that uses calculus to estimate the roots (zeros or x-intercepts) of a function.

In a previous lab, we considered the find_root and solve commands in Sage, which can be used to find roots. The purpose of this lab is to give us some idea of the math behind these commands. There is no practical reason to use Newton's Method (unless you really need to find a root without a computer/calculator). However, it is instructive to see how calculus (which is not directly related to roots) can be used to estimate roots.

Newton's basic idea is that a differentiable function is "close" to its tangent lines (at least near the point of tangency). So if we know the root of a tangent line (easy algebra), we suppose that the root of the function is somewhere close.

Example 1

Estimate the roots of f(x)=x2+2x1f(x)=x^2+2x-1. [Of course, we can find the roots of a quadratic function with the Quadratic Formula, but it's just an example.] Let's look at a graph:

f(x)=x^2+2*x-1
plot(f,xmin=-5,xmax=5)

From this graph, it appears that there is a root near x=12x=\frac{1}{2}.

Let's find the tangent line at this point. Remember, an equation for the tangent line at x=x0x=x_0 is y=f(x0)+f(x0)(xx0)y=f(x_0)+f'(x_0)(x-x_0).

f(x)=x^2+2*x-1
df(x)=derivative(f,x)        #find the derivative
TL(x)=f(1/2)+df(1/2)*(x-1/2) #find the tangent line
print 'The tangent line is y =',TL(x)
plot(f,xmin=0,xmax=1)+plot(TL,xmin=0,xmax=1,color='red')+point((1/2,f(1/2)),color='black',size=25)
The tangent line is y = 3*x - 5/4

Notice that the tangent line and the graph of ff are close together near x=12x=\frac{1}{2}, and we can see that the root of the tangent line is close to the root of ff.

Sage gives the equation of the tangent line as y=3x54y=3x-\frac{5}{4}. To find the root of this line, just plug in 0 for yy and solve for xx:

0=3x540=3x-\frac{5}{4}

54=3x\frac{5}{4}=3x

543=x\frac{\frac{5}{4}}{3}=x

x=5120.41666x=\frac{5}{12}\approx 0.41666

Let's zoom in on the graph and plot this root.

f(x)=x^2+2*x-1
df(x)=derivative(f,x)
TL(x)=f(1/2)+df(1/2)*(x-1/2)
plot(f,xmin=0.4,xmax=.5)+plot(TL,xmin=0.4,xmax=.5,color='red')+point((5/12,0),size=25,color='black')+point((1/2,f(1/2)),color='black',size=25)

On this window, we can see that the root of ff is not exactly 512\frac{5}{12}, but this is a good approximation.

If we want a better approximation, we can repeat the process, starting with a point of tangency that's closer to the root. Remember we started this example using x=12x=\frac{1}{2} as our point of tangency. This time, we'll use x=512x=\frac{5}{12}, the approximation we got last time. We can see from the graph that this is closer to the root of ff than 12\frac{1}{2}, so the tangent line we get from x=512x=\frac{5}{12} should be even closer to the actual root of ff.

So let's find a new tangent line:

f(x)=x^2+2*x-1
df(x)=derivative(f,x)
TL(x)=f(5/12)+df(5/12)*(x-5/12)
print 'The tangent line is y =',TL(x)
plot(f,xmin=.41,xmax=.42)+plot(TL,xmin=.41,xmax=.42,color='red',linestyle='--')+point((5/12,f(5/12)),size=25,color='black')
The tangent line is y = 17/6*x - 169/144

We can see that the tangent line and ff are very close together on this window, so the root of the tangent line is very close to the root of ff.

Let's find the root of the tangent line:

0=176x1691440=\frac{17}{6}x-\frac{169}{144}

169144=176x\frac{169}{144}=\frac{17}{6}x

169144176=x\frac{\frac{169}{144}}{\frac{17}{6}}=x

x=1694080.414216x=\frac{169}{408}\approx0.414216

Just for comparison, let's find the exact root of ff.

f(x)=x^2+2*x-1
solve(f(x)==0,x)
[x == -sqrt(2) - 1, x == sqrt(2) - 1]

There are two solutions, x=21x=-\sqrt{2}-1 and x=21x=\sqrt{2}-1. We have been trying to estimate the postive root, which is x=210.414214x=\sqrt{2}-1\approx0.414214.

This is very close to our estimate, which is off by only 169408(21)0.0000021239\frac{169}{408}-(\sqrt{2}-1)\approx0.0000021239 (about 0.0005%0.0005\% error).

If we wanted to improve our estimate even more, then we could repeat the process again, using x=169408x=\frac{169}{408} as our new point of tangency.

I hope you see from this short example that this process involves repetition of the same steps over and over with different numbers. It lends itself very easily to be made into an algorithm that can be automated.

Before we discuss this algorithm, I have two comments:

  • We have to start the process with a point of tangency (we used 12\frac{1}{2} in our example above). The closer this point is to the actual root, the better. In fact, if our point of tangency is chosen poorly, Newton's Method might never get close to the actual root, or it might take more repetitions than we want to use.

  • Newton's Method only estimates one root at a time. If we had chosen a different point of tangency, we may have ended up at the other root of ff (in our example, ff has two roots).

Now let's develop the algorithm.

Generalizing Newton's Method

Given a function ff with at least one root, the first step is to get a ballpark estimate for the root to use as the x-coordinate of our first point of tangency.

We'll use a graph to get this initial estimate, which we'll call x0x_0.

Now we find the tangent line to ff when x=x0x=x_0:

TL(x)=f(x0)+f(x0)(xx0)TL(x)=f(x_0)+f'(x_0)(x-x_0)

Then we find the root of the tangent line by setting it equal to 0 and solving for xx:

0=f(x0)+f(x0)(xx0)0=f(x_0)+f'(x_0)(x-x_0)

First, distribute:

0=f(x0)+f(x0)xf(x0)x00=f(x_0)+f'(x_0)\cdot x-f'(x_0)\cdot x_0

Then gather like terms:

f(x0)x0f(x0)=f(x0)xf'(x_0)\cdot x_0-f(x_0)=f'(x_0)\cdot x

Now divide both sides by f(x0)f'(x_0):

f(x0)x0f(x0)f(x0)=x\frac{f'(x_0)\cdot x_0-f(x_0)}{f'(x_0)}=x

And simplify:

x0f(x0)f(x0)=xx_0-\frac{f(x_0)}{f'(x_0)}=x

This is the root of the tangent line, so this becomes our new estimate for the root of ff. Let's call it x1x_1. Thus,

x1=x0f(x0)f(x0)x_1=x_0-\frac{f(x_0)}{f'(x_0)}

To improve our approximation, use x=x1x=x_1 as our new point of tangency and repeat the process. The algebra is exactly the same. Simply replace all the x0x_0 in the above equations with x1x_1.

What we get is a new approximation, x2x_2, given by

x2=x1f(x1)f(x1)x_2=x_1-\frac{f(x_1)}{f'(x_1)}

In general, if we have repeated this process nn times to get an estimate xnx_n, then repeating the process again will give a new estimate, xn+1x_{n+1} given by:

xn+1=xnf(xn)f(xn)x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

Example 2

Now we're ready to lay out an algorithm for Newton's Method.

As I go through the steps, I'll use the example of f(x)=x3+2x1f(x)=x^3+2x-1.

Step 1 Define the function

f(x)=x^3+2*x-1

Step 2 Choose an initial value by graphing the function

plot(f,xmin=-5,xmax=5,ymin=-10,ymax=10)

From this graph, we can see that ff has one (real) root, somewhere in the neighborhood of x=12x=\frac{1}{2}.

We'll use this as our initial guess. I'm going to abandon subscripts to make it easier to do this in Sage. Instead, I'll call our estimate "rr," and we'll replace rr with better approximations as we go:

r=1/2  #initial guess

Step 3 Define the number of repetitions

Now we have to decide how many times we want to repeat the process. I'll call this number "nn."

n=10    #number of repetitions

Step 4 Define the derivative

df(x)=derivative(f,x)    #find the derivative

Step 5 Iterate Newton's Method

for i in range(n):
    r=N(r-f(r)/df(r),digits=30)  #I will use N() to convert to decimal - it makes it run faster
    r
0.454545454545454545454545454545 0.453398336679093776885574992831 0.453397651516647791761892775232 0.453397651516403767644746569954 0.453397651516403767644746539000 0.453397651516403767644746539000 0.453397651516403767644746539000 0.453397651516403767644746539000 0.453397651516403767644746539000 0.453397651516403767644746539000

Our estimates get better as we go down the list, so our best estimate for the root of ff is x0.453397651516404x\approx0.453397651516404 (rounded to 15 decimal places).

For comparison, the exact solution is (118593+12)1/323(118593+12)1/30.453397651516404\left(\frac{1}{18}\sqrt{59}\sqrt{3} + \frac{1}{2}\right)^{1/3} - \frac{\frac{2}{3}}{\left(\frac{1}{18}\sqrt{59}\sqrt{3} + \frac{1}{2}\right)^{1/3}}\approx0.453397651516404

To 15 decimal places, our estimate matches the actual root exactly - and that's after only 5 repetitions.

Note 1 It is possible to use Newton's Method to find complex roots as well (just start with a complex guess), but we will not do this in this class.

Note 2 As mentioned above, Newton's Method does not always work. An obvious example is when our initial guess is a local minimum or maximum. Then we get a horizontal tangent line, which has no roots. Another problem arises if the function has a limited domain; in this case, the root of the tangent line may end up outside the domain of the function.

Example 3

Here's an example of a poor choice of initial guess. I have chosen an initial guess of x=0.05x=0.05 which is close to the minimum at x=0x=0. Although the root of the function is around 1.51.5, the root of the tangent line is about 2020.

f(x)=x^2-2
df(x)=derivative(f,x)
TL(x)=f(.05)+df(.05)*(x-.05)
plot(f,xmin=-5,xmax=10,ymax=10)+plot(TL,xmin=-5,xmax=20,color='red')

In this case, Newton's Method will eventually get back to the root, but it requires a few more repetitions.

f(x)=x^2-2
df(x)=derivative(f,x)
r=.05
n=10
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
20.0250000000000021316282072803 10.0624375780274667284134369772 5.13059828749466002723375156069 2.76020818631209362823054307491 1.74239560615097680133916978411 1.44512027881510583654836619356 1.41454406258654789561762947488 1.41421360098284711712728104594 1.41421356237309557584829389812 1.41421356237309504880168872421
N(sqrt(2),digits=30) #This is the exact answer for comparison.
1.41421356237309504880168872421

Example 4

Here's another poor choice of initial guess. Consider the function f(x)=ln(3x)f(x)=\ln(3x). I'm going to choose an initial guess of x=1x=1.

Below is a graph of ff with the tangent line at x=1x=1.

Notice that the root of the tangent line is a negative number, which is not in the domain of ff. That means we can't find a new tangent line using that root for the point of tangency.

f(x)=ln(3*x)
df(x)=derivative(f,x)
TL(x)=f(1)+df(1)*(x-1)
plot(f,xmin=0,xmax=3)+plot(TL,xmin=-1,xmax=3,color='red')

If we try Newton's Method, we end up with complex numbers!

f(x)=ln(3*x)
df(x)=derivative(f,x)
r=1
n=10
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
-0.0986122886681096913952452369225 -0.218716840163818113193686201622 + 0.309799641633409422971177200538*I 0.486573380029164645618123813404 + 0.747851351286347649636703893786*I 0.750863295538137729101406222465 - 0.472094422407517673683036387351*I 0.281017333228159361240218105002 + 0.411360592167483474474911835095*I 0.567723052665054031979560084438 - 0.0269370880736900299716388852836*I 0.266054602422544386407234133316 + 0.0143538498950625204601308076185*I 0.326421345527530860589818981361 + 0.00322899526648179221542540133414*I 0.333277139958278835460959346320 + 0.0000676076180662915803991932452887*I 0.333333335453865043806626867823 + 1.13977978161361846019880300147e-8*I

Newton's Method in One Cell

You can copy and paste the input below to implement Newton's Method.

Change the function ff, then hit run.

If no root is visible, adjust the plot window and run until you find a root.

Once you see a root, choose an initial guess in the ballpark of the root, change the number rr, and hit run.

Example 5

Find the roots of f(x)=x22f(x)=x^2-2.

f(x)=x^2-2                                 #change this (original function)
plot(f,xmin=-10,xmax=10,ymin=-10,ymax=10)  #you may need to adjust the plot window
r=1.5                                        #change this based on the graph (initial guess)
n=10                                       #number of iterations
df(x)=derivative(f,x)                      #only change the rest if you change the function's name
print ''
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
1.41666666666666674068153497501 1.41421568627450980404962203283 1.41421356237468991062629577120 1.41421356237309504880168962350 1.41421356237309504880168872421 1.41421356237309504880168872421 1.41421356237309504880168872421 1.41421356237309504880168872421 1.41421356237309504880168872421 1.41421356237309504880168872421

We found a root near 1.41421.4142.

Notice how the answers stabilize after a few repetitions. This suggests we have the right answer.

For this example, the graph shows two roots. We have found the positive root; to find the negative root, change the initial guess and run it again.

f(x)=x^2-2                                 #change this (original function)
plot(f,xmin=-10,xmax=10,ymin=-10,ymax=10)  #you may need to adjust the plot window
r=-1.5                                        #change this based on the graph (initial guess)
n=10                                       #number of iterations
df(x)=derivative(f,x)                      #only change the rest if you change the function's name
print ''
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
-1.41666666666666674068153497501 -1.41421568627450980404962203283 -1.41421356237468991062629577120 -1.41421356237309504880168962350 -1.41421356237309504880168872421 -1.41421356237309504880168872421 -1.41421356237309504880168872421 -1.41421356237309504880168872421 -1.41421356237309504880168872421 -1.41421356237309504880168872421

Now we have found the second root, around 1.4142-1.4142.

Using Newton's Method to find a point of intersection

If you have two functions, f1f_1 and f2f_2, and you want to know where they cross, then define a new function f=f1f2f=f_1-f_2 and find the roots of ff. This gives you the x-coordinate of a point of intersection. Plug this into either f1f_1 or f2f_2 to get the corresponding y-coordinate.

Example 6

Find the points of intersection of f1(x)=2sin(x)3f_1(x)=\sqrt[3]{2\sin(x)} and f2(x)=e2xf_2(x)=e^{-2x} for 0xπ0\le x \le \pi.

f1(x)=(2*sin(x))^(1/3)
f2(x)=e^(-2*x)
f(x)=f1(x)-f2(x)       #find the root of f1-f2
plot(f,xmin=0,xmax=pi) #I have adjusted the window
r=.2                   #Based on the graph, one root is close to 0.2
n=10
df(x)=derivative(f,x)
print ''
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
0.174564614046446152487753889929 0.175422261883443201580220134114 0.175423347113948424255917239990 0.175423347115679360373198129103 0.175423347115679360373202532593 0.175423347115679360373202532593 0.175423347115679360373202532593 0.175423347115679360373202532593 0.175423347115679360373202532593 0.175423347115679360373202532593

This tells us that the x-coordinate of one point of intersection is approximately 0.1754233471156790.175423347115679.

To find the y-coordinate, plug this into f1f_1 or f2f_2. I'll plug it into both as a check (I'd better get the same answer).

f1(x)=(2*sin(x))^(1/3)
f2(x)=e^(-2*x)
N(f1(0.175423347115679360373202532593),digits=30)
N(f2(0.175423347115679360373202532593),digits=30)
0.704091686899284448083383801316 0.704091686899284448083383801316

So our point of intersection is approximately (0.175423347115679,0.704091686899284)(0.175423347115679,0.704091686899284).

The second root is hard to approximate with Newton's Method (at least in Sage), since Sage does not allow fractional powers of negative numbers (this happens when xx goes past π\pi). We need a very good initial guess to avoid ending up with complex numbers.

f1(x)=(2*sin(x))^(1/3)
f2(x)=e^(-2*x)
f(x)=f1(x)-f2(x)
plot(f,xmin=0,xmax=pi)
r=3.14159265
n=10
df(x)=derivative(f,x)
print ''
for i in range(n):
    r=N(r-f(r)/df(r),digits=30)
    r
3.14159265034448884600068852662 3.14159265033359929124197408890 3.14159265033358710682062112553 3.14159265033358710680542340382 3.14159265033358710680542340382 3.14159265033358710680542340382 3.14159265033358710680542340382 3.14159265033358710680542340382 3.14159265033358710680542340382 3.14159265033358710680542340382

This is the x-coordinate; now find the y-coordinate:

N(f1(3.14159265033358710680542340382),digits=30)
N(f2(3.14159265033358710680542340382),digits=30)
0.00186744274386954580104325212193 0.00186744274386954580104327322816

So our second point of intersection is approximately (3.141592650333587,0.001867442743870)(3.141592650333587,0.001867442743870).

Here's a graph of the two functions:

f1(x)=(2*sin(x))^(1/3)
f2(x)=e^(-2*x)
plot(f1,xmin=0,xmax=pi)+plot(f2,xmin=0,xmax=pi,color='red')+point((0.175423347115679,0.704091686899284),size=25,color='black')+point((3.141592650333587,0.001867442743870),size=25,color='black')