Newton’s Method.

We consider Newton’s method for attempting to calculate digits of (Sqrt[2])as the root of f(x) = (x^(2 ))-2 at x = (Sqrt[2]). 

Remember first that the Newton iteration is 

(Subscript[x, i+1])= (Subscript[x, i])+((f((Subscript[x, i])))/(f'((Subscript[x, i])))) = g((Subscript[x, i]))

so the first step is to code the Newton iteration g(x).

g[x_] := x - (((x^2 -2))/(2x))

g[x]

Iteration and Nest.

Now that we’ve coded the Newton iteration, we need to apply it. In most computer languages, one would use a loop. Mathematica has constructs which are more precisely suited for this kind of iteration.

Block[{g,x},
Nest[g,x,6]
]

and

Block[{g,x},
NestList[g,x,6]
]

So let’s try this for the Newton iteration we’ve defined.

N[NestList[g,2,6],50] // ColumnForm

N[(Sqrt[2]),50]

Counting correct digits.

That’s rather impressive! We now want to measure the number of correct digits at each iteration step and see how fast the method is converging.

We will use a useful trick (that you’ll see over and over again). The exponent of a normalized floating point decimal number is given by its logarithm in base 10. So if we’re measuring relative error, we should divide

(error/(true value))=(((Subscript[mantissa, error])×(10^(Subscript[exponent, error])))/((Subscript[mantissa, true value])× (10^(Subscript[exponent, true value]))))≃ (10^((Subscript[exponent, error])-(Subscript[exponent, true value])))

since the mantissas are both between 1 and 2. The number of correct digits is then given by

CorrectDigits[approximation_,truevalue_] := Block[{$MaxExtraPrecision=10000},N[Abs[Log10[Abs[(((N[approximation,10000] - truevalue))/truevalue)]]]]];

though we note that we have to do the computation in an absurd 10,000 digits of accuracy to be certain of the results!

CorrectDigits[3.1415926,π]

Graphing the number of correct digits in Newton’s method

We can now look at this more systematically:

NewtonMethodData = MapIndexed[{#2[[1]],CorrectDigits[#1,(Sqrt[2])]}&, NestList[g,2,6]];
NewtonMethodData // TableForm

That looks like a pattern!

ListLogPlot[NewtonMethodData,PlotMarkers->{Automatic, Medium}]

Interpreting the LogPlot

Now the thing to remember about Log plots is that the LogPlot of a function in the form 

f(n) = a (e^kn)= a ((((e^k)))^n)

is the plot of

log f(n) = k n + log(a)

which is to say, a linear function. We can recover the constants k and a by fitting the data:

a ((Exp[k ])^n) /. FindFit[NewtonMethodData,a Exp[k n],{k,a},n]

which is to say, the number of correct digits seems to go up by a factor of about 2 at each step!

An instructive example.

Of course, it’s completely reasonable to wonder whether this always works! So let’s try 
two more instructive examples. First, let’s try to find the root of

sin(x) = 0

at x = π.

g[x_]:= x - ((Sin[x])/(Cos[x]))

NewtonMethodData = MapIndexed[{#2[[1]],CorrectDigits[#1,π]}&, NestList[g,3,6]];
NewtonMethodData // TableForm

ListLogPlot[NewtonMethodData,PlotMarkers->{Automatic, Medium}]

a ((Exp[k ])^n) /. FindFit[NewtonMethodData,a Exp[k n],{k,a},n]

Another instructive example

Now let’s consider the root of the function

f(x) = (((x-1))^3)

at x=1.

g[x_]:= x - ((((x-1))^3)/(3(((x-1))^2)))

NestList[g,2,6]

NewtonMethodData = MapIndexed[{#2[[1]],CorrectDigits[#1,1]}&, NestList[g,2,20]];
NewtonMethodData // TableForm

NewtonMethodData = MapIndexed[{#2[[1]],CorrectDigits[#1,1]}&, NestList[g,2,20]];

ListPlot[NewtonMethodData,PlotMarkers->{Automatic, Medium}]

m n + b /. FindFit[NewtonMethodData,m n +b,{m,b},n]