//////////////////////////////////////////////////////////////////////////////// // File: regula_falsi.c // // Routines: // // int Regula_Falsi // //////////////////////////////////////////////////////////////////////////////// #include // required for fabs() //////////////////////////////////////////////////////////////////////////////// // int Regula_Falsi( 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, unsigned int, // // void**,int*), // // int *msg) // // // // Description: // // Estimate the root (zero) of f(x) using the regula falsi 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 regula falsi method is a variant of the secant method in which // // the initial estimates must bracket a root such that the function // // evaluated at the initial estimates have opposite signs. Given initial // // estimates of a root by a and b for which f(a) f(b) < 0, a new esimate // // is given by finding the intersection c of the line joining (a,f(a)) and// // (b, f(b)) with the x-axis. The successive interval which brackets // // a root is obtained by relacing that endpoint for which the function // // evaluated at that endpoint has the same sign as that of the function // // evaluated at c with c. // // // // 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 endpoint. // // // // 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 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 Regula_Falsi() // // terminates immediately with a return value SETUP_ERROR, 0. // // 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 // // Regula_Falsi and data and msg are arguments to // // Regula_Falsi. // // An example of a verify_setup() procedure is // // Regula_Falsi_Verify_Setup // // which is associated with the example terminate procedure // // Regula_Falsi_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, // // unsigned int iteration_count, void **data, int *msg) // // and is called as follows: // // (*terminate)(*a, *fa, *b, *fb, iteration_count, data, msg) // // where *a, *fa, *b, *fb, are the dereferenced argments to the // // Regula_Falsi and data and msg are arguments to // // Regula_Falsi and iteration_count is an internal variable to // // Regula_Falsi. // // An example of a terminate procedure is // // Regula_Falsi_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 the intersection // // of the secant with the x-axis. // // 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. // // // // The return value is either 3 (ZERO) if f(x) evaluates to 0 at // // the intersection of the secant with the x-axis, 2 (NO_POSSIBLE_NUMBER) // // if there is no machine representable number between *a and *b, // // 1 (TERMINATED_BY_USER) if the subprogram is terminated by the user's // // function 'terminate', or 0 (SETUP_ERROR) if there is an error with // // the setup data. // // // // Example: // // // // extern double f(double); // // extern int Regula_Falsi_Verify_Setup( double, double, double, // // double, void **, int *); // // extern int Regula_Falsi_Terminate( double, double, double, // // unsigned int, double, void **, int *); // // double a, b, fa, fb; // // int msg, k; // // unsigned int max_number_of_iterations; // // void *data[5]; // // // // (* determine a, b, fa, fb such that fa * fb < 0.0 *) // // (* setup the array data - abs and rel tolerance and test region *) // // // // k = Regula_Falsi( f, &a, &fa, &b, &fb, data, // // Regula_Falsi_Verify_Setup, Regula_Falsi_Terminate, &msg);// //////////////////////////////////////////////////////////////////////////////// // // enum {ILLEGAL_POINTER, SETUP_ERROR, TERMINATED_BY_USER, NO_POSSIBLE_NUMBER, FOUND_ZERO}; int Regula_Falsi( 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, unsigned int, void**,int*), int *msg) { unsigned int iteration_count = 0; double *swap; double r,s; double c, fc; 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 the f(*a) is positive. // if ( *fa < 0.0 ) { swap = fa; fa = fb; fb = swap; swap = a; a = b; b = swap; } // Continue finding the intersection c of the line joining the two // // points (*a, f(*a)) and (*b, f(*b)) with the x-axis, and replacing // // the endpoint for which the value of the function has the same sign // // as the function evaluated at c until terminated by the user's // // supplied function terminates the iteration or either there is no // // machine representable number between *a and *b or the function // // evaluates to zero at c. // while ( !terminate(*a, *fa, *b, *fb, iteration_count++, data, msg) ) { // Estimate the location of the root. // if( *fa > - *fb ) { s = *fb / *fa; r = 1.0 - s; c = *b / r - (s * *a) / r; } else { s = *fa / *fb; r = 1.0 - s; c = *a / r - (s * *b) / r; } // Check if the new estimate of a root corresponds to an existing // // endpoint and if it is, check if the halving the interval also // // yields an existing endpoint, it so then return NO_POSSIBLE_NUMBER // // otherwise use the midpoint as the new estimate. // if ( (c == *a) || (c == *b) ) { c = 0.5 * *a + 0.5 * *b; if ( (c == *a) || (c == *b) ) { return NO_POSSIBLE_NUMBER; } } // If the new estimate of a root is a root, return. // fc = f(c); if (fc == 0.0) {*a = c; *b = c; *fa = fc; *fb = fc; return FOUND_ZERO;} // If not, update the interval endpoints. // ( fc > 0.0 ) ? ( *a = c, *fa = fc ) : ( *b = c, *fb = fc ); } return TERMINATED_BY_USER; }