Newton-Raphson Method

Description

The Newton-Raphson method is an iterative procedure to estimate a root of an equation f(x) = 0 where the user gives an initial estimate, a tolerance and programs the derivative of the function f(x). The procedure then calculates a new estimate until the distance between the two estimates is less than the preassigned tolerance. If x_i-1 is the current estimate, the new estimate x_i is that point on the x-axis which is the intersection of the x-axis and the tangent to the curve y = f(x) at the point ( x_i-1, f(x_i-1) ). The subsequent estimate then uses the point ( x_i, f(x_i) ).

It is best to avoid using the Newton-Raphson method if there are local extrema near the root and if the root has order > 1. The convergence of the Newton-Raphson method depends on the initial estimate and on the order of the root. For roots with order 1, convergence is quadratic near the root, however Newton-Raphson requires calculating not only the function but also its derivative at each step. For this reason, Newton-Raphson may require fewer steps, but more function calls than other methods.

Function List

Functions:

double Newton_Raphson( double (*f)(double), double (*df)(double), double a, double tolerance, int max_iteration_count, int *err)

Find a root of f(x) near a. The absolute difference of the returned value and the actual root will be less than tolerance, a user specified number specifying the desired accuracy of the result. The input argument (*df)(double) is the value of the derivative of f(x) evaluated at the argument. The input argument max_iteration_count allows the user to control the maximum number of iterations to attempt. The method returns an err of 0 if the iteration was successful, -1 if the number of iterations exceeded the maximum allowable number of iterations as specified by the user, and -2 if an attempt was made to divide by zero, which is indicative of a local minimum or a local maximum. Even if the an err return of 0 occurs, it is possible that the result is erroneous, especially if there is a local minimum or local maximum near the root. For this reason, if it is not known whether there is a local extremum in a small neighborhood of the root, the result should be checked for its validity.

C Source

File: newton_raphson.c

////////////////////////////////////////////////////////////////////////////////
// File: newton_raphson.c                                                     //
// Routines:                                                                  //
//    double Newton_Raphson                                                   //
////////////////////////////////////////////////////////////////////////////////

#include < math.h >                                      // required for fabs()

////////////////////////////////////////////////////////////////////////////////
//  double Newton_Raphson( double (*f)(double), double (*df)(double),         //
//            double a, double tolerance, int max_iteration_count, int *err)  //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the Newton-Raphson method       //
//     with 'a' as an initial estimate.  The process terminates successfully  //
//     when either the function vanishes or when the change in the estimate   //
//     of the root is less than 'tolerance' or unsuccessfully when the        //
//     number of iterations reaches 'max_iteration_count'.                    //
//                                                                            //
//     The Newton-Raphson method uses the intersection of the x-axis and the  //
//     tangent to the curve at an estimate of the root to obtain a new        //
//     estimate of the location of the root.  The process is continued until  //
//     the change of locations of the estimates is within a preassigned       //
//     tolerance.  The Newton-Raphson method requires the user not only       //
//     furnish the function but also its derivative.  The method can fail if  //
//     the root is near a local minimum or a local maximum.  If the order     //
//     of the zero of f(x) is greater than one, one can find a root of        //
//     f(x) / f'(x).  The Newton-Raphson algorithm is quadratically           //
//     convergent near a simple root.                                         //
//                                                                            //
//  Arguments:                                                                //
//     double *f         Pointer to function of a single variable of type     //
//                       double.                                              //
//     double *df        Pointer to the derivative of the function *f.        //
//     double  a         Initial estimate.                                    //
//     double  tolerance Desired accuracy of the zero.                        //
//     int     max_iteration_count Terminate iteration if the number of       //
//                                 iterations meets or exceed this number.    //
//     int     *err      0 if successful, -1 if not, i.e. the number of       //
//                       iterations meets or exceeds the max_iteration_count, //
//                       or -2 if the derivative vanishes at an evaluation    //
//                       point.                                               //
//                                                                            //
//  Return Values:                                                            //
//     A root somewhere in a neighborhood of a.                               //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), df(double),zero, a, tolerance = 1.e-6;            //
//        int max_iteration = 10;                                             //
//        int err;                                                            //
//                                                                            //
//        (determine the initial estimate a of the root.)                     //
//        zero = Newton_Raphson( f, df, a, tolerance, max_iteration, &err);   //
//        ...                                                                 //
//     }                                                                      //
//     double f(double x) { define f }                                        //
//     double df(double x) { define df }                                      //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Newton_Raphson( double (*f)(double), double (*df)(double), double a,
                       double tolerance, int max_iteration_count, int *err) {

   double delta, fa, dfa;
   int i;
   
   *err = 0;
   for (i = 0; i < max_iteration_count; i++) {
      if ( ( fa = (*f)(a) ) == 0.0 ) return a;

         // If an attempt to divide by zero, return with an error code -2. //

      if ( ( dfa = (*df)(a) ) == 0.0 ) { *err = -2; return 0.0; }

         // Estimate the location of the root relative to a. //

      delta = - fa / dfa;

         // If the location of the root relative to a is less than //
         // the tolerance, return the estimate of the root.        //

      if (fabs(delta) < tolerance ) return a + delta;

         // Otherwise, update the estimate of the root. //

      a += delta;
   }

   *err = -1;
   return a;
}