//////////////////////////////////////////////////////////////////////////////// // File: bisection_method.c // // Routines: // // int Bisection_Method // //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// // int Bisection_Method( double (*f)(double), double *a, double *fa, // // double *b, double *fb, void **data, // // int (*verify_setup)(double, double, double, double, void**, int*), // // int (*terminate)(double, double, double, double, void**, int*), *int)) // // // // Description: // // Estimate the root (zero) of f(x) using the bisection method where // // 'a' and 'b' are initial estimates where either f(a) > 0 and f(b) < 0 // // or f(a) < 0 and f(b) > 0. // // // // 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 // // iterative process is continued until the terminated by the user or // // either there is no machine representable number between the two // // endpoints or the function evaluates to zero at the new midpoint. // // // // Arguments: // // double *f // // Pointer to a user supplied function of a single variable of type // // double. // // double *a // // The address of one of the bounds of an interval which contains a // // zero of f(x). On input, *a is the value of either the lower or // // upper bound of the initial interval. On output, *a is the value // // of the final bound of that endpoint of the reduced interval. // // f(*a) * f(*b) < 0. // // double *fa // // The address of the value of f(x) evaluated at *a. On input, *fa // // is the value f(*a) where *a is the initial value of the bound and // // on output *fa will be the value of f(*a) evaluated at the final // // value *a. // // double *b // // The address of the other bound of the interval which contains a // // zero of f(x). On input, *b is the value of that bound of the // // initial interval. On output, *b is the value of the final bound // // of that endpoint of the shrunken interval. f(*a) * f(*b) < 0. // // double *fb // // The address of the value of f(x) evaluated at *b. On input, *fb // // is the value f(*b) where *b is the initial value of the bound and // // on output *fb will be the value of f(*b) evaluated at the final // // value *b. // // void *data[] // // This is an array of pointers containing data the user wants to // // pass to the user's terminate procedure and consequently to the // // user's verify procedure. This array is not used directly in // // this routine, but only passed through to the user's routines. // // int (*verify_setup)() // // The user supplied function which checks that the initial input // // data, a, fa, b, fb, and the array data, comply with the // // requirements that the user supplied function terminate() needs // // so that it functions correctly. Since the array data is an // // array of pointers, one can only check that the pointers are // // non-zero and if they are non-zero then the objects that they // // point to meet the requirements set forth by the terminate() // // procedure. // // The return values of verify_setup() are: // // 0 - The setup is ok. // // Otherwise - There is a problem with the setup. // // If there is a problem with the setup, this Bisection_Method() // // terminates immediately with a return value SETUP_ERROR, 1. // // The function verify_setup() is declared as follows: // // int verify_setup(double a, double fa, double b, double fb, // // void **data, int *msg) // // and is called as follows: // // (*verify_setup)(*a, *fa, *b, *fb, data, msg); // // where *a, *fa, *b, *fb, are the dereferenced argments to the // // Bisection_Method and data and msg are arguments to // // Bisection_Method. // // An example of a verify_setup() procedure is // // Bisection_Method_Verify_Setup // // which is associated with the example terminate procedure // // Bisection_Method_Terminate. // // int (*terminate)() // // The user supplied function to control whether of not to terminate // // the iterative process. The function, terminate(), returns a 0 if // // the process is to continue or a 1 if the process is to terminate. // // The function terminate() is declared as follows: // // int terminate(double a, double fa, double b, double fb, // // void **data, int *msg) // // and is called as follows: // // (*terminate)(*a, *fa, *b, *fb, data, msg) // // where *a, *fa, *b, *fb, are the dereferenced argments to the // // Bisection_Method and data and msg are arguments to // // Bisection_Method. // // An example of a terminate procedure is // // Bisection_Method_Terminate. // // int *msg // // An additional parameter passed to the user routines verify_setup // // and terminate to allow the user to pass back an error code // // in the case of verify_setup or a halt code in the case of // // terminate. // // // // Return Values: // // The return values are: // // 4 (FOUND_ZERO) if f(x) evaluates to 0 at (*a + *b) / 2. // // 3 (NO_POSSIBLE_NUMBER) if there is no machine representable number // // between *a and *b. // // 2 (TERMINATED_BY_USER) if the subprogram is terminated by the // // user's function 'terminate().' // // 1 (SETUP_ERROR) if the user's functions 'verify_setup()' // // returned a 0 (ERR). // // 0 (ILLEGAL_POINTER) at least one of the argument pointers is // // NULL. // // // // Example: // // // // extern double f(double); // // extern int Bisection_Method_Verify_Setup( double, double, double, // // double, void **, int *); // // extern int Bisection_Method_Terminate( double, double, double, // // double, void **, int *); // // double a, b, fa, fb; // // int msg, k; // // void *data[4]; // // // // (* determine a, b, fa, fb such that fa * fb < 0.0 *) // // (* setup the array data - abs and rel tolerance and test region *) // // // // k = Bisection_Method( f, &a, &fa, &b, &fb, data, // // Bisection_Method_Verify_Setup, Bisection_Method_Terminate, &msg);// //////////////////////////////////////////////////////////////////////////////// // // enum {ILLEGAL_POINTER, SETUP_ERROR, TERMINATED_BY_USER, NO_POSSIBLE_NUMBER, FOUND_ZERO}; int Bisection_Method( double (*f)(double), double *a, double *fa, double *b, double *fb, void *data[], int (*verify_setup)(double, double, double, double, void**, int*), int (*terminate)(double, double, double, double, void**,int*), int *msg) { double c, fc; double *swap; if ( (f == 0) || (a == 0) || (fa == 0) || (b == 0) || (fb == 0) || (data == 0) || (verify_setup == 0) || (terminate == 0) || (msg == 0) ) return ILLEGAL_POINTER; if (!verify_setup( *a, *fa, *b, *fb, data, msg) ) return SETUP_ERROR; // Insure that a points to the endpoint where the function is positive. // if ( *fa < 0.0 ) { swap = fa; fa = fb; fb = swap; swap = a; a = b; b = swap; } // Continue halving the interval [*a, *b] if *a < *b or [*b, *a] if // // *b < *a until terminated by the user's supplied function terminate // // or either there is no machine representable number between *a and *b// // or the function evaluates to zero at the midpoint. // while ( !terminate(*a, *fa, *b, *fb, data, msg) ) { c = 0.5 * *a + 0.5 * *b; if ( (c == *a) || (c == *b) ) { return NO_POSSIBLE_NUMBER; } fc = (*f)(c); if (fc == 0.0) { return FOUND_ZERO; } if ( fc > 0.0 ) { *a = c; *fa = fc; } else {*b = c; *fb = fc; } } return TERMINATED_BY_USER; }