# Newton-Raphson's method for root finding (ALGOL 68)

### From LiteratePrograms

**Other implementations**:**ALGOL 68**| C

The Newton-Raphson method finds approximations of zeroes of differentiable functions. We present the use of this method for efficient calculation of square roots (in fact, it can be easily generalized to any roots) of positive numbers.

We can use Newton's method to find the zero of the function *f*(*x*) = *x*^{2} − *N*. The *x* where *f*(*x*) = 0 is the
square root of *N*. Newton's method starts with a guess
*x*_{0} and iteratively refines this guess. We can control
the quality of the approximation by the number of iterations made;
the fewer we do the faster we get the result but the less accurate it
is. The outline of the algorithm thus looks like this:

<<Newton's method>>=define local variables make an initial guess for root of n; DO refine guess; IF guess sufficient THEN stop iteration FI OD

First we want to look at how the refinement works. Given *x*_{i},
the next *x*_{i + 1} is chosen as the zero in the tangent
at *f*(*x*_{i}). The following figure illustrates this for
the first iteration step in calculating the square root of 2 with an
initial guess of *x*_{0} = 1.5. The red line is the tangent
at *f*(1.5). Where it crosses the x axis is our refined
approximation *x*_{1}. With each iteration, we make a step
towards the true zero of *f*(*x*) = *x*^{2} − 2, marked as
*r*.

Given *x*_{i}, we calculate *x*_{i + 1} by . For *f*(*x*) = *x*^{2} − *N*
we have *f*'(*x*) = 2*x* and thus assuming `xn`

was the
result of our last iteration we calculate our new guess `xn`

as

<<refine guess>>=x := xn; xn := x - (x * x - n) / (2 * x)

This uses variables `x`

and `xn`

which we need to define.

<<define local variables>>=LONG REAL x, xn;

How long do we iterate? This method iterates so fast on square roots so lets just let it converge.

If at any time *x*_{n} = *x*_{n + 1} then it is obvious that we
need not iterate any longer. Our guess is also sufficient then.

<<guess sufficient>>=x = xn

Now all that is left to do is finding a suitable initial guess. We choose
a very simple method and just calculate the values of *f*(*x*)
for all integer numbers . For *n* > 0 these will initially be negative. So we know when we get to
the first *f*(*x*) > 0 then the square root must be between
*x* and *x* − 1. If *f*(*x*) = 0, we have
found the square root and can simply return it.

<<make an initial guess for root of n>>=IF n = 0 OR n = 1 THEN xn := n; stop iteration ELIF n < 0 THEN xn := n/0 ELSE xn := n FI

For taking square roots of very large or very small fractional values
of *n*, we could instead use an approach based on
where we successively double or halve the guess until it changes sign.

We're done and can place our algorithm into a function definition and wrap it up testing by generating a sequence of square roots of non squares.

<<newton_sqrt.a68>>=PROC newton sqrt = (LONG REAL n)LONG REAL: ( Newton's method; stop iteration: xn ); main: ( FOR n TO 27 DO LONG REAL value = n + ENTIER (newton sqrt(n) + 0.5); printf(($"Square root of "g(-2)": "gl$, value, newton sqrt(value))) OD )

Output:

Square root of 2: +1.414213562373095048801688724e +0 Square root of 3: +1.732050807568877293527446342e +0 Square root of 5: +2.236067977499789696409173669e +0 Square root of 6: +2.449489742783178098197284075e +0 Square root of 7: +2.645751311064590590501615754e +0 Square root of 8: +2.828427124746190097603377448e +0 Square root of 10: +3.162277660168379331998893545e +0 Square root of 11: +3.316624790355399849114932737e +0 Square root of 12: +3.464101615137754587054892683e +0 Square root of 13: +3.605551275463989293119221268e +0 Square root of 14: +3.741657386773941385583748732e +0 Square root of 15: +3.872983346207416885179265400e +0 Square root of 17: +4.123105625617660549821409856e +0 Square root of 18: +4.242640687119285146405066173e +0 Square root of 19: +4.358898943540673552236981984e +0 Square root of 20: +4.472135954999579392818347338e +0 Square root of 21: +4.582575694955840006588047194e +0 Square root of 22: +4.690415759823429554565630114e +0 Square root of 23: +4.795831523312719541597438064e +0 Square root of 24: +4.898979485566356196394568149e +0 Square root of 26: +5.099019513592784830028224109e +0 Square root of 27: +5.196152422706631880582339025e +0 Square root of 28: +5.291502622129181181003231507e +0 Square root of 29: +5.385164807134504031250710492e +0 Square root of 30: +5.477225575051661134569697828e +0 Square root of 31: +5.567764362830021922119471299e +0 Square root of 32: +5.656854249492380195206754897e +0

As mentioned, the algorithm can be generalized to the *n*th
root of *N* by using *f*(*x*) = *x*^{n} − *N* and
*f*'(*x*) = *n**x*^{n − 1}.

Download code |