Secant Method

Description

The secant method is an iterative procedure to estimate a root of an equation f(x) = 0 where the user gives two initial estimates and a tolerance. One of the two initial estimates is regarded as being older than the other, the procedure then calculates a new estimate and discards the oldest estimate until the distance between the two estimates is less than the preassigned tolerance. If x_i-1 and x_i are the two estimates with x_i-1 being older than x_i, the new estimate x_i+1 is that point on the x-axis which is the intersection of the x-axis and the line joining the points ( x_i-1, f(x_i-1) ) and ( x_i, f(x_i) ) on the curve y = f(x). The subsequent estimate then uses the points ( x_i, f(x_i) ) and ( x_i+1, f(x_i+1) ) to estimate x_i+1.

It is best to avoid using the secant method if there are local extrema near the root. The convergence of the secant method depends on the initial estimates. It is generally better if the two initial conditions are close so that the secant is approximately equal to the tangent.

Function List

Functions:

double Secant_Method( double (*f)(double), double a, double b, double tolerance, int max_iteration_count, int *err)

Find a root of f(x) in near a and b. 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 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 exceed 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: secant_method.c

////////////////////////////////////////////////////////////////////////////////
// File: secant_method.c                                                      //
// Routines:                                                                  //
//    double Secant_Method                                                    //
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
//  double Secant_Method( double (*f)(double), double a, double b,            //
//                      double tolerance, int max_iteration_count, int *err)  //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the secant method where 'a'     //
//     and 'b' are two initial estimates.  The process terminates successfully//
//     when the change in the estimate of the root is less than 'tolerance'   //
//     or unsuccessfully when the number of iterations reaches 'max_iteration_//
//     count' or if an attempt is made to divide by zero.                     //
//                                                                            //
//     The secant method uses two points on the curve f(x) to estimate the    //
//     root as the intersection of the x-axis with the line joining the       //
//     two points (the secant).  The process then replaces the oldest of the  //
//     two estimates with the newest estimate.  The process continues until   //
//     the distance between the two estimates is less than the tolerance.     //
//     The process can fail if there is a local maximum or minimum near the   //
//     the root and is sensitive to having good initial estimates.  If there  //
//     is a local maximum or minimum in the neigborhood of the root, the      //
//     return value should be verified before use.  The secant method is      //
//     similar to the Newton-Raphson method in which the derivative is        //
//     replaced by a numerical estimate of the derviative.                    //
//                                                                            //
//  Arguments:                                                                //
//     double *f         Pointer to function of a single variable of type     //
//                       double.                                              //
//     double  a         Initial estimate.                                    //
//     double  b         Initial estimate.                                    //
//     double  tolerance Desired accuracy of the zero, tolerance > 0.         //
//     int     max_iteration_count Maximum allowable number of iterations.    //
//     int     *err      0 if successful, -1 if not, i.e. the iteration       //
//                       count meets or exceeds the max_iteration_count,      //
//                       -2 if an attempt is made to divide by zero.          //
//                                                                            //
//  Return Values:                                                            //
//     A zero somewhere in a neighborhood of a and b.                         //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), zero, a, b, tolerance = 1.e-6;                    //
//        int max_iteration = 10;                                             //
//        int err;                                                            //
//                                                                            //
//        (determine nearby values a and b of a zero)                         //
//        zero = Secant_Method( f, a, b, tolerance, max_iterations, &err);    //
//        ...                                                                 //
//     }                                                                      //
//     double f(double x) { define f }                                        //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Secant_Method( double (*f)(double), double a, double b, 
                        double tolerance, int max_iteration_count, int *err) {

   double fa = (*f)(a), fb = (*f)(b), delta;
   int i;

   *err = 0;
   for (i = 0; i < max_iteration_count; i++) {

         // Insure that the f(b) has smaller magnitude than f(a). //

      if ( fabs(fa) < fabs(fb) ) {
         delta = a; a = b; b = delta; delta = fa; fa = fb; fb = delta;
      }
    
         // If an attempt to divide by zero, return with an error code -2. //

      if ( fb == fa ) { *err = -2; return 0.0; }

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

      delta = fb * (a - b) / (fb - fa);

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

      if ( delta == 0.0 ) return b;

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

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

      a = b + delta;
      fa = (*f)(a);
   }
   
   *err = -1;
   return a;   
}