00001 // This may look like C code, but it is really -*- C++ -*- 00002 /* 00003 ************************************************************************ 00004 * 00005 * Numerical Math Package 00006 * 00007 * Brent's root finder 00008 * obtains a zero of a function of one variable 00009 * 00010 * Synopsis 00011 * double zeroin(ax,bx,f,tol=EPSILON) 00012 * const double ax The root is to be sought within 00013 * const double bx the interval [ax,bx] 00014 * UnivariateFunctor& f The function under consideration 00015 * const double tol Acceptable tolerance for the root 00016 * position. It is an optional parameter 00017 * with default value DBL_EPSILON 00018 * 00019 * Zeroin returns an approximate location for the root with accuracy 00020 * 4*DBL_EPSILON*abs(x) + tol 00021 * 00022 * Algorithm 00023 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 00024 * computations. M., Mir, 1980, p.180 of the Russian edition 00025 * 00026 * The function makes use of a bissection procedure combined with 00027 * a linear or quadratic inverse interpolation. 00028 * At each step the code operates three abscissae - a, b, and c: 00029 * b - the last and the best approximation to the root 00030 * a - the last but one approximation 00031 * c - the last but one or even an earlier approximation such that 00032 * 1) |f(b)| <= |f(c)| 00033 * 2) f(b) and f(c) have opposite signs, i.e. b and c encompass 00034 * the root 00035 * Given these abscissae, the code computes two new approximations, one by the 00036 * bissection procedure and the other one from interpolation (if a,b, and c 00037 * are all different the quadratic interpolation is used, linear otherwise). 00038 * If the approximation obtained by the interpolation looks 00039 * reasonable (i.e. falls within the current interval [b,c], not too close 00040 * to the end points of the interval), the point is accepted as a new 00041 * approximation to the root. Otherwise, the result of the bissection is used. 00042 * Therefore, the range of uncertainty is guaranteed to tighten at 00043 * least by a factor of 1.6 00044 * 00045 * $Id: zeroin.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $ 00046 * 00047 ************************************************************************ 00048 */ 00049 00050 #ifdef __GNUC__ 00051 #pragma implementation "math_num.h" 00052 #endif 00053 #include "math_num.h" 00054 00055 namespace linalg 00056 { 00057 using namespace linalg; 00058 00059 00060 double zeroin( // An estimate to the root 00061 const double ax, // Specify the interval the root 00062 const double bx, // to be sought in 00063 UnivariateFunctor& f, // Function under investigation 00064 const double tol) // Acceptable tolerance 00065 { 00066 assure( tol > 0, "Tolerance must be positive"); 00067 assure( bx > ax, 00068 "Left end point of the interval should be strictly less than the " 00069 "right one" ); 00070 00071 double b = bx; // the last and the best approx to the root 00072 double fb = f(b); 00073 double a = ax; // the last but one approximation 00074 double fa = f(a); 00075 double c = a; // the last but one or even an earlier approx 00076 double fc = fa; // see the condition above 00077 00078 for(;;) // Main iteration loop 00079 { 00080 const double prev_step = b-a; // Step from the previous iteration 00081 00082 if( fabs(fc) < fabs(fb) ) 00083 { // Swap data so that b would be the 00084 a = b; b = c; c = a; // best approximation found so far 00085 fa=fb; fb=fc; fc=fa; 00086 } 00087 // Estimate the effective tolerance 00088 const double tol_act = 2*DBL_EPSILON*fabs(b) + tol/2; 00089 double new_step = (c-b)/2; // Bissection step for this iteration 00090 00091 if( fabs(new_step) <= tol_act || fb == 0 ) 00092 return b; // Acceptable approximation is found 00093 00094 // Figuring out if the interpolation can be tried 00095 if( fabs(prev_step) >= tol_act // If prev_step was large enough 00096 && fabs(fa) > fabs(fb) ) // and was in true direction, 00097 { // Interpolatiom may be tried 00098 00099 double p; // Interpolation step is calcu- 00100 double q; // lated in the form p/q; divi- 00101 // sion operations is delayed 00102 // until the last moment 00103 const double cb = c-b; 00104 00105 if( a==c ) // If we've got only two distinct 00106 { // points linear interpolation 00107 register const double t1 = fb/fa; // can only be applied 00108 p = cb*t1; 00109 q = 1.0 - t1; 00110 } 00111 else // Quadratic inverse interpolation 00112 { 00113 register const double t1=fb/fc, t2=fb/fa; 00114 q = fa/fc; 00115 p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); 00116 q = (q-1.0) * (t1-1.0) * (t2-1.0); 00117 } 00118 00119 if( p > 0 ) // Formulas above computed new_step 00120 q = -q; // = p/q with the wrong sign (on purpose). 00121 else // Correct this, but in such a way so 00122 p = -p; // that p would be positive 00123 00124 if( 2*p < (1.5*cb*q-fabs(tol_act*q)) // If b+p/q falls in [b,c] 00125 && 2*p < fabs(prev_step*q) ) // and isn't too large 00126 new_step = p/q; // it is accepted 00127 // If p/q is too large then the 00128 // bissection procedure can 00129 // reduce [b,c] to a larger 00130 // extent 00131 } 00132 00133 if( fabs(new_step) < tol_act ) // Adjust the step to be not less 00134 new_step = new_step > 0 ? // than the tolerance 00135 tol_act : -tol_act; 00136 00137 a = b; fa = fb; // Save the previous approximation 00138 b += new_step; fb = f(b); // Do step to a new approximation 00139 if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) ) 00140 { // Adjust c for it to have the sign 00141 c = a; fc = fa; // opposite to that of b 00142 } 00143 } 00144 00145 } 00146 }