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² − N. The x where f(x) = 0 is the square root of N. Newton's method starts with a guess x₀ 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₀ = 1.5. The red line is the tangent at f(1.5). Where it crosses the x axis is our refined approximation x₁. With each iteration, we make a step towards the true zero of f(x) = x² − 2, marked as r.

Given x_i, we calculate x_{i + 1} by . For f(x) = x² − N we have f'(x) = 2x 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 nth root of N by using f(x) = xⁿ − N and f'(x) = nx^{n − 1}.

Download code