4.7 Newton's Method

Mathematica script by Chris Parrish,
   cparrish@sewanee.edu

Problems from James Stewart,
  "Calculus: Concepts and Contexts,"
  Second Edition, Brooks/Cole, 2001

Finding roots

Hughes-Hallett, Gleason, et al, Exercise 5.7.1, page 300

Let's find the real root of the equation x^3+3x^2+3x - 6 = 0.
First, graph the function so that we know what to expect.

In[495]:=

Needs["Graphics`Colors`"] Clear[f, x]  f[x_] := x^3 + 3x^2 + 3x - 6 Plot[f[x], {x, -6, 6}, PlotStyleBlue] ;

[Graphics:HTMLFiles/4.7_Newton_4.gif]

Calculate the derivative of f.

In[499]:=

f'[x]

Plot[f'[x],{x,-6,6},
     PlotStyle->ForestGreen];

Out[499]=

3 + 6 x + 3 x^2

[Graphics:HTMLFiles/4.7_Newton_6.gif]

In[501]:=

Solve[f'[x] == 0,x]

Out[501]=

{{x -1}, {x -1}}

In[502]:=

f'[-1]

Out[502]=

0

Let's let x = 0 be our initial guess, and use Newton to get a better guess.

In[503]:=

newt[x_] = x - f[x]/f'[x]

{0, newt[0], newt[newt[0]], newt[newt[newt[0]]], newt[newt[newt[newt[0]]]]} //N

Out[503]=

x - (-6 + 3 x + 3 x^2 + x^3)/(3 + 6 x + 3 x^2)

Out[504]=

RowBox[{{, RowBox[{0., ,, 2., ,, 1.25926, ,, 0.963308, ,, 0.914213}], }}]

Let's package this up in a Mathematica routine.

Newton's Method

Hughes-Hallett, Gleason, et al, Exercise 5.7.1, page 300

The Newton-Raphson formula indicates how to improve
a reasonable approximation to the root of a function f(x).
The "N[...,12]"  specifies that we wish to use 12-digit accuracy.

In[505]:=

Clear[newtonRaphson,xn]

newtonRaphson[xn_] := N[xn - f[xn]/f'[xn],12]

Now let's apply the method to a specific function f(x), the same one we were working on previously.

In[507]:=

Clear[f,x,x0]

f[x_] := x^3 + 3x^2 + 3x - 6;

x0 = 0;

x1 = newtonRaphson[x0]
x2 = newtonRaphson[x1]
x3 = newtonRaphson[x2]

Out[510]=

2.00000000000

Out[511]=

1.25925925926

Out[512]=

0.9633080182

A built-in Mathematics function called NestList will do this iteration for us automatically.

In[513]:=

nrList = NestList[newtonRaphson,x0,3]

Out[513]=

RowBox[{{, RowBox[{0, ,, 2.00000000000, ,, 1.25925925926, ,, 0.9633080182}], }}]

We can even arrange for these numbers to be displayed in a table.

In[514]:=

nrTable = Transpose[{Range[0, 3], NestList[newtonRaphson, x0, 3]}] ;  TableFor ...                                                                                                 n

Out[515]//TableForm=

n p  n
Newton-Raphson 0 0
1 2.00000000000
2 1.25925925926
3 0.9633080182

Approximating 50^(1/3)

Hughes-Hallett, Gleason, et al, Exercise 5.7.2, page 300

In[516]:=

Clear[f, x] <br /> f[x_] := x^3 - 50 ;  Plot[f[x], {x, 0, 6}, PlotStyleBlue] ;

[Graphics:HTMLFiles/4.7_Newton_19.gif]

In[519]:=

Clear[newtonRaphson, x0, xn] <br /> newtonRaphson[xn_] := N[xn - f[xn]/f '[xn], 12] <br /> x0  ...                                                                                                 n

Out[523]//TableForm=

n p  n
Newton-Raphson 0 4
1 3.70833333333
2 3.68419040806
3 3.6840315055

Solving Exp[-x] = Log[x]

Hughes-Hallett, Gleason, et al, Exercise 5.7.5, page 300

In[524]:=

Clear[g,h,x]

g[x_] := Exp[-x];
h[x_] := Log[x];

Plot[{g[x],h[x]},{x,0.01,4},
     PlotStyle->{Blue,Red}];

[Graphics:HTMLFiles/4.7_Newton_22.gif]

In[528]:=

Clear[f,x]

f[x_] := Exp[-x] - Log[x];

Plot[f[x],{x,0.01,4},
     PlotStyle->{Indigo}];

[Graphics:HTMLFiles/4.7_Newton_23.gif]

In[531]:=

Clear[newtonRaphson, x0, xn] <br /> newtonRaphson[xn_] := N[xn - f[xn]/f '[xn], 12] <br /> x0  ...                                                                                                 n

Out[535]//TableForm=

n p  n
Newton-Raphson 0 1
1 1.26894142137
2 1.30910840327
3 1.30979938867

Solutions of Log[x] = 1/x

Hughes-Hallett, Gleason, et al, Exercise 5.7.9, page 300

In[536]:=

Clear[g,h,x]

g[x_] := Log[x];
h[x_] := 1/x;

Plot[{g[x],h[x]},{x,0.1,3},
     PlotStyle->{Blue,Red}];

[Graphics:HTMLFiles/4.7_Newton_26.gif]

In[540]:=

Clear[f,x]

f[x_] := Log[x] - 1/x;

Plot[f[x],{x,0.1,3},
     PlotStyle->{Indigo}];

[Graphics:HTMLFiles/4.7_Newton_27.gif]

In[543]:=

Clear[newtonRaphson, x0, xn] <br /> newtonRaphson[xn_] := N[xn - f[xn]/f '[xn], 12] <br /> x0  ...                                                                                                 n

Out[548]//TableForm=

n p  n
Newton-Raphson 0 1
1 1.50000000000
2 1.73508140270
3 1.76291539065
4 1.7632227978
5 1.7632228344

Zeros of  x^3 + x - 1

Hughes-Hallett, Gleason, et al, Exercise 5.7.10, page 300

In[549]:=

Clear[f, x] <br /> f[x_] := x^3 + x - 1 ; <br /> Plot[f[x], {x, -2, 2}, PlotStyleIndigo] ;

[Graphics:HTMLFiles/4.7_Newton_32.gif]

In[552]:=

Clear[newtonRaphson, x0, xn] <br /> newtonRaphson[xn_] := N[xn - f[xn]/f '[xn], 12] <br /> x0  ...                                                                                                 n

Out[557]//TableForm=

n p  n
Newton-Raphson 0 1
1 0.750000000000
2 0.68604651163
3 0.68233958260
4 0.68232780395
5 0.6823278038

Let's check on the accuracy of that last result.

In[558]:=

xNR = 0.682328;

f[xNR]

Out[559]=

4.70168*10^-7

In[560]:=

f[xNR - 0.000001]

Out[560]=

RowBox[{-, 1.92654*10^-6}]


Created by Mathematica  (March 16, 2004)