Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 86
Текст из файла (страница 86)
Ridders’ method does a much better job, but it too cansometimes be fooled. Is there a way to combine superlinear convergence with thesureness of bisection?360Chapter 9.Root Finding and Nonlinear Sets of EquationsYes. We can keep track of whether a supposedly superlinear method is actuallyconverging the way it is supposed to, and, if it is not, we can intersperse bisectionsteps so as to guarantee at least linear convergence. This kind of super-strategyrequires attention to bookkeeping detail, and also careful consideration of howroundoff errors can affect the guiding strategy. Also, we must be able to determinereliably when convergence has been achieved.An excellent algorithm that pays close attention to these matters was developedin the 1960s by van Wijngaarden, Dekker, and others at the Mathematical Centerin Amsterdam, and later improved by Brent [1]. For brevity, we refer to the finalform of the algorithm as Brent’s method. The method is guaranteed (by Brent)to converge, so long as the function can be evaluated within the initial intervalknown to contain a root.Brent’s method combines root bracketing, bisection, and inverse quadraticinterpolation to converge from the neighborhood of a zero crossing.
While the falseposition and secant methods assume approximately linear behavior between twoprior root estimates, inverse quadratic interpolation uses three prior points to fit aninverse quadratic function (x as a quadratic function of y) whose value at y = 0 istaken as the next estimate of the root x.
Of course one must have contingency plansfor what to do if the root falls outside of the brackets. Brent’s method takes care ofall that. If the three point pairs are [a, f(a)], [b, f(b)], [c, f(c)] then the interpolationformula (cf. equation 3.1.1) isx=[y − f(a)][y − f(b)]c[y − f(b)][y − f(c)]a+[f(c) − f(a)][f(c) − f(b)] [f(a) − f(b)][f(a) − f(c)][y − f(c)][y − f(a)]b+[f(b) − f(c)][f(b) − f(a)](9.3.1)Setting y to zero gives a result for the next root estimate, which can be written asx = b + P/Q(9.3.2)where, in terms ofR ≡ f(b)/f(c),S ≡ f(b)/f(a),T ≡ f(a)/f(c)(9.3.3)we haveP = S [T (R − T )(c − b) − (1 − R)(b − a)]Q = (T − 1)(R − 1)(S − 1)(9.3.4)(9.3.5)In practice b is the current best estimate of the root and P/Q ought to be a “small”correction. Quadratic methods work well only when the function behaves smoothly;they run the serious risk of giving very bad estimates of the next root or causingmachine failure by an inappropriate division by a very small number (Q ≈ 0).Brent’s method guards against this problem by maintaining brackets on the rootand checking where the interpolation would land before carrying out the division.When the correction P/Q would not land within the bounds, or when the boundsare not collapsing rapidly enough, the algorithm takes a bisection step.
Thus,9.3 Van Wijngaarden–Dekker–Brent Method361Brent’s method combines the sureness of bisection with the speed of a higher-ordermethod when appropriate. We recommend it as the method of choice for generalone-dimensional root finding where a function’s values only (and not its derivativeor functional form) are available.#include <math.h>#include "nrutil.h"#define ITMAX 100#define EPS 3.0e-8Maximum allowed number of iterations.Machine floating-point precision.float zbrent(float (*func)(float), float x1, float x2, float tol)Using Brent’s method, find the root of a function func known to lie between x1 and x2. Theroot, returned as zbrent, will be refined until its accuracy is tol.{int iter;float a=x1,b=x2,c=x2,d,e,min1,min2;float fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm;if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))nrerror("Root must be bracketed in zbrent");fc=fb;for (iter=1;iter<=ITMAX;iter++) {if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {c=a;Rename a, b, c and adjust bounding intervalfc=fa;d.e=d=b-a;}if (fabs(fc) < fabs(fb)) {a=b;b=c;c=a;fa=fb;fb=fc;fc=fa;}tol1=2.0*EPS*fabs(b)+0.5*tol;Convergence check.xm=0.5*(c-b);if (fabs(xm) <= tol1 || fb == 0.0) return b;if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {s=fb/fa;Attempt inverse quadratic interpolation.if (a == c) {p=2.0*xm*s;q=1.0-s;} else {q=fa/fc;r=fb/fc;p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));q=(q-1.0)*(r-1.0)*(s-1.0);}if (p > 0.0) q = -q;Check whether in bounds.p=fabs(p);min1=3.0*xm*q-fabs(tol1*q);min2=fabs(e*q);if (2.0*p < (min1 < min2 ? min1 : min2)) {e=d;Accept interpolation.d=p/q;} else {d=xm;Interpolation failed, use bisection.e=d;}} else {Bounds decreasing too slowly, use bisection.d=xm;e=d;362Chapter 9.Root Finding and Nonlinear Sets of Equations}a=b;fa=fb;if (fabs(d) > tol1)b += d;elseb += SIGN(tol1,xm);fb=(*func)(b);Move last best guess to a.Evaluate new trial root.}nrerror("Maximum number of iterations exceeded in zbrent");return 0.0;Never get here.}CITED REFERENCES AND FURTHER READING:Brent, R.P.
1973, Algorithms for Minimization without Derivatives (Englewood Cliffs, NJ: PrenticeHall), Chapters 3, 4. [1]Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §7.2.9.4 Newton-Raphson Method Using DerivativePerhaps the most celebrated of all one-dimensional root-finding routines is Newton’s method, also called the Newton-Raphson method. This method is distinguishedfrom the methods of previous sections by the fact that it requires the evaluationof both the function f(x), and the derivative f (x), at arbitrary points x. TheNewton-Raphson formula consists geometrically of extending the tangent line at acurrent point xi until it crosses zero, then setting the next guess xi+1 to the abscissaof that zero-crossing (see Figure 9.4.1).
Algebraically, the method derives from thefamiliar Taylor series expansion of a function in the neighborhood of a point,f(x + δ) ≈ f(x) + f (x)δ +f (x) 2δ + ....2(9.4.1)For small enough values of δ, and for well-behaved functions, the terms beyondlinear are unimportant, hence f(x + δ) = 0 impliesδ=−f(x).f (x)(9.4.2)Newton-Raphson is not restricted to one dimension.
The method readilygeneralizes to multiple dimensions, as we shall see in §9.6 and §9.7, below.Far from a root, where the higher-order terms in the series are important, theNewton-Raphson formula can give grossly inaccurate, meaningless corrections. Forinstance, the initial guess for the root might be so far from the true root as to letthe search interval include a local maximum or minimum of the function. This canbe death to the method (see Figure 9.4.2).
If an iteration places a trial guess nearsuch a local extremum, so that the first derivative nearly vanishes, then NewtonRaphson sends its solution off to limbo, with vanishingly small hope of recovery.3639.4 Newton-Raphson Method Using Derivativef(x)123xFigure 9.4.1. Newton’s method extrapolates the local derivative to find the next estimate of the root. Inthis example it works well and converges quadratically.f(x)321xFigure 9.4.2. Unfortunate case where Newton’s method encounters a local extremum and shoots off toouter space.
Here bracketing bounds, as in rtsafe, would save the day.364Chapter 9.Root Finding and Nonlinear Sets of Equationsf (x)1x2Figure 9.4.3. Unfortunate case where Newton’s method enters a nonconvergent cycle. This behavioris often encountered when the function f is obtained, in whole or in part, by table interpolation.
Witha better initial guess, the method would have succeeded.Like most powerful tools, Newton-Raphson can be destructive used in inappropriatecircumstances. Figure 9.4.3 demonstrates another possible pathology.Why do we call Newton-Raphson powerful? The answer lies in its rate ofconvergence: Within a small distance of x the function and its derivative areapproximately:f (x)+ · · ·,2f (x + ) = f (x) + f (x) + · · ·f(x + ) = f(x) + f (x) + 2(9.4.3)By the Newton-Raphson formula,xi+1 = xi −f(xi ),f (xi )(9.4.4)i+1 = i −f(xi ).f (xi )(9.4.5)so thatWhen a trial solution xi differs from the true root by i , we can use (9.4.3) to expressf(xi ), f (xi ) in (9.4.4) in terms of i and derivatives at the root itself. The result isa recurrence relation for the deviations of the trial solutionsi+1 = −2if (x).2f (x)(9.4.6)9.4 Newton-Raphson Method Using Derivative365Equation (9.4.6) says that Newton-Raphson converges quadratically (cf.
equation 9.2.3). Near a root, the number of significant digits approximately doubleswith each step. This very strong convergence property makes Newton-Raphson themethod of choice for any function whose derivative can be evaluated efficiently, andwhose derivative is continuous and nonzero in the neighborhood of a root.Even where Newton-Raphson is rejected for the early stages of convergence(because of its poor global convergence properties), it is very common to “polishup” a root with one or two steps of Newton-Raphson, which can multiply by twoor four its number of significant figures!For an efficient realization of Newton-Raphson the user provides a routine thatevaluates both f(x) and its first derivative f (x) at the point x.
The Newton-Raphsonformula can also be applied using a numerical difference to approximate the truelocal derivative,f (x) ≈f(x + dx) − f(x).dx(9.4.7)This is not, however, a recommended procedure for the following reasons: (i) Youare doing two function evaluationsper step, so at best the superlinear order of√convergence will be only 2. (ii) If you take dx too small you will be wipedout by roundoff, while if you take it too large your order of convergence will beonly linear, no better than using the initial evaluation f (x0 ) for all subsequentsteps. Therefore, Newton-Raphson with numerical derivatives is (in one dimension)always dominated by the secant method of §9.2.
(In multidimensions, where thereis a paucity of available methods, Newton-Raphson with numerical derivatives mustbe taken more seriously. See §§9.6–9.7.)The following function calls a user supplied function funcd(x,fn,df) whichsupplies the function value as fn and the derivative as df. We have includedinput bounds on the root simply to be consistent with previous root-finding routines:Newton does not adjust bounds, and works only on local information at the pointx. The bounds are used only to pick the midpoint as the first guess, and to rejectthe solution if it wanders outside of the bounds.#include <math.h>#define JMAX 20Set to maximum number of iterations.float rtnewt(void (*funcd)(float, float *, float *), float x1, float x2,float xacc)Using the Newton-Raphson method, find the root of a function known to lie in the interval[x1, x2].
The root rtnewt will be refined until its accuracy is known within ±xacc. funcdis a user-supplied routine that returns both the function value and the first derivative of thefunction at the point x.{void nrerror(char error_text[]);int j;float df,dx,f,rtn;rtn=0.5*(x1+x2);Initial guess.for (j=1;j<=JMAX;j++) {(*funcd)(rtn,&f,&df);dx=f/df;rtn -= dx;if ((x1-rtn)*(rtn-x2) < 0.0)nrerror("Jumped out of brackets in rtnewt");366Chapter 9.Root Finding and Nonlinear Sets of Equationsif (fabs(dx) < xacc) return rtn;Convergence.}nrerror("Maximum number of iterations exceeded in rtnewt");return 0.0;Never get here.}While Newton-Raphson’s global convergence properties are poor, it is fairlyeasy to design a fail-safe routine that utilizes a combination of bisection and NewtonRaphson. The hybrid algorithm takes a bisection step whenever Newton-Raphsonwould take the solution out of bounds, or whenever Newton-Raphson is not reducingthe size of the brackets rapidly enough.#include <math.h>#define MAXIT 100Maximum allowed number of iterations.float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,float xacc)Using a combination of Newton-Raphson and bisection, find the root of a function bracketedbetween x1 and x2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
