Regula Falsi Method

Description

The regula falsi method is a variant of the secant method in which the initial estimates bracket the root in such a manner that the function evaluated at the two initial estimates has opposite signs. The regula falsi iterative procedure uses the two estimates as in the secant method to estimate a new estimate, but the estimate which is replaced is the one for which the sign of the function is the same as the sign of the function for the new estimate. The regula falsi method is guaranteed to converge, but convergence may be extremely slow.

It is best to avoid using the regula falsi method for final convergence to a root. It can be used as a good starting method. The restriction that the two estimates bracket the root means that regula falsi avoids the problems associated with having local extrema since it is no longer possible to attempt to divide by zero. Since regula falsi requires initial estimates which bracket the root in which the function evaluated at the initial estimates has opposite signs, it is impossible to use regula falsi to find the zeros of strictly nonnegative functions or nonpositive functions.

Function List

Functions:

double Regula_Falsi( double (*f)(double), double a, double b, double tolerance, int max_iterations, 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 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 initial estimates fail to satify the requirement that the function evaluated at the two values must have opposite signs, and -2 if the number of iterations exceed the maximum allowable number of iterations as specified by the user.

C Source

File: regula_falsi.c

////////////////////////////////////////////////////////////////////////////////
// File: regula_falsi.c                                                       //
// Routines:                                                                  //
//    double Regula_Falsi                                                     //
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
//  double Regula_Falsi( double (*f)(double), double a, double b,             //
//                           double tolerance, int max_iterations, int *err)  //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the regula falsi 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       //
//     iteration terminates when the zero is constrained to be within an      //
//     interval of width < 'tolerance', in which case the value returned is   //
//     the point on the x-axis corresponding to the intersection of the       //
//     x-axis and the line joining the last pair of points, or when the       //
//     number of iterations meets or exceeds 'max_iterations'.                //
//                                                                            //
//     The regula falsi method is a variant of the secant method in which     //
//     the initial estimates must bracket a root.  Each successive estimate   //
//     replaces that endpoint so that the resulting interval still brackets   //
//     a root.  This method is slow and is considered a good starting method  //
//     but not a method to use for the final convergence to the root.         //
//                                                                            //
//  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     max_iterations  The maximum allowable number of iterations.    //
//     int     *err      0 if successful, -1 if not, i.e. if f(a)*f(b) > 0,   //
//                       -2 if the number of iterations >= max_iterations.    //
//                                                                            //
//  Return Values:                                                            //
//     A root contained within the interval (a,b), if a < b, or (b, a),       //
//     if a > b.                                                              //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), zero, a, b, tolerance = 1.e-6;                    //
//        int err, max_iterations = 100;                                      //
//                                                                            //
//        (determine lower bound, a, and upper bound, b, of a zero)           //
//        zero = Regula_Falsi( f, a, b, tolerance, max_iterations, &err);     //
//     }                                                                      //
//     double f(double x) { define f }                                        //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Regula_Falsi( double (*f)(double), double a, double b, 
                             double tolerance, int max_iterations, int *err) {

   double fa = (*f)(a), fb = (*f)(b), c, fc = fa * fb;
   int i;
 
      // 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.     //
   
   *err = 0;
   if ( fc >= 0.0 )
     if ( fc > 0.0 ) { *err = -1; return 0.0; }
     else return ( fa == 0.0 ? a : b );

      // Insure that the f(a) is positive. //

   if ( fa < 0.0 ) { c = a; a = b; b = c; c = fa; fa = fb; fb = c; }

      // While the length of the interval containing a root is greater than  //
      // the preassigned tolerance, estimate a new endpoint by the           //
      // intersection of the line through points on the curve at the         //
      // endpoints (the secant) and the x-axis.  The endpoint which is       //
      // replaced is the one for which the new estimate still brackets the   //
      // root.                                                               //

   for ( i = 0; i < max_iterations; i++ ) {
         
         // Estimate the location of the root. //

      fc = ( fa > -fb ) ? (*f)( c = b - fb * (b - a) / (fb - fa) ) :
                          (*f)( c = a - fa * (b - a) / (fb - fa) );
       
         // If the new estimate of a root is a root, return. //
 
      if ( fc == 0.0 ) return c;

         // If not, update the interval endpoints. //
 
      ( fc > 0.0 ) ?  ( a = c, fa = fc ) : ( b = c, fb = fc );

         // If the length of the interval containing a root is less //
         // than the preassigned tolerance, return the point where  //
         // secant crosses the x-axis.                              //

      if ( fabs( a - b ) < tolerance ) 
         return ( fa > -fb ) ? b - fb * (b - a) / (fb - fa) :
                               a - fa * (b - a) / (fb - fa) ;
   }

   *err = -2;
   return c;
}