Newton’s Method

Andy Long, Spring 2025

The Newton Raphson method is designed to make convergence as quick as possible among  so-called “fixed point methods” for finding roots (that is, values of x for which f(x)=0).

The central problem is that we require a derivative in order to proceed. But if we have the derivative, then the method approximates the function f with its tangent line at x=a, and finds its root instead (which is an easy problem); then it supposes that it’s gotten closer to the true root, and does it again.... In principle, it could do it forever, producing a sequence that will converge to the true root.

The following function will illustrate this process with a given function f, on an interval [left,right]:

NewtonShow[f_,{left_,right_}]:=
Manipulate[
Module[{a1=a-f[a]/f'[a],y0=f[a]},
Show[
Plot[{f[x],f[a]+f'[a](x-a)},{x,left,right},PlotLabel->{a1,0}],
ListPlot[{{a,y0},{a1,0}},PlotStyle->{PointSize->Large,Red}]
]
],
{a,left,right}
]

For example: let’s check out a root of sin(x) on the interval [-Pi/2,Pi/2] (shrunk slightly for reasons I may explain...).

g[x_]:=Sin[x]
NewtonShow[g,{-Pi/2.1,Pi/2.1}]

In truth, however, the method may fail for a variety of reasons. Let’s look at some.

NewtonShow[g,{-Pi/2,Pi/2}]

Here is an improved (but much more complicated) version of the algorithm. It stores the sequence of iterates, and produces  some nice output.

NewtonsMethod[f_,a0_,eps_,nits_] :=
Module[
{a=N[a0], (* Use local variables (so as not to alter passed arguments) *)
fa = f[N[a0]],
i=0,
iterates, xvalues, yvalues,ymin,ymax,
outputSequences (* Let's collect the output, since we've gone to the trouble of computing it. *)
},
outputSequences={{i,a,fa}}; (* Sequences of iteration, x-value, and y-value *)
While[Abs[fa] > eps && i<=nits,
a = a-(fa/(f'[a]));
fa=f[a];
i=i+1;
outputSequences=Append[outputSequences,{i,a,fa}];  (* once again, collect output. *)
]; 
(* Now let's augment the naive function to include printing some nice formatted output: *)
Print[NumberForm[TableForm[outputSequences,
TableHeadings->{None,{"i","(Subscript[x, k])","f[(Subscript[x, k])]"}}],16]]; 
Print["  x  = ",NumberForm[a,16] ]; 
Print["f(x) = ",NumberForm[fa,16] ];
outputSequences=Transpose[outputSequences];
xvalues=outputSequences[[2]];
yvalues=outputSequences[[3]];
ymin=Min[yvalues];ymax=Max[yvalues];
Print[
Show[
Plot[f[x],{x,Min[xvalues],Max[xvalues]},
PlotRange->{ymin-.1(ymax-ymin),ymax+.1(ymax-ymin)}],
ListPlot[Transpose[{xvalues,yvalues}],PlotStyle->{PointSize-> Large,Red}]
]
];
(* Now I'll actually return the output, so that it can be plotted, etc. *)
Return[xvalues]; 
]

This function’s roots lead us to the Babylonian square root algorithm:

The quadratic  (x^2)-2   allows Newton’s method to converge to a square root of 2.

f[x_]=(x^2)-2;
Plot[f[x],{x,-2,2}]

This one converges to a root for any starting value (except for one particular starting value -- can you tell which initial value of x will cause trouble?): do you see why, and how it will converge for all other values?

Clear[a]
Simplify[a-f[a]/f'[a]]

Perhaps it’s clear now which value will cause trouble?:)

x0=.2;
eps=10^(-6);
nits=100;
NewtonsMethod[f,x0,eps,nits];

NewtonShow[f,{x0,5}]

NewtonsMethod[f,-x0,eps,nits];
Sqrt[2.0]

We can use Newton’s method for all kinds of interesting jobs. For example, Division -- without division! Well, it looks like there’s division --

f[x_]=1/x-2; (* The root of this function is x=1/2 *)
f[x_]=1/x-3; (* The root of this function is x=1/3 *)

but if you work out a = a-((f[a])/(f'[a])) for a quotient like 1/n, the expression loses the quotient, and nothing is left but a product:

Clear[a,n]
f[x_]=1/x-n; (* The root of this function is x=1/n *)
Simplify[a-f[a]/f'[a]]

So division can be done in a computer with nothing but addition and multiplication! If you have a quotient like m/n, then multiply m times the output of Newton’s method for 1/n:

Clear[a,m,n]
f[x_]=1/x-n; 
m=34;n=21;
x0=.05
m*NewtonsMethod[f,x0,eps,nits]

n=3
Plot[f[x],{x,-2,2}]
x0=.2;
eps=10^(-6);
nits=10;
NewtonsMethod[f,x0,eps,nits];

Question: what function would you use to find the quotient 89/55? Try it, and see if you get something near the golden mean.

f[x_]=1/x-55; 
x0=.02
89*NewtonsMethod[f,x0,eps,nits]
N[GoldenRatio]

But sometimes it doesn’t converge....

Once again, the starting value is very important:

f[x_]=1/x-3; 
x0=2;
NewtonsMethod[f,x0,eps,nits];
NewtonShow[f,{-3,3}]

Now try these problems from your text:

#54:

f[x_]:=x^2-2
x0=1;

#55:

f[x_]:=(x-1)^2-2
x0=2;

#56:

f[x_]:=E^x-2
x0=1;

#57:

f[x_]:=Log[x]-1
x0=2;