//////////////////////////////////////////////////////////////////////////////// // File: regula_falsi_terminate.c // // Routines: // // int Regula_Falsi_Terminate // //////////////////////////////////////////////////////////////////////////////// #include // required for fabs() //////////////////////////////////////////////////////////////////////////////// // int Regula_Falsi_Terminate( double a, double fa, double b, double fb, // // unsigned int iteration_count, void *data[], int *halt_code) // // // // Description: // // This routine is an example of a routine to determine whether to // // continue or terminate the routine Regula_Falsi. This routine // // returns a 1 for termination or a 0 for continuation. Further if the // // routine returns a 1, halt, it also returns the reason to halt in // // *halt_code. // // // // In the argument list, a and b, are the endpoints of an interval which // // brackets a zero of the function f(x). The numbers fa and fb are the // // values of the function at a and b respectively. For a continuous // // function f(x) one is assured that the interval contains a zero if the // // signs of fa and fb are different. If the signs of fa and fb are the // // same the function may or may not contain a zero, consider the function // // f(x) = 1 and the function f(x) = x^2 on the interval [-1,1]. // // For Regula Falsi one requires that the signs of fa and fb be different.// // // // The array 'data' contains 5 pointers, the first to the relative // // tolerance, the second to the absolute tolerance, the next two // // to the lower bound and upper bound respectively of the test region // // which is used to determine whether to use the relative tolerance or // // to use the absolute tolerance to test for convergence, and the last // // pointer to the maximum number of iterations allowed. The relative // // tolerance, the absolute tolerance and the upper bound of the test // // region must be strictly positive while the lower bound of the test // // region must be strictly negative and the maximum number of iterations // // allowed is non-negative. // // // // If the interval containing a zero of f(x) is [x0, x1], where x0 < x1, // // and the interval [mu, eta] is the test region where mu < 0 < eta, // // then the iteration process is halted (return 1) with *halt_code set // // to 1 (CONVERGED_ABS) when the [x0, x1] intersects the interval // // [mu, eta] and |x1 - x0| <= absolute tolerance, or with *halt_code set // // to 0 (CONVERGED_REL) when [x0, x1] does not intersect the interval // // [mu,eta] and |x1 - x0| <= min(|x0|,|x1|) * relative tolerance. // // Iteration is also halted when the function value at either endpoint of // // the interval [x0,x1] is zero with *halt_code set to 2 (ZERO) or when // // the iteration count meets or exceeds the maximum allowable number of // // iterations with a *halt_code set to 3 (ITERATION_COUNT_TO_EXCEED_MAX), // // provided terminate() has not already returned a halt. // // // // Arguments: // // double a // // Either the upper or lower bound which brackets a zero of f(x), // // f(a) * f(b) < 0. // // double fa // // The value of f(x) at a. // // double b // // The other bound which brackets a zero of f(x), f(a) * f(b) < 0. // // double fb // // The value of f(x) at b. // // unsigned int iteration_count // // The number of iterations performed. // // void* data[] // // For this routine the array data is defined as follows: // // data[0] is a pointer to the relative tolerance. // // data[1] is a pointer to the absolute tolerance. // // data[2] is a pointer to the test_region_lower_bound. // // data[3] is a pointer to the test_region_upper_bound. // // data[4] is a pointer to the maximum number of iterations. // // // // Return Values: // // The return value is either 0 if the iterative process is to continue // // or 1 if the iterative process is to terminate. // // // // If a 1 is returned, then, for this example, the return codes, // // *halt_code, are: // // 0 - Converged - relative error is less than rel_tolerance. // // 1 - Converged - absolute error is less than abs_tolerance. // // 2 - Either fa or fb is zero. // // 3 - The iteration count meets or exceeds the maximum allowable // // number of iterations. // // // // Example: // // { // // double a,fa,b,fb; // // void *data[]; // // int halt; // // unsigned int iteration_count; // // . // // . // // if ( Regula_Falsi_Terminate(a,fa,b,fb,iteration_count,data,&halt) ) // // return; // // . // // . // // ... // // } // //////////////////////////////////////////////////////////////////////////////// // // enum {CONVERGED_REL, CONVERGED_ABS, ZERO, ITERATION_COUNT_TO_EXCEED_MAX, UNDECIDED}; enum {CONTINUE_ITERATING, HALT}; enum {USE_ABSOLUTE_TOLERANCE, USE_RELATIVE_TOLERANCE}; int Regula_Falsi_Terminate( double a, double fa, double b, double fb, unsigned int iteration_count, void *data[], int *halt_code) { double rel_tolerance = *((double*) (data[0])); double abs_tolerance = *((double*) (data[1])); double test_region_lower_bound = *((double*) (data[2])); double test_region_upper_bound = *((double*) (data[3])); double d; double x0, x1; unsigned int max_number_of_iterations = *((unsigned int*) (data[4])); int convergence_criterion; // In the event that fa = 0.0 or fb = 0.0 set err to ZERO and Return HALT // // In the Regula Falsi routine, some compilers (especially if the // // optimization flags are set) may compute f and leave the result in a // // floating point register, then do the comparison to 0.0 as a long double// // which fails (is not zero), then store the result as a double which is // // then zero. // if ( (fa == 0.0) || (fb == 0.0) ) { *halt_code = ZERO; return HALT; } // Determine convergence criterion - if max(a,b) < test_region_lower_bound // // or if min(a,b) > test_region upper bound then use the relative tolerance // // criterion otherwise use the absolute tolerance criterion. // if (a < b) { x0 = a; x1 = b; } else { x0 = b; x1 = a; } if ((x0 <= test_region_upper_bound) && (x1 >= test_region_lower_bound)) convergence_criterion = USE_ABSOLUTE_TOLERANCE; else convergence_criterion = USE_RELATIVE_TOLERANCE; // Test for convergence // *halt_code = UNDECIDED; switch (convergence_criterion) { case USE_ABSOLUTE_TOLERANCE: if (x0 >= 0.0) { if ((x1 - abs_tolerance) <= x0) *halt_code = CONVERGED_ABS; } else if (x1 <= (abs_tolerance + x0)) *halt_code = CONVERGED_ABS; break; case USE_RELATIVE_TOLERANCE: d = (fabs(a) > fabs(b)) ? fabs(b) : fabs(a); if (d >= 1.0) { if ((fabs( a - b ) / d) <= rel_tolerance) *halt_code = CONVERGED_REL; } else if (fabs( a - b ) <= rel_tolerance * d) *halt_code = CONVERGED_REL; break; default: break; } if (*halt_code == UNDECIDED) if (iteration_count >= max_number_of_iterations) *halt_code = ITERATION_COUNT_TO_EXCEED_MAX; if (*halt_code != UNDECIDED) return HALT; return CONTINUE_ITERATING; }