Bisection Method

Description

The bisection method is one of the simplest to understand and program and one of the slowest. In order to apply the bisection method, two approximations to a root must be obtained which bracket the root. Moreover, the function evaluated at the two approximations must be of opposite signs, i.e. if a and b are two numbers such that a < b, then f(a)·f(b) < 0. Note that this technique cannot be applied to functions such as f(x) = x².

Function List

Functions:

double Bisection_Method( double (*f)(double), double a, double b, double tolerance, int *err)

Find a root of f(x) in the interval ( min(a,b), max(a,b) ) where f(a)·f(b) < 0. 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 method returns an err of 0 if the iteration was successful and -1 if the initial estimates fail to satify the requirement that the function evaluated at the two values must have opposite signs.

C Source

File: bisection_method.c

////////////////////////////////////////////////////////////////////////////////
// File: bisection_method.c                                                   //
// Routines:                                                                  //
//    double Bisection_Method                                                 //
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
//  double Bisection_Method( double (*f)(double), double a, double b,         //
//                                          double tolerance, int *err)       //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the bisection method where      //
//     'a' and 'b' are initial estimates which bracket the root i.e. either   //
//     f(a) > 0 and f(b) < 0 or f(a) < 0 and f(b) > 0.  The iterative process //
//     terminates when the root is contained within an interval of length <   //
//     'tolerance', in which case the value returned is the midpoint of that  //
//     interval, or if during the iterative process the function vanishes at  //
//     a test point, the test point is returned.                              //
//                                                                            //
//     The bisection method evaluates the function at the midpoint of the     //
//     interval, then the endpoint of the interval where evaluation of the    //
//     function has the same sign as the function evaluated at the midpoint   //
//     is replaced with the midpoint, thus halving the interval.  The         //
//     process is continued until the length of the interval is less than     //
//     the preassigned tolerance.  The bisection method always converges,     //
//     but is slow.                                                           //
//                                                                            //
//  Arguments:                                                                //
//     double *f         Pointer to function of a single variable of type     //
//                       double.                                              //
//     double  a         Initial estimate, f(a)*f(b) < 0.                     //
//     double  b         Initial estimate. f(a)*f(b) < 0.                     //
//     double  tolerance Desired accuracy of the zero, tolerance > 0.         //
//     int     *err      0 if successful, -1 if not, i.e. if f(a)*f(b) > 0.   //
//                                                                            //
//  Return Values:                                                            //
//     The root contained within the interval (a,b) if a < b or (b,a) if not. //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), zero, a, b, tolerance = 1.e-6;                    //
//        int err;                                                            //
//                                                                            //
//        (determine lower bound, a and upper bound, b of a zero)             //
//        zero = Bisection_Method( f, a, b, tolerance, &err);                 //
//        ...                                                                 //
//     }                                                                      //
//     double f(double x) { define f }                                        //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Bisection_Method( double (*f)(double), double a, double b, 
                                                 double tolerance, int *err) {
   double fa = (*f)(a), fc = fa * (*f)(b), c;

      // If the initial estimates do not bracket a root, set the err flag and //
      // return.  If an initial estimate is a root, then return the root.     //

   if ( fc > 0.0 ) { *err = -1; return 0.0; }
   *err = 0;
   if ( fc  == 0.0 ) return ( fa == 0.0 ? a : b );

      // While the interval length is greater than the preassigned tolerance,//
      // continue halving the interval.  When the interval length becomes    //
      // less than the tolerance, return the midpoint of that interval.      //

   if ( fa > 0.0 ) 
      while ( fabs( b - a ) > tolerance )
         if ( (fc = (*f)( c = 0.5 * (a + b) ) ) == 0.0 ) return c; 
         else ( fc > 0.0 ) ?  (a = c)  : (b = c);
   else 
      while ( fabs( b - a ) > tolerance ) 
         if ( (fc = (*f)( c = 0.5 * (a + b) ) ) == 0.0 ) return c; 
         else ( fc > 0.0 ) ?  (b = c)  : (a = c);

      // Return the midpoint of the final interval. //

   return 0.5 * ( a + b );
}