《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 09.4 Root Finding and Nonlinear Sets of Equations 9.4 Newton-Raphson Method Using Derivative

362 Chapter 9.Root Finding and Nonlinear Sets of Equations 2 a=b: Move last best guess to a. fa-fb; if (fabs(d)>tol1) Evaluate new trial root. b+=d; else b+=SIGN(tol1,xm); fb=(*func)(b); nrerror("Maximum number of iterations exceeded in zbrent"); return 0.0; Never get here. 83g CITED REFERENCES AND FURTHER READING: Brent,R.P.1973,Algorithms for Minimization without Derivatives(Englewood Cliffs,NJ:Prentice- 11800 Hall),Chapters 3,4.[1] Forsythe,G.E.,Malcolm,M.A.,and Moler,C.B.1977,Computer Methods for Mathematical Computations (Englewood Cliffs,NJ:Prentice-Hall),87.2. 9.4 Newton-Raphson Method Using Derivative America Press. Perhaps the most celebrated ofall one-dimensional root-finding routines is New- ton's method,also called the Newton-Raphson method.This method is distinguished from the methods of previous sections by the fact that it requires the evaluation SCIENTIFIC( of both the function f(x),and the derivative f'(z),at arbitrary points x.The Newton-Raphson formula consists geometrically of extending the tangent line at a current pointri until it crosses zero,then setting the next guess to the abscissa of that zero-crossing(see Figure 9.4.1).Algebraically,the method derives from the familiar Taylor series expansion of a function in the neighborhood of a point, fe+)≈f回+fa5+"回+. (9.4.1) 2 Numerica 10621 For small enough values of 6,and for well-behaved functions,the terms beyond 43126 linear are unimportant,hence f(x+6)=0 implies 6s、x) (9.4.2) f'(x) North Newton-Raphson is not restricted to one dimension.The method readily generalizes to multiple dimensions,as we shall see in 89.6 and 89.7.below. Far from a root,where the higher-order terms in the series are important,the Newton-Raphson formula can give grossly inaccurate,meaningless corrections.For instance,the initial guess for the root might be so far from the true root as to let the search interval include a local maximum or minimum of the function.This can be death to the method (see Figure 9.4.2).If an iteration places a trial guess near such a local extremum,so that the first derivative nearly vanishes,then Newton- Raphson sends its solution off to limbo,with vanishingly small hope of recovery

362 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 machine￾readable 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). } a=b; Move last best guess to a. fa=fb; if (fabs(d) > tol1) Evaluate new trial root. b += d; else b += SIGN(tol1,xm); fb=(*func)(b); } 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: Prentice￾Hall), Chapters 3, 4. [1] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §7.2. 9.4 Newton-Raphson Method Using Derivative Perhaps the most celebrated of all one-dimensional root-finding routines is New￾ton’s method, also called the Newton-Raphson method. This method is distinguished from the methods of previous sections by the fact that it requires the evaluation of both the function f(x), and the derivative f
(x), at arbitrary points x. The Newton-Raphson formula consists geometrically of extending the tangent line at a current point xi until it crosses zero, then setting the next guess xi+1 to the abscissa of that zero-crossing (see Figure 9.4.1). Algebraically, the method derives from the familiar 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 beyond linear are unimportant, hence f(x + δ)=0 implies δ = − f(x) f
(x) . (9.4.2) Newton-Raphson is not restricted to one dimension. The method readily generalizes 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, the Newton-Raphson formula can give grossly inaccurate, meaningless corrections. For instance, the initial guess for the root might be so far from the true root as to let the search interval include a local maximum or minimum of the function. This can be death to the method (see Figure 9.4.2). If an iteration places a trial guess near such a local extremum, so that the first derivative nearly vanishes, then Newton￾Raphson sends its solution off to limbo, with vanishingly small hope of recovery

9.4 Newton-Raphson Method Using Derivative 363 f(x) 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.4.1.Newton's method extrapolates the local derivative to find the next estimate of the root.In only),or this example it works well and converges quadratically. own 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) Figure 9.4.2.Unfortunate case where Newton's method encounters a local extremum and shoots off to outer space.Here bracketing bounds,as in rtsafe,would save the day

9.4 Newton-Raphson Method Using Derivative 363 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 machine￾readable 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). 1 2 3 x f(x) Figure 9.4.1. Newton’s method extrapolates the local derivative to find the next estimate of the root. In this example it works well and converges quadratically. f(x) x 1 2 3 Figure 9.4.2. Unfortunate case where Newton’s method encounters a local extremum and shoots off to outer space. Here bracketing bounds, as in rtsafe, would save the day

364 Chapter 9. Root Finding and Nonlinear Sets of Equations f(x) http://www.nr.com or call 1-800-872- Permission is read able files (including this one) granted for from NUMERICAL RECIPES IN C: 19881992 by Cambridge to any server computer, -7423 (North America tusers to make one paper e University Press. THE Figure 9.4.3.Unfortunate case where Newton's method enters a nonconvergent cycle.This behavior ART is often encountered when the function f is obtained,in whole or in part,by table interpolation.With 是 a better initial guess,the method would have succeeded. ictly proh Programs Like most powerful tools,Newton-Raphson can be destructive used in inappropriate circumstances.Figure 9.4.3 demonstrates another possible pathology. Why do we call Newton-Raphson powerful?The answer lies in its rate of to dir convergence:Within a small distance e of z the function and its derivative are approximately: 1881892 OF SCIENTIFIC COMPUTING(ISBN f+)=fa)+efa)+e2f"四 十· 2 (9.4.3) f'(x+e)=f(x)+ef"(x)+… Numerical Recipes 10-521 43108 By the Newton-Raphson formula. f(xi) (outside i+1=i一 f'(xi) (9.4.4) North Software. so that f(xi) Ei+1=Ei- visit website f'(x) (9.4.5) machine When a trial solution x;differs from the true root by ei,we can use(9.4.3)to express f(),f'()in (9.4.4)in terms of e;and derivatives at the root itself.The result is a recurrence relation for the deviations of the trial solutions 2f"(x) 6+1=-6f2f(m (9.4.6)

364 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 machine￾readable 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). x f(x) 2 1 Figure 9.4.3. Unfortunate case where Newton’s method enters a nonconvergent cycle. This behavior is often encountered when the function f is obtained, in whole or in part, by table interpolation. With a better initial guess, the method would have succeeded. Like most powerful tools, Newton-Raphson can be destructive used in inappropriate circumstances. Figure 9.4.3 demonstrates another possible pathology. Why do we call Newton-Raphson powerful? The answer lies in its rate of convergence: Within a small distance  of x the function and its derivative are approximately: f(x + ) = f(x) + f
(x) +  2 f

(x) 2 + ··· , f
(x + ) = f
(x) + f

(x) + ··· (9.4.3) By the Newton-Raphson formula, xi+1 = xi − f(xi) f
(xi) , (9.4.4) so that i+1 = i − f(xi) f
(xi) . (9.4.5) When a trial solution xi differs from the true root by i, we can use (9.4.3) to express f(xi), f
(xi) in (9.4.4) in terms of i and derivatives at the root itself. The result is a recurrence relation for the deviations of the trial solutions i+1 = −2 i f

(x) 2f
(x) . (9.4.6)

9.4 Newton-Raphson Method Using Derivative 365 Equation (9.4.6)says that Newton-Raphson converges quadratically (cf.equa- tion 9.2.3).Near a root,the number of significant digits approximately doubles with each step.This very strong convergence property makes Newton-Raphson the method of choice for any function whose derivative can be evaluated efficiently,and whose 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"polish up"a root with one or two steps of Newton-Raphson,which can multiply by two or four its number of significant figures! For an efficient realization of Newton-Raphson the user provides a routine that 81 evaluates both f(z)and its first derivative f'()at the point z.The Newton-Raphson formula can also be applied using a numerical difference to approximate the true local derivative, f'(x)≈f+)-fa) (9.4.7) dr This is not,however,a recommended procedure for the following reasons:(i)You are doing two function evaluations per step,so at best the superlinear order of convergence will be only v2.(ii)If you take dz too small you will be wiped out by roundoff,while if you take it too large your order of convergence will be only 墨会d的将 Press. linear,no better than using the initial evaluation f'(ro)for all subsequent steps. Therefore,Newton-Raphson with numerical derivatives is(in one dimension)always dominated by the secant method of 89.2.(In multidimensions,where there is a paucity of available methods,Newton-Raphson with numerical derivatives must be taken more seriously.See $89.6-9.7.) The following function calls a user supplied function funcd(x,fn,df)which supplies the function value as fn and the derivative as df.We have included input 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 point 是 x.The bounds are used only to pick the midpoint as the first guess,and to reject the solution if it wanders outside of the bounds 10621 #include <math.h> #define JMAX 20 Set to maximum number of iterations. Fuunrggoleioh Numerical Recipes 43106 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 (outside [x1,x2.The root rtnewt will be refined until its accuracy is known within +xacc.funcd is a user-supplied routine that returns both the function value and the first derivative of the North Software. function at the point x. f void nrerror(char error_text []) int i; 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; 1f((x1-rtn)*(rtn-x2)<0.0) nrerror("Jumped out of brackets in rtnewt");

9.4 Newton-Raphson Method Using Derivative 365 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 machine￾readable 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). Equation (9.4.6) says that Newton-Raphson converges quadratically (cf. equa￾tion 9.2.3). Near a root, the number of significant digits approximately doubles with each step. This very strong convergence property makes Newton-Raphson the method of choice for any function whose derivative can be evaluated efficiently, and whose 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 “polish up” a root with one or two steps of Newton-Raphson, which can multiply by two or four its number of significant figures! For an efficient realization of Newton-Raphson the user provides a routine that evaluates both f(x) and its first derivative f
(x) at the point x. The Newton-Raphson formula can also be applied using a numerical difference to approximate the true local derivative, f
(x) ≈ f(x + dx) − f(x) dx . (9.4.7) This is not, however, a recommended procedure for the following reasons: (i) You are doing two function evaluations per step, so at best the superlinear order of convergence will be only √2. (ii) If you take dx too small you will be wiped out by roundoff, while if you take it too large your order of convergence will be only linear, no better than using the initial evaluation f
(x0) for all subsequent steps. Therefore, Newton-Raphson with numerical derivatives is (in one dimension) always dominated by the secant method of §9.2. (In multidimensions, where there is a paucity of available methods, Newton-Raphson with numerical derivatives must be taken more seriously. See §§9.6–9.7.) The following function calls a user supplied function funcd(x,fn,df) which supplies the function value as fn and the derivative as df. We have included input 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 point x. The bounds are used only to pick the midpoint as the first guess, and to reject the solution if it wanders outside of the bounds. #include <math.h> #define JMAX 20 Set 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. funcd is a user-supplied routine that returns both the function value and the first derivative of the function 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");

366 Chapter 9.Root Finding and Nonlinear Sets of Equations if (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 fairly easy to design a fail-safe routine that utilizes a combination of bisection and Newton- Raphson.The hybrid algorithm takes a bisection step whenever Newton-Raphson would take the solution out of bounds,or whenever Newton-Raphson is not reducing the size of the brackets rapidly enough #include <math.h> 18881892 #define MAXIT 100 Maximum 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 bracketed from NUMERICAL RECIPES I between x1 and x2.The root,returned as the function value rtsafe,will be refined until its accuracy is known within +xacc.funcd is a user-supplied routine that returns both the function value and the first derivative of the function. void nrerror(char error.-text[☐)； University Press. 令 THE int j; Ameri computer one paper float df,dx,dxold,f,fh,fl; ART float temp,xh,xl,rts; (*funcd)(x1,&f1,&df) (*funcd)(x2,&fh,&df) 1f((f1>0.02fh>0.0)I1(f1<0.02fh<0.0)) nrerror("Root must be bracketed in rtsafe"); if (fl =0.0)return x1; if (fh ==0.0)return x2; to dir 1f(f1<0.0){ Orient the search so that f(x1)<0. x1=x1; xh=x2; SCIENTIFIC COMPUTING(ISBN else 18881920 xh=x1; x1=x2; v@cam 10-621 rts=0.5*(x1+x2): Initialize the guess for root. dxold=fabs(x2-x1); the "stepsize before last," Numerical Recipes 43108 dx=dxold; and the last step. (*funcd)(rts.&f,&df): for (i=1;i<=MAXIT;++){ Loop over allowed iterations. if ((((rts-xh)*df-f)*((rts-xl)*df-f)>0.0) Bisect if Newton out of range, (outside 膜 II (fabs(2.0*f)>fabs(dxold*df))){ or not decreasing fast enough Software. dxold=dx: North dx=0.5*(xh-x1); rts=xl+dx; Ame if (xl =rts)return rts; Change in root is negligible. else Newton step acceptable.Take it. dxold=dx; dx=f/df; temp=rts; rts -dx; if (temp =rts)return rts; if (fabs(dx)<xacc)return rts; Convergence criterion. (*funcd)(rts,&f,&df); The one new function evaluation per iteration

366 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 machine￾readable 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). if (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 fairly easy to design a fail-safe routine that utilizes a combination of bisection and Newton￾Raphson. The hybrid algorithm takes a bisection step whenever Newton-Raphson would take the solution out of bounds, or whenever Newton-Raphson is not reducing the size of the brackets rapidly enough. #include <math.h> #define MAXIT 100 Maximum 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 bracketed between x1 and x2. The root, returned as the function value rtsafe, will be refined until its accuracy is known within ±xacc. funcd is a user-supplied routine that returns both the function value and the first derivative of the function. { void nrerror(char error_text[]); int j; float df,dx,dxold,f,fh,fl; float temp,xh,xl,rts; (*funcd)(x1,&fl,&df); (*funcd)(x2,&fh,&df); if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) nrerror("Root must be bracketed in rtsafe"); if (fl == 0.0) return x1; if (fh == 0.0) return x2; if (fl < 0.0) { Orient the search so that f(xl) < 0. xl=x1; xh=x2; } else { xh=x1; xl=x2; } rts=0.5*(x1+x2); Initialize the guess for root, dxold=fabs(x2-x1); the “stepsize before last,” dx=dxold; and the last step. (*funcd)(rts,&f,&df); for (j=1;j<=MAXIT;j++) { Loop over allowed iterations. if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) Bisect if Newton out of range, || (fabs(2.0*f) > fabs(dxold*df))) { or not decreasing fast enough. dxold=dx; dx=0.5*(xh-xl); rts=xl+dx; if (xl == rts) return rts; Change in root is negligible. } else { Newton step acceptable. Take it. dxold=dx; dx=f/df; temp=rts; rts -= dx; if (temp == rts) return rts; } if (fabs(dx) < xacc) return rts; Convergence criterion. (*funcd)(rts,&f,&df); The one new function evaluation per iteration