////////////////////////////////////////////////////////////////////////////////
// File: illinois_algorithm.c //
// Routines: //
// double Illinois_Algorithm //
////////////////////////////////////////////////////////////////////////////////
#include < math.h > // required for fabs()
////////////////////////////////////////////////////////////////////////////////
// double Illinois_Algorithm( double (*f)(double), double a, double b, //
// double tolerance, int max_iterations, int *err) //
// //
// Description: //
// Estimate the root (zero) of f(x) using the Illinois algorithm 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 //
// length < '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. //
// //
// The Illinois algorithm is a variant of the regula falsi method in //
// which during the iterative process if a change to the same endpoint //
// occurs twice in succession, then the ordinate corresponding to the //
// other endpoint is halved. The Illinois algorithm always converges. //
// //
// 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. //
// 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; //
// //
// (determine lower bound, a and upper bound, b of a zero) //
// zero = Illinois_Algorithm( f, a, b, tolerance, &err); //
// } //
// double f(double x) { define f } //
////////////////////////////////////////////////////////////////////////////////
// //
double Illinois_Algorithm( 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 ab = 1, 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 ab = 1, then the endpoint a 'should have' changed, //
// and if ab = 0, then the endpoint b 'should have changed. //
// If the endpoint that should have changed didn't, then //
// halve its ordinate value. //
if ( ab ) ( fc > 0.0 ) ? ( ab = 0 ) : ( fa *= 0.5 );
else ( fc > 0.0 ) ? ( fb *= 0.5 ) : ( ab = 1 );
// 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 a - fa * (b - a) / (fb - fa);
}