//////////////////////////////////////////////////////////////////////////////// // File: illinois_algorithm.c // // Routines: // // int Illinois_Algorithm // //////////////////////////////////////////////////////////////////////////////// #include // required for fabs() //////////////////////////////////////////////////////////////////////////////// // int Illinois_Algorithm( 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: // // // // The Illinois algorithm is a variant of the regula falsi 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. If the right or left endpoint is successively // // replaced the ordinate of the other endpoint is halved. Initially, the // // ordinate is the value of the function, f(x), at the endpoint. // // // // 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 always converges. // // // // 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 Illinois_Algorithm() // // 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 // // Illinois_Algorithm and data and msg are arguments to // // Illinois_Algorithm. // // An example of a verify_setup() procedure is // // Illinois_Algorithm_Verify_Setup // // which is associated with the example terminate procedure // // Illinois_Algorithm_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 // // Illinois_Algorithm and data and msg are arguments to // // Illinois_Algorithm and iteration_count is an internal variable// // to Illinois_Algorithm. // // An example of a terminate procedure is // // Illinois_Algorithm_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. // // // // Example: // // // // extern double f(double); // // extern int Illinois_Algorithm_Verify_Setup( double, double, double, // // double, void **, int *); // // extern int Illinois_Algorithm_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 = Illinois_Algorithm( f, &a, &fa, &b, &fb, data, // // Illinois_Algorithm_Verify_Setup, Illinois_Algorithm_Terminate,// // &msg);// //////////////////////////////////////////////////////////////////////////////// // // enum {ILLEGAL_POINTER, SETUP_ERROR, TERMINATED_BY_USER, NO_POSSIBLE_NUMBER, FOUND_ZERO}; enum {REPLACED_A, REPLACED_B, STARTING}; int Illinois_Algorithm( 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; int ab = STARTING; 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 and if this endpoint (right or left)// // is successively replaced then halve the ordinate of the other // // endpoint 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) ) { if (ab == REPLACED_A) *fb = f(*b); if (ab == REPLACED_B) *fa = f(*a); 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 ); // If the same endpoint (right or left) has been // // successively replaced then halve the ordinate // // of the other endpoint. // switch (ab) { case REPLACED_A: if (fc > 0.0) *fb *= 0.5; else ab = REPLACED_B; break; case REPLACED_B: if (fc < 0.0) *fa *= 0.5; else ab = REPLACED_A; break; default: ab = (fc > 0.0) ? REPLACED_A : REPLACED_B; } } if (ab == REPLACED_A) *fb = f(*b); if (ab == REPLACED_B) *fa = f(*a); return TERMINATED_BY_USER; }