Newton-Raphson Method, Part 1

Newton-Raphson Method, Part 1

We give a procedure for solving a nonlinear equation in one unknown x, namely, we want to solve the equation f(x) = 0 for the unknown x, where f(x) is nonlinear in x.

For example, we want to solve the equation f(x) = x^2 - 4*sin(x) = 0 for some value of x > 0.

The Newton-Raphson (NR) method is an iterative procedure that is started with an initial guess x_0, and successively computes new guesses x_1, x_2, x_3, ... , x_k, ... , using the Newton-Raphson formula

x_{k+1} = x_k - f(x_k) / f'(x_k)

Under appropriate conditions on f(), in an interval [a,b] the sequence {x_k} will converge to a solution in that interval.

The problem with the NR method is that for each new function f(), we need to write a new program to compute the three entities f(x), f'(x), and {x_k}. Here we give a procedure that uses a fixed main module, main.c, that reads the initial value x_0 and uses externally defined functions f(x) and f'(x) in an external module funderiv.c, to compute the series {x_k} and provide reasonable stopping conditions. We also provide a glue script, nrsolve.sh, that compiles main.c and funderiv.c to give the executable nrsolve, and then runs nrsolve, in order to produce the sequence {x_k}.

We can pre-debug both the script nrsolve.sh and the main module main.c, and make sure that these are both working properly.  Thus the user only needs to write and debug the module funderiv.c that contains f(x) and f'(x), to produce a working NR-solution for that function f(x).

Here is our script nrsolve.sh

*** nrsolve.sh ***
#!/bin/bash

if [ $# -ne 2 ]; then
  echo "Run program as:\n  nrsolve.sh funderiv.c x_0"
  exit 1
fi

gcc main.c $1 -lm -o nrsolve

if [ $? ]; then
  ./nrsolve $2
fi
*** end ***

This script produces the executable nrsolve from the two modules main.c and funderiv.c (which may actually have another name), and then if compilation is successful, runs the resulting nrsolve executable, using the provided starting value x_0.

And here is our main.c module

*** main.c ***
#include <stdio.h>
#include <stdlib.h>
#define EPS 1.0e-12L
#define MAXITER 1000
#define ABS(x) ((x) < 0.0? -(x):(x))

extern long double fun(long double x);
extern long double deriv(long double x);
long double sol;
int j;

void newtonraphson(void)
{
  long double new;
  long double fsol;
  for( j = 0; j < MAXITER ; ++j, sol = new) {
    fsol = fun(sol);
    fprintf(stderr, "x=%.14Lf f(x)=%.14Lf\n", sol, fsol);
    if(ABS(fsol) < EPS) break;
    new = sol - fsol / deriv(sol);
  }
}

int main(int argc, char *argv[])
{
  long double a;
  if( argc != 2) exit(1);
  sscanf(argv[1], "%Lf", &sol);
  fprintf(stderr, "Starting value supplied: %Lf\n", sol);
  newtonraphson();
  fprintf(stderr, "Solution: %.14Lf, error %.14Lf after %d iters\n", sol, EPS, j);
  exit(0);
}
*** end ***

To complete our NR solution for the function f(x) = x^2 - 4*sin(x) = 0, see Newton-Raphson Method, Part 2