Newton’s Method
March 18, 2012
Math 151 (Pelsmajer)

Warning: “newton[x_]” is redefined as the lab progresses: sometimes it uses f and sometimes it uses p.  Also, the definition of f will change.
When you do the exercises, you should invent new function names, to be safe.

Introduction

Newton’s Method finds approximate roots of a function, when you start out with a rough idea of where the solutions are located.

This is helpful when you want approximate solutions to an equation.  Say you want to solve A(x)=B(x).  
Solve A(x)-B(x)=0 instead.  (It is the same problem).  
In other words, let f(x)=A(x)-B(x) and find the zeroes (or "roots") of the function f(x)... which you can do using Newton's Method.

The idea behind Newton’s Method is that if we make a fairly good initial guess (say (Subscript[x, 0])) to the zero of a function f(x), then the x-intercept of the tangent line to the graph of f at ((Subscript[x, 0]),f ((Subscript[x, 0]))) called (Subscript[x, 1]), is usually a better estimate of the zero than was (Subscript[x, 0]).

Run the following, then drag the slider to the right to demonstrate:

f[x_]:=(x^3)-4*x-1;
x0=3;
frames=ConstantArray[0,2];
thecurve=Plot[f[x],{x,2,4},PlotRange->{{2,3.1},{-5,15}},PlotStyle->Green];
firstguess=Graphics[{PointSize[Medium],Blue,Point[{x0,0}]}];
frames[[1]]=Show[thecurve,firstguess];
xline=Graphics[{Dashed,Line[{{x0,0},{x0,f[x0]}}]}];
pointoncurve=Graphics[{PointSize[Medium],Green,Point[{x0,f[x0]}]}];
tangent=p'[x0]*(x-x0)+f[x0];
tplot=Plot[tangent,{x,2,4},PlotRange->{{2,3.1},{-5,15}},PlotStyle->Red];
frames[[2]]=Show[frames[[1]],xline,pointoncurve,tplot];

Manipulate[frames[[n]],{n,1,2,1},FrameLabel->"Newton's Method"]

In the plot above, there is a zero near (Subscript[x, 0])=2.1,  (Subscript[x, 0])=3 is our initial guess, and (Subscript[x, 1]) is roughly 2.4.

Repeating the trick starting with (Subscript[x, 1]), we find an even better estimate of the zero, (Subscript[x, 2]).
Repeating it with (Subscript[x, 2]), we find (Subscript[x, 3]), etc.

Soon we will see that when (Subscript[x, n]) is the n th approximate root, the (n+1)st approximation is given by (Subscript[x, n+1])=(Subscript[x, n])-((f((Subscript[x, n])))/(f'((Subscript[x, n])))).  Which allows us to do this:

Run the following cells, and drag the slider to the right to show the successive approximations found by Newton’s Method:

newton[x_]:=x-f[x]/f'[x]

f[x_]:=(x^3)-4*x-1
xn=3;
numsteps=4;
frames=ConstantArray[0,numsteps];
thecurve=Plot[f[x],{x,2,4},PlotRange->{{2,3.1},{-5,15}},PlotStyle->Green];
firstguess=Graphics[{PointSize[Medium],Blue,Point[{xn,0}]}];
frames[[1]]=Show[thecurve,firstguess];
For[n=2,n≤numsteps,n++,
xline=Graphics[{Dashed,Line[{{xn,0},{xn,f[xn]}}]}];
xnpointoncurve=Graphics[{PointSize[Medium],Green,Point[{xn,f[xn]}]}];
tangent=f'[xn]*(x-xn)+f[xn];
tplot=Plot[tangent,{x,2,4},PlotRange->{{2,3.1},{-5,15}},PlotStyle->Red];
xintercept=Graphics[{PointSize[Medium],Blue,Point[{newton[xn],0}]}];
frames[[n]]=Show[frames[[n-1]],xnpointoncurve,tplot, xintercept];
xn=newton[xn];
]
Manipulate[frames[[n]],{n,1,numsteps,1},FrameLabel->"Newton's Method"]

Where does the formula come from?

The tangent line at (Subscript[x, n]) is y = f'((Subscript[x, n]))(x-(Subscript[x, n]))+f((Subscript[x, n])).

To get the x-intercept, set y=0 and solve for x .

0 = f'((Subscript[x, n]))(x-(Subscript[x, n]))+f((Subscript[x, n]))

-f((Subscript[x, n])) = f'((Subscript[x, n]))(x-(Subscript[x, n]))

-f((Subscript[x, n]))/f'((Subscript[x, n])) =x-(Subscript[x, n])

(Subscript[x, n]) -f((Subscript[x, n]))/f'((Subscript[x, n])) = x

x= (Subscript[x, n]) -f((Subscript[x, n]))/f'((Subscript[x, n])) . 

So we set (Subscript[x, n+1]) = (Subscript[x, n]) -f((Subscript[x, n]))/f'((Subscript[x, n])) .

newton[x_]:=x-f[x]/f'[x]

Examples

Example 1

Find all numbers x in the interval  [0,1]  for which   (([cos(((3 π x)/2))] )^2)=x.

To determine the number of solutions to this problem, we plot a graph of the functions
F(x)=(([cos(((3 π x)/2))] )^2) and G(x)=x for 0<x<1. 
Then the solutions to this problem are all the values of x at which F(x)=G(x) .

funcF[x_]:=(Cos[3*Pi*x/2])^2
funcG[x_]:=x
Plot[{funcF[x], funcG[x]},{x,0,1},PlotStyle->{Green,Blue}]

We can see that there are three points of intersection: one near 0.25, one near 0.5 and one near 0.75.

The solutions are also the zeroes (or roots) of the function
                                              f(x)=G(x)-F(x).

We plot a graph of the functions f, F, and G so we can see indeed that f has zeroes where F and G intersect (near 0.25, near 0.5 and near 0.75).

f[x_]:=funcG[x]-funcF[x]
Plot[{funcF[x], funcG[x], f[x]},{x,0,1},PlotStyle->{Green, Blue,Red}]

In order to estimate these three zeros of f in the interval [0,1] using Newton's method, we begin by defining a function “newton”, which will be used to compute the successive approximations. This function corresponds to the formula for the Newton iteration mentioned in the Introduction.

newton[x_]:=x-f[x]/f'[x]

(When we iterate with this function, it does not remember all intermediate values (Subscript[x, n]), but instead constantly overwrites the value of x as soon as a better approximate solution has been found - see the line “xn=newton[xn]” in the code below.)

The following code iterates Newton's method until we have found an x for which f(x) is less than (10^(-11)) in absolute value.

xn=0.25;    (* our initial guess *)
Print["(Subscript[x, 0])=", PaddedForm[xn,{12,12}],", f((Subscript[x, 0]))=", NumberForm[f[xn],12]];
k=0;             (* a counter *)
While[Abs[f[xn]]>10^(-11),            (* the desired accuracy *)
xn=newton[xn];
k=k+1;
Print[(Subscript["x", k]),"=", PaddedForm[xn,{12,12}],", f(",(Subscript["x", k]),")=", NumberForm[f[xn],12]];
]

(“PaddedForm” is used so that every decimal is written with the same number of decimal places shown.  “NumberForm” is similar but different.  Neither one is so important.)

Optional alternative: There is a more compact---but rather cryptic---Mathematica construct that does the same job, which is based on pure functions and nested lists. (If you are interested, you can find more information on these concepts in the Help.) We format the output in the same way as above using PaddedForm.

NestWhileList[newton[#]&,0.25,Abs[f[#]]>10^(-11)&];
PaddedForm[%,{12,12}]

Since |f((Subscript[x, 4])) | is less than (10^(-11)), it approximates zero with accuracy (10^(-11)).  Thus, the first zero of f in the interval [0,1] has x-value 0.2277581316 , rounded to ten decimal places.

The second zero in the interval [0,1] is approximately x=0.5. 
Actually it is exactly 0.5,  since (([Cos(((3π)/4))])^2)=((((1/(Sqrt[2]))))^2)=(1/2).  Done.

To approximate the final zero of f in [0,1], we make the initial guess (Subscript[x, 0])=0.75 (which we found using the graph of f earlier) then apply Newton's method to obtain the successive approximations.

xn=0.75;    (* our initial guess *)
Print["(Subscript[x, 0])=", PaddedForm[xn,{12,12}],", f((Subscript[x, 0]))=", NumberForm[f[xn],12]];
k=0;             (* a counter *)
While[Abs[f[xn]]>10^(-11),            (* the desired accuracy *)
xn=newton[xn];
k=k+1;
Print[(Subscript["x", k]),"=", PaddedForm[xn,{12,12}],", f(",(Subscript["x", k]),")=", NumberForm[f[xn],12]];
]

You can see that the final zero (rounded to ten decimal places) is x=0.7722418684.

Example 2

Estimate all roots of the third degree polynomial function p(x)=(x^3)-4x-1.

We begin by graphing the polynomial y=p(x)  in the plane.

p[x_]:=(x^3)-4*x-1
Plot[p[x],{x,-3,3},PlotStyle -> Blue]

Clearly, p has exactly three real zeros. 
Let's use Newton's method, as before, to estimate these roots.

We can choose anything as the initial value; let's try (Subscript[x, 0])=1.0.  (You need to use "1.0" rather than "1" because of a technicality.)

As in Example 1, we use a function “newton” to obtain succesive estimates:

newton[x_]:=x-p[x]/p'[x]

xn=1.0;    (* our initial guess *)
Print["(Subscript[x, 0])=", PaddedForm[xn,{12,12}],", p((Subscript[x, 0]))=", NumberForm[p[xn],12]];
k=0;             (* a counter *)
While[Abs[p[xn]]>10^(-10),            (* the desired accuracy *)
xn=newton[xn];
k=k+1;
Print[(Subscript["x", k]),"=", PaddedForm[xn,{12,12}],", p(",(Subscript["x", k]),")=", NumberForm[p[xn],12]];
]

That’s one of the roots (the x-value). 

To understand what’s going on, we graph the approximations and the tangent lines at each pair ((Subscript[x, n]),p((Subscript[x, n]))) for  n=1,2,. . . ,9, so that we can easily see how (Subscript[x, n+1]) is attained graphically from (Subscript[x, n]) at each stage (and slide the slider, as usual, to see it progress).

xn=1.0;
numsteps=5;
frames=ConstantArray[0,numsteps];
thecurve=Plot[p[x],{x,-3,3},PlotRange->{{-3,3},{-15,15}},PlotStyle->Green];
firstguess=Graphics[{PointSize[Medium],Blue,Point[{xn,0}]}];
frames[[1]]=Show[thecurve,firstguess];
For[n=2,n≤numsteps,n++,
xline=Graphics[{Dashed,Line[{{xn,0},{xn,p[xn]}}]}];
xnpointoncurve=Graphics[{PointSize[Medium],Green,Point[{xn,p[xn]}]}];
tangent=p'[xn]*(x-xn)+p[xn];
tplot=Plot[tangent,{x,-3,3},PlotRange->{{-3,3},{-15,15}},PlotStyle->Red];
xintercept=Graphics[{PointSize[Medium],Blue,Point[{newton[xn],0}]}];
frames[[n]]=Show[frames[[n-1]],xnpointoncurve,tplot, xintercept];
xn=newton[xn];
]
Manipulate[frames[[n]],{n,1,numsteps,1},FrameLabel->"Newton's Method"]

There are two more roots to be found.

Exercises

Exercise 0: Finish Example 2

Find the other two zeros.

Exercise 1: Finding a Fixed Point

(a) Make a plot that shows the graphs of y = cos(x) and y = x, for x between 0 and 2.  

At the intersection of these curves, cos(x) = x, which is called a “fixed point”.

(b) Find it, using Newton’s method. (Hint: re-read the beginning of the Introduction of this lab.)

(c) Compute cos(x) at the fixed point (you should get the same point again).

Exercise 2: Calculate  π

The function t(x)=tan x has a zero at x=π.
(a) Use Newton’s method with tan x to obtain an approximation of π, starting with (Subscript[x, 0])=3.
(b) Do it again but obtain 50 decimal places of accuracy.  Ordinarily Mathematica does not use numbers of that precision, so you’ll need a little something extra:

1. Rather than x0=3, start out with 100 decimal places of precision by using the following:

x0=SetPrecision[3,100]

2. The demonstrations quit when they acheived 11 decimal places of accuracy (in the y-value).  You want to change “11” to be 50 or maybe a bit more.
    (Look for the “11” and change it.)
3. So that we can see all that precision, change every 12 in PaddedForm and NumberForm to be 51 or more.

Exercise 3: Calculate e

The natural logarithm ln (x)=1 when x=e (“Euler's number”).
(Mathematica writes ln(x) simply as “Log[x]”.)

Use Newton's method to obtain the value of  e to 50 significant digits of accuracy.