Newton-Raphson's method for root finding (C)

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
while (!(guess sufficient))
{
  refine guess
}

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>>=
double x = 0.0;
double xn = 0.0;

How long do we iterate? Let's say 100 steps are sufficient for our purposes. We could iterate shorter, getting a less accurate result. We could iterate longer, getting a more accurate one should we really need it.

<<guess sufficient>>=
iters++ >= 100

iters is another variable we must define.

<<define local variables>>=
int iters = 0;

However, we don't have to do all these iterations. 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. In practice, you'll notice that the algorithm will often terminate on this condition long before the 100 iterations are reached.

<<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>>=
int i;
for (i = 0; i <= (int)n; ++i)
{
  double val = i*i-n;
  if (val == 0.0)
    return i;
  if (val > 0.0)
  {
    xn = (i+(i-1))/2.0;
    break;
  }
}

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 with some small test code.

<<newton_sqrt.c>>=
#include <stdio.h>
double newton_sqrt(double n)
{
  Newton's method
  return xn;
}
int main(void)
{
  printf("Square root of 5: %f\n", newton_sqrt(5));
  printf("Square root of 25: %f\n", newton_sqrt(25));
  printf("Square root of 20: %f\n", newton_sqrt(20));
  return 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