354 Chapter 9. Root Finding and Nonlinear Sets of Equations #include #define JMAX 40 Maximum allowed number of bisections. float rtbis(float (*func)(float),float xi,float x2,float xacc) Using bisection,find the root of a function func known to lie between x1 and x2.The root, returned as rtbis,will be refined until its accuracy is txacc. void nrerror(char error-text☐)； int ji float dx,f,fmid,xmid,rtb; f=(*func)(x1): 三 fmid=(*func)(x2); if (f*fmid >=0.0)nrerror("Root must be bracketed for bisection in rtbis"); rtb=f0 for (j=1;j<=JMAX;++){ lies at x+dx. fmid=(*func)(xmid=rtb+(dx *0.5)); Bisection loop. if (fmid <=0.0)rtb-xmid; if (fabs(dx)<xacc II fmid ==0.0)return rtb; ICAL nrerror("Too many bisections in rtbis"); return 0.0; Never get here. RECIPES I America computer, Press. 9 9.2 Secant Method,False Position Method, Programs and Ridders'Method IENTIFIC For functions that are smooth near a root,the methods known respectively to dir as false position (or regula falsi)and secant method generally converge faster than bisection.In both of these methods the function is assumed to be approximately linear in the local region of interest,and the next improvement in the root is taken as the point where the approximating line crosses the axis.After each iteration one of the previous boundary points is discarded in favor of the latest estimate of the root. Recipes 10621 The only difference between the methods is that secant retains the most recent Numerica of the prior estimates(Figure 9.2.1;this requires an arbitrary choice on the first 43106 iteration),while false position retains that prior estimate for which the function value Recipes has opposite sign from the function value at the current best estimate of the root, (outside 腿 so that the two points continue to bracket the root(Figure 9.2.2).Mathematically, the secant method converges more rapidly near a root of a sufficiently continuous North function.Its order of convergence can be shown to be the"golden ratio"1.618..., so that lim+const x (9.2.1) The secant method has,however,the disadvantage that the root does not necessarily remain bracketed.For functions that are not sufficiently continuous.the algorithm can therefore not be guaranteed to converge:Local behavior might send it off towards infinity
354 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #define JMAX 40 Maximum allowed number of bisections. float rtbis(float (*func)(float), float x1, float x2, float xacc) Using bisection, find the root of a function func known to lie between x1 and x2. The root, returned as rtbis, will be refined until its accuracy is ±xacc. { void nrerror(char error_text[]); int j; float dx,f,fmid,xmid,rtb; f=(*func)(x1); fmid=(*func)(x2); if (f*fmid >= 0.0) nrerror("Root must be bracketed for bisection in rtbis"); rtb = f 0 for (j=1;j<=JMAX;j++) { lies at x+dx. fmid=(*func)(xmid=rtb+(dx *= 0.5)); Bisection loop. if (fmid <= 0.0) rtb=xmid; if (fabs(dx) < xacc || fmid == 0.0) return rtb; } nrerror("Too many bisections in rtbis"); return 0.0; Never get here. } 9.2 Secant Method, False Position Method, and Ridders’ Method For functions that are smooth near a root, the methods known respectively as false position (or regula falsi) and secant method generally converge faster than bisection. In both of these methods the function is assumed to be approximately linear in the local region of interest, and the next improvement in the root is taken as the point where the approximating line crosses the axis. After each iteration one of the previous boundary points is discarded in favor of the latest estimate of the root. The only difference between the methods is that secant retains the most recent of the prior estimates (Figure 9.2.1; this requires an arbitrary choice on the first iteration), while false position retains that prior estimate for which the function value has opposite sign from the function value at the current best estimate of the root, so that the two points continue to bracket the root (Figure 9.2.2). Mathematically, the secant method converges more rapidly near a root of a sufficiently continuous function. Its order of convergence can be shown to be the “golden ratio” 1.618 ... , so that lim k→∞ |k+1| ≈ const × |k| 1.618 (9.2.1) The secant method has, however, the disadvantage that the root does not necessarily remain bracketed. For functions that are not sufficiently continuous, the algorithm can therefore not be guaranteed to converge: Local behavior might send it off towards infinity.
9.2 Secant Method.False Position Method.and Ridders'Method 355 f(x http://www.nr.com or call 1-800-872-7423 (North America Permission is granted for internet users to make one paper copy for their Figure 9.2.1.Secant method.Extrapolation or interpolation lines (dashed)are drawn through the two only),or most recently evaluated points,whether or not they bracket the function.The points are numbered in the order that they are used 2 f八x) rsend email to directcustserv@cambridge.org(outside North America) Copyright (C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) 3 Figure 9.2.2.False position method.Interpolation lines (dashed)are drawn through the most recent points that bracket the root.In this example,point I thus remains "active"for many steps.False position converges less rapidly than the secant method,but it is more certain
9.2 Secant Method, False Position Method, and Ridders’ Method 355 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). f(x) 2 3 4 1 x Figure 9.2.1. Secant method. Extrapolation or interpolation lines (dashed) are drawn through the two most recently evaluated points, whether or not they bracket the function. The points are numbered in the order that they are used. f(x) x 4 3 2 1 Figure 9.2.2. False position method. Interpolation lines (dashed) are drawn through the most recent points that bracket the root. In this example, point 1 thus remains “active” for many steps. False position converges less rapidly than the secant method, but it is more certain
356 Chapter 9. Root Finding and Nonlinear Sets of Equations f(x) http://www.nr. Permission is read able files .com or call (including this one) granted fori 19881992 11-800-872 134 (North America to any server computer, one paper Cambridge University Press. tusers to make from NUMERICAL RECIPES IN THE Figure 9.2.3.Example where both the secant and false position methods will take many iterations to 是 ART arrive at the true root.This function would be difficult for mamy other root-finding methods. Programs False position,since it sometimes keeps an older rather than newer function evaluation,has a lower order of convergence.Since the newer function value will sometimes be kept,the method is often superlinear,but estimation of its exact order is not so easy. to dir Here are sample implementations of these two related methods.While these methods are standard textbook fare,Ridders'method,described below,or Brent's method.in the next section.are almost always better choices.Figure 9.2.3 shows the 1788-1982 OF SCIENTIFIC COMPUTING(ISBN behavior of secant and false-position methods in a difficult situation. #include #define MAXIT 30 Set to the maximum allowed number of iterations. Further reproduction Numerical Recipes 10-621 43108 float rtflsp(float (*func)(float),float x1,float x2,float xacc) Using the false position method,find the root of a function func known to lie between x1 and x2.The root,returned as rtflsp.is refined until its accuracy is +xacc. (outside void nrerror(char error._text[☐)； North Software. int j; float fl,fh,xl,xh,swap,dx,del,f,rtf; f1=(*func)(x1); visit website fh=(*func)(x2) Be sure the interval brackets a root. machine if(f1*fh>0.0) nrerror("Root must be bracketed in rtflsp"); 1f(f1<0.0) Identify the limits so that xl corresponds to the low X=X1: side. xh=x2; else x1=x2; xh=x1; swap=fl;
356 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 2 f(x) 134 x Figure 9.2.3. Example where both the secant and false position methods will take many iterations to arrive at the true root. This function would be difficult for many other root-finding methods. False position, since it sometimes keeps an older rather than newer function evaluation, has a lower order of convergence. Since the newer function value will sometimes be kept, the method is often superlinear, but estimation of its exact order is not so easy. Here are sample implementations of these two related methods. While these methods are standard textbook fare, Ridders’ method, described below, or Brent’s method, in the next section, are almost always better choices. Figure 9.2.3 shows the behavior of secant and false-position methods in a difficult situation. #include #define MAXIT 30 Set to the maximum allowed number of iterations. float rtflsp(float (*func)(float), float x1, float x2, float xacc) Using the false position method, find the root of a function func known to lie between x1 and x2. The root, returned as rtflsp, is refined until its accuracy is ±xacc. { void nrerror(char error_text[]); int j; float fl,fh,xl,xh,swap,dx,del,f,rtf; fl=(*func)(x1); fh=(*func)(x2); Be sure the interval brackets a root. if (fl*fh > 0.0) nrerror("Root must be bracketed in rtflsp"); if (fl < 0.0) { Identify the limits so that xl corresponds to the low xl=x1; side. xh=x2; } else { xl=x2; xh=x1; swap=fl;
9.2 Secant Method,False Position Method,and Ridders'Method 357 fl=fh; fh=swapi dx=xh-xl; for (j=1;j #define MAXIT 30 Maximum allowed number of iterations 是 Programs float rtsec(float (*func)(float),float x1,float x2,float xacc) Using the secant method,find the root of a function func thought to lie between x1 and x2. The root,returned as rtsec,is refined until its accuracy is +xacc. void nrerror(char error._text[☐)： to dir int j; float fl,f,dx,swap,xl,rts; f1=(*func)(x1); ectcustser f=(*func)(x2); ART OF SCIENTIFIC COMPUTING(ISBN 0-521 if (fabs(fl)<fabs(f)){ Pick the bound with the smaller function value as rtsexI the most recent guess. X1=X2; swap=fl; 8nwg f1=f; Numerical Recipes books or 1988-1992 by Numerical Recipes f=swapi -431085 else x1=x1; rts=x2; Software. for (j=1;j<=MAXIT;j++){ Secant loop. dx=(xl-rts)*f/(f-fl) Increment with respect to latest value (outside North America) xlerts; f1=f; rts +dx; f=(*func)(rts): if (fabs(dx)<xacc Il f ==0.0)return rts; Convergence. nrerror("Maximum number of iterations exceeded in rtsec"); return 0.0; Never get here
9.2 Secant Method, False Position Method, and Ridders’ Method 357 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). fl=fh; fh=swap; } dx=xh-xl; for (j=1;j #define MAXIT 30 Maximum allowed number of iterations. float rtsec(float (*func)(float), float x1, float x2, float xacc) Using the secant method, find the root of a function func thought to lie between x1 and x2. The root, returned as rtsec, is refined until its accuracy is ±xacc. { void nrerror(char error_text[]); int j; float fl,f,dx,swap,xl,rts; fl=(*func)(x1); f=(*func)(x2); if (fabs(fl) < fabs(f)) { Pick the bound with the smaller function value as rts=x1; the most recent guess. xl=x2; swap=fl; fl=f; f=swap; } else { xl=x1; rts=x2; } for (j=1;j<=MAXIT;j++) { Secant loop. dx=(xl-rts)*f/(f-fl); Increment with respect to latest value. xl=rts; fl=f; rts += dx; f=(*func)(rts); if (fabs(dx) < xacc || f == 0.0) return rts; Convergence. } nrerror("Maximum number of iterations exceeded in rtsec"); return 0.0; Never get here. }
358 Chapter 9.Root Finding and Nonlinear Sets of Equations Ridders'Method A powerful variant on false position is due to Ridders [1].When a root is bracketed between z1 and z2,Ridders'method first evaluates the function at the midpoint x3=(x1+x2)/2.It then factors out that unique exponential function which turns the residual function into a straight line.Specifically,it solves for a factor e that gives f(z1)-2f(a3)eQ+f(x2)e2Q=0 (9.2.2) This is a quadratic equation in e,which can be solved to give eQ=f(a)+signlf(2)vf(r3)2-f(1)f(2) (9.2.3) f(c2) Now the false position method is applied,not to the values f(1),f(3),f(2),but to the values f(1),f(x3)e,f(2)e20,yielding a new guess for the root,4.The overall updating formula (incorporating the solution 9.2.3)is 、是只0个 令 4=3+(-r)sgm）-f2fa） (9.2.4) Vf(3)2-f(1)f(x2) Equation(9.2.4)has some very nice properties.First,x4 is guaranteed to lie in the interval(1,2),so the method never jumps out of its brackets.Second, the convergence of successive applications of equation(9.2.4)is guadratic,that is, m 2 in equation (9.1.4).Since each application of (9.2.4)requires two function evaluations,the actual order of the method is v2,not 2;but this is still quite 可 respectably superlinear:the number of significant digits in the answer approximately doubles with each two function evaluations.Third,taking out the function's"bend" via exponential(that is,ratio)factors,rather than via a polynomial technique(e.g., fitting a parabola),turns out to give an extraordinarily robust algorithm.In both reliability and speed,Ridders'method is generally competitive with the more highly developed and better established(but more complicated)method of Van Wijngaarden, Dekker,and Brent,which we next discuss. Fuurggoglrion Numerical Recipes 10.621 43106 #include #include "nrutil.h" #define MAXIT 60 (outside #define UNUSED (-1.11e30) Software. float zriddr(float (*func)(float),float x1,float x2,float xacc) Using Ridders'method,return the root of a function func known to lie between x1 and x2 The root,returned as zriddr,will be refined to an approximate accuracy xacc. int j; float ans,fh,f1,fm,fne,s,xh,x1,x知，Xnew; f1=(*func)(x1); fh=(*func)(x2) if(f1>0.0fh0.0)[ x1=x1; xh=x2; ans-UNUSED; Any highly unlikely value,to simplify logic below
358 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Ridders’ Method A powerful variant on false position is due to Ridders [1]. When a root is bracketed between x1 and x2, Ridders’ method first evaluates the function at the midpoint x3 = (x1 + x2)/2. It then factors out that unique exponential function which turns the residual function into a straight line. Specifically, it solves for a factor eQ that gives f(x1) − 2f(x3)eQ + f(x2)e2Q =0 (9.2.2) This is a quadratic equation in eQ, which can be solved to give eQ = f(x3) + sign[f(x2)]f(x3)2 − f(x1)f(x2) f(x2) (9.2.3) Now the false position method is applied, not to the values f(x1), f(x3), f(x2), but to the values f(x1), f(x3)eQ, f(x2)e2Q, yielding a new guess for the root, x4. The overall updating formula (incorporating the solution 9.2.3) is x4 = x3 + (x3 − x1) sign[f(x1) − f(x2)]f(x3) f(x3)2 − f(x1)f(x2) (9.2.4) Equation (9.2.4) has some very nice properties. First, x4 is guaranteed to lie in the interval (x1, x2), so the method never jumps out of its brackets. Second, the convergence of successive applications of equation (9.2.4) is quadratic, that is, m = 2 in equation (9.1.4). Since each application of (9.2.4) requires two function evaluations, the actual order of the method is √ 2, not 2; but this is still quite respectably superlinear: the number of significant digits in the answer approximately doubles with each two function evaluations. Third, taking out the function’s “bend” via exponential (that is, ratio) factors, rather than via a polynomial technique (e.g., fitting a parabola), turns out to give an extraordinarily robust algorithm. In both reliability and speed, Ridders’ method is generally competitive with the more highly developed and better established (but more complicated) method of Van Wijngaarden, Dekker, and Brent, which we next discuss. #include #include "nrutil.h" #define MAXIT 60 #define UNUSED (-1.11e30) float zriddr(float (*func)(float), float x1, float x2, float xacc) Using Ridders’ method, return the root of a function func known to lie between x1 and x2. The root, returned as zriddr, will be refined to an approximate accuracy xacc. { int j; float ans,fh,fl,fm,fnew,s,xh,xl,xm,xnew; fl=(*func)(x1); fh=(*func)(x2); if ((fl > 0.0 && fh 0.0)) { xl=x1; xh=x2; ans=UNUSED; Any highly unlikely value, to simplify logic below
9.3 Van Wiingaarden-Dekker-Brent Method 359 for (j=1;j=fh?1.0:-1.0)*fm/s); Updating formula. if (fabs(xnew-ans)<=xacc)return ans; ans=xnew; fnew=(*func)(ans); Second of two function evaluations per if (fnew ==0.0)return ans; iteration. if (SIGN(fm,fnew)!fm){ Bookkeeping to keep the root bracketed X1=X; on next iteration. fl=fm; xh=ans; fh=fnew; else if (SIGN(fl,fnew)!=fl){ 19881992 xh=ans; fh=fnew; 1.800 else if (SIGN(fh,fnew)!fh){ xl=ans; 872 fl-fnew; Cambridge from NUMERICAL RECIPESI else nrerror("never get here."); if (fabs(xh-xl)<=xacc)return ans; nrerror("zriddr exceed maximum iterations"); else (North America server computer, sto make one paper University Press. IN C: THE if (fl==0.0)return x1; ART if (fh =0.0)return x2; nrerror("root must be bracketed in zriddr."); Programs return 0.0; Never get here. strictly proh to dir Copyright (C) CITED REFERENCES AND FURTHER READING: OF SCIENTIFIC COMPUTING(ISBN Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerica/Analysis,2nd ed.(New York: 198819920 McGraw-Hill),88.3. Ostrowski,A.M.1966,Solutions of Equations and Systems of Equations,2nd ed.(New York: Academic Press),Chapter 12. Numerical Recipes 10-621 Ridders,C.J.F.1979,IEEE Transactions on Circuits and Systems,vol.CAS-26,pp.979-980.[1] 43108 (outside 9.3 Van Wijngaarden-Dekker-Brent Method North Software. Ame While secant and false position formally converge faster than bisection,one visit website finds in practice pathological functions for which bisection converges more rapidly. machine These can be choppy,discontinuous functions,or even smooth functions if the second derivative changes sharply near the root.Bisection always halves the interval, while secant and false position can sometimes spend many cycles slowly pulling distant bounds closer to a root.Ridders'method does a much better job,but it too can sometimes be fooled.Is there a way to combine superlinear convergence with the sureness of bisection?
9.3 Van Wijngaarden–Dekker–Brent Method 359 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (j=1;j= fh ? 1.0 : -1.0)*fm/s); Updating formula. if (fabs(xnew-ans) <= xacc) return ans; ans=xnew; fnew=(*func)(ans); Second of two function evaluations per if (fnew == 0.0) return ans; iteration. if (SIGN(fm,fnew) != fm) { Bookkeeping to keep the root bracketed xl=xm; on next iteration. fl=fm; xh=ans; fh=fnew; } else if (SIGN(fl,fnew) != fl) { xh=ans; fh=fnew; } else if (SIGN(fh,fnew) != fh) { xl=ans; fl=fnew; } else nrerror("never get here."); if (fabs(xh-xl) <= xacc) return ans; } nrerror("zriddr exceed maximum iterations"); } else { if (fl == 0.0) return x1; if (fh == 0.0) return x2; nrerror("root must be bracketed in zriddr."); } return 0.0; Never get here. } CITED REFERENCES AND FURTHER READING: Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York: McGraw-Hill), §8.3. Ostrowski, A.M. 1966, Solutions of Equations and Systems of Equations, 2nd ed. (New York: Academic Press), Chapter 12. Ridders, C.J.F. 1979, IEEE Transactions on Circuits and Systems, vol. CAS-26, pp. 979–980. [1] 9.3 Van Wijngaarden–Dekker–Brent Method While secant and false position formally converge faster than bisection, one finds in practice pathological functions for which bisection converges more rapidly. These can be choppy, discontinuous functions, or even smooth functions if the second derivative changes sharply near the root. Bisection always halves the interval, while secant and false position can sometimes spend many cycles slowly pulling distant bounds closer to a root. Ridders’ method does a much better job, but it too can sometimes be fooled. Is there a way to combine superlinear convergence with the sureness of bisection?