Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 30
Текст из файла (страница 30)
Compare with Exercise 4.1.4.3 Geometrically estimate the root of the function F(x)whose graph is given below.(a) For an initial bracket of [0.0, 1.0] what are the nextthree brackets using the bisection method on this func-150CHAPTER 4 ROOTS OF NONLINEAR EQUATIONStion?(b) If xl = 0.0 and x2 = 1.0, mark on the graph theapproximate location of x3 using one step of the secant method.(c) If x1 = 0.5, mark on the graph the approximatelocation of x2 and x3 using two steps of Newton’smethod.linear polynomial that interpolates G(y) at f(xi) andf(xi-l), show that P(0) = xi+1 is given by (4.6).4.7 In the convergence proof for the secant rule, it wasstated that if ε = max(ε0, ε1) < 1, then the inequalityimplied34.4 The polynomial f(x) = x - 2x - 5 has a root a in[2,3].(a) Show that [2,3] is a bracket for f(x).and(b) Apply four steps of the bisection method to reduce the width of the bracket to 1/16.(c) Calculate x3 and x4 by the secant method startingwith xl = 3 and x2 = 2.Establish this.4.8 The special function(d) Calculate x2, x3, and x4 using Newton’s methodwith x1 = 2.4.5 To find where sinx = x/2 for x > 0,(a) find a bracket for an appropriate function f(x).(b) Apply four steps of the bisection method to reduce the width of the bracket by 1/16.(c) Calculate x3 and x4 by the secant method startingwith x1 and x2 equal to the bracket values.(d) Calculate x2, x3, and x4 using Newton’s methodwith x1 the midpoint of your bracket.4.6 Given x0 and x1, show that the inverse secant rule,discussed at the end of Section 4.1, and the direct secant rule (4.6) produce the same iterates.
In particular,with G(y) the inverse function of f(x), and P(y) the4.2is important in statistics and in many areas of scienceand engineering. Because the integrand is positive forall t, the function is strictly increasing and so has aninverse x = erf-1(y). The inverse error function is animportant function in its own right and can be evaluated for given y by solving the equation y = erf(x).The algorithm of the MATLAB function erfinv.mfirst forms a rational approximation to y that is accurate to about six figures.
Two Newton iterationsare then done to get a result to full accuracy. Whatis Newton’s method for solving this equation? Whywould you expect two iterations to be enough? (Don’tforget to consider the multiplicity of the root.)AN ALGORITHM COMBINING BISECTION AND THE SECANT RULEIt is a challenging task to fuse several methods into an efficient computational scheme.This section is devoted to a code, Zero, based on one written by Dekker [6] that doesthis.
Roughly speaking, the code uses the secant rule unless bisection appears advantageous. A very similar code is found in the NAG library. Brent [3] added the use ofquadratic inverse interpolation to Dekker’s ideas. Brent’s code is the basis for codes inMATLAB and the IMSL library.Normal input to Zero is a subprogram for the evaluation of a continuous functionf(x) and arguments B and C for which f(B)f(C) < 0. Throughout the computationB and C are end points of an interval with f(B)f(C) < 0 that is decreasing in length.In favorable circumstances B is computed with the secant rule and is a much betterapproximate root than either C or the midpoint M = ( C + B)/2 of the interval. To4.2 AN ALGORITHM COMBINING BISECTION AND THE SECANT RULE151deal with unfavorable circumstances, the code interchanges the values of B and C asnecessary so that |f(B)| < |f(C)| holds.
If at any time the computed f(B) is zero, thecomputation is terminated and B is reported as a root.The convergence test is a mixed relative-absolute error test. Two parameters ABSERR and RELERR are input and it is asked of each iterate whether(4.13)For reasons discussed in Chapter 1, the code will not permit RELERR to be smallerthan 10 units of roundoff, nor ABSERR to be zero. However, to understand what thetest means, first suppose that RELERR is zero.
The test is then asking if an intervalbelieved to contain a root has a length no more than 2 × ABSERR. If so, the midpointM is no farther than ABSERR from a root and this is a pure absolute error test onM as an approximate root. However, it is believed that the quantity B reported as theanswer is closer to a root than M is. Even if it is not, the test implies that B is within2 × ABSERR of a root. Similarly, if the parameter ABSERR is zero and if the testwerethe test would be a pure relative error test for the approximate root M.
Because it isbelieved that B is a better approximate root, it is used in the test rather than M. Thefzero.m function in MATLAB has a similar, but somewhat simpler, test. The codes inthe NAG and lMSL libraries have convergence tests that are broadly similar, but theyalso test the size of the residual and convergence can occur either because the root hasbeen located to a specified accuracy or because the magnitude of the residual is smallerthan a specified value.Unless there is a reason to do otherwise, Zero uses the secant rule. A variable A isinitialized to C. The two variables A, B are the two iterates used by the secant rule tocalculateA danger with the secant rule (and Newton’s method) is an interpolant that is horizontalor nearly so.
The extreme case is a division by zero in this formula. This danger isavoided by requiring D to lie in the interval [B,C] known to contain a root and checkingthis without performing the division. Pursuing the tactic further, the code requires thatD lie in [B,M] on the grounds that B ought to be a better approximation to the root thanC, so if the secant rule is working properly, D ought to be closer to B than to C.
If Ddoes not lie in [B,M], the midpoint M is used as the next iterate.The performance of the code can be improved in some circumstances by selectingan iterate in a different way. If D is too close to B, a better tactic is to move a minimumdistance away from B. The quantity max[ABSERR, |B| × RELERR] is called TOL inthe code.
If |D - B| < TOL, then the value B + TOL × sign(C - B) is used instead ofD. This choice cannot result in an iterate outside the interval [B,C] since |B - C| >2 × TOL (or else the error test would have been passed). If the root a is further from152CHAPTER 4ROOTS OF NONLINEAR EQUATIONSB than TOL, the iterate chosen in this way is closer to the root than D. If it is closer,this iterate and B will bracket the root and the code will converge at the next test onthe error because the length of the bracket is TOL.There are circumstances in which the current iterate B is converging to a root, butthe end point C is fixed. Because convergence is judged by the length of the interval[B,C] and because the rate of convergence of the secant rule depends on using valuesfrom the neighborhood of the root, the code monitors the length of the interval.
If fouriterations have not resulted in a reduction by a factor of l/8, the code bisects threetimes. This guarantees that the code will reduce the length of an interval containing aroot by a factor of l/8 in a maximum of seven function evaluations.In summary, if the value D of the secant rule lies outside [B,M] or if the overallreduction in interval length has been unsatisfactory, the code bisects the interval. IfD is too close to B, a minimum change of TOL is used.
Otherwise D is used. Afterdeciding how to compute the next iterate, it is formed explicitly and replaces B. Iff(B) = 0, the code exits. Otherwise, quantities are updated for the next iteration: theold B replaces A. The old C is kept or is replaced by the old B, whichever results inf(B)f(C) < 0.If the code is given normal input [f(x) continuous, f(B)f(C) < 01, then on normal return, either the computed f(B) = 0, or the computed f(B) and f(C) satisfyf(B)f(C) < 0, |f(B)| < |f(C)| and the output values of B and C satisfy (4.13). Inthe latter case there is either a root of f(x) in the interval [B,C] or else one of the endpoints is so close to a root that the sign has been computed incorrectly in the workingprecision.EXERCISES4.9 The algorithm described combining the bisectionmethod with the secant method is very efficient.
Suppose that the initial B and C satisfy |B - C| = 1010,and a root is sought with an absolute error of at most10-5(a) How many function evaluations does the bisection method use?(b) What is the maximum number of function evaluations needed by the combined algorithm?4.3 ROUTINES FOR ZERO FINDINGThe algorithm of the preceding section has been implemented in a routine called Zerodesigned to compute a root of the nonlinear equation F(x) = 0.
A typical invocationof Zero in C++ isflag = Zero(f, b, c, abserr, relerr, residual);in FORTRAN it isCALL ZERO(F, B, C, ABSERR, RELERR, RESIDL, FLAG)and4.3 ROUTINES FOR ZERO FINDING153flag = Zero(f, &b, &c, abserr, relerr, &residual);in C. In FORTRAN F, or f in C and C++, is the name of the function subprogramfor evaluating F(x).
In FORTRAN it must be declared in an EXTERNAL statementin the program that calls ZERO. Normal input consists of a continuous function F(x)and values B and C such that F(B)F(C) < 0. Both B and C are also output quantities,so they must be variables in the calling program. On output it is always the case that|F(B)| < |F(C)|The code attempts to bracket a root between B and C, with B being the betterapproximation, so that the convergence testis satisfied.
It makes no sense to allow RELERR < u, the unit roundoff of the computer used, because this is asking for a more accurate result than the correctly roundedtrue result. To provide a little protection near limiting precision, it is required thatRELERR > 10u. If the desired root should be zero, or very close to zero, a purerelative error test is not appropriate. For this reason it is required that ABSERR > 0.Normal output has either F(B)F(C) < 0 and the convergence test met, or F(B) =0. This is signaled by FLAG = 0. At most 500 evaluations of F are allowed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
