Nash - Compact Numerical Methods for Computers (523163), страница 35
Текст из файла (страница 35)
Many interpolantsexist, but the simplest, a straight line, is perhaps the easiest and most reliable touse. This inverse linear interpolation seeks the zero of the liney - f(u) = [f (v)- f(u) ] (b - u) / (v - u)(13.23)which by substitution of y = 0 givesb = [uf(v ) - uf(u)]/[f(v ) - f(u) ] .(13.24)Once f(b) has been evaluated, the interval can be reduced by testing for condition(13.21) as in the bisection process.All the above discussion is very much a matter of common sense. However, incomputer arithmetic the formulae may not give the expected results. For instance,consider the function depicted in figure 13.1(a). Suppose that the magnitude of f(v)is very much greater than that of f(u). Then formula (13.24) will return a value bvery close to u.
In fact, to the precision used, it may happen that b and u areidentical and an iteration based on (13.24) may never converge. For this reason itis suggested that this procedure, known as the method of False Position, becombined with bisection to ensure convergence. In other words, after every fewiterations using the False Position formula (13.24), a step is made using thebisection formula (13.20). In all of this discussion the practical problems ofevaluating the function correctly have been discretely ignored.
The algorithmswill, of course, only work correctly if the function is accurately computed. Acton(1970, chap 1) discusses some of the difficulties involved in computing functions.A final possibility which can be troublesome in practice is that either of theformulae (13.20) or (13.24) may result in a value for b outside the interval [u, v]when applied in finite-precision arithmetic, a clearly impossible situation when thecalculations are exact.
This has already been discussed in §1.2 (p 6). Appropriate tests must be included to detect this, which can only occur when u and v areclose to each other, that is, when the iteration is nearing convergence. The authorhas had several unfortunate experiences with algorithms which, lacking theappropriate tests, continued to iterate for hundreds of steps.Some readers may wonder where the famous Newton’s algorithm has disappeared in this discussion. In a manner similar to False Position, Newton’s methodseeks a root at the zero of the line(13.25)y - f(u) = f´(u)(b - u)where the point (u, f(u)) has been used with f'(u), the derivative at that point, to162Compact numerical methods for computersgenerate the straight line. Thusb = u-f ( u) /f '(u )(13.26)which gives the zero of the line, is suggested as the next point approximating theroot of f and defines an iteration if b replaces u.
The iteration converges veryrapidly except in a variety of delightful cases which have occupied many authors(see, for instance, Acton 1970, chap 2, Henrici 1964, chap 4) and countlesscareless programmers for many hours. The difficulty occurs principally when f´(u)becomes small. The bisection/False Position combination requires no derivativesand is thoroughly reliable if carefully implemented. The only situation which mayupset it is that in which there is a discontinuity in the function in the interval[u, v], when the algorithm may converge to the discontinuity if f( b- )f (b +) < 0,where b - and b + refer to values of b as approached from below and above.It should be noted that the False Position formula (13.24) is an approximationto Newton’s formula (13.26) by approximatingf ' (u) = [f (v) - f (u) / (v - u) .(13.27)The root-finding algorithm based on (13.24) with any two points u, v instead of apair which straddle at least one root is called the secant algorithm.Algorithm 18.
Root-finding by bisection and False Positionprocedure root1d(var 1bound, ubound: real; {range in whichroot is to be found -- refined by procedure}var ifn: integer; {to count function evaluations}tol : real; {the width of the final interval[lbound, ubound] within which a root is to belocated.
Zero is an acceptable value.}var noroot: boolean {to indicate that intervalon entry may not contain a root since bothfunction values have the same sign});{alg18.pas == a root of a function of one variableCopyright 1988 J.C.Nash}varnbis: integer;b, fb, flow, fup : real;notcomp: boolean;beginwriteln(‘alg18.pas -- root of a function of one variable’);{STEP 0 -- partly in the procedure call}notcomp := false; {to set flag or else the ‘known’ root will be displayedby the function routine}ifn := 2; {to initialize the function evaluation count}nbis := 5; {ensure a bisection every 5 function evaluations}fup := fn1d(ubound,notcomp);if notcomp then halt;flow := fn1d(lbound,notcomp);if notcomp then halt; {safety check}writeln(‘f(‘,lbound:8:5,‘)=‘,flow,’f(‘,ubound:8:5,‘)=‘,fup);One-dimensional problems163Algorithm 18.
Root-finding by bisection and False Position (cont.)if fup*flow>0 then noroot := true else noroot := false; {STEP 1}while (not noroot) and ((ubound-lbound)>tol) dobegin {main body of procedure to find root. Note that the test onnoroot only occurs once.}{STEP 9b -- now placed here in Pascal version}if (nbis*((ifn - 2) div nbis) = (ifn - 2)) thenbegin {STEP 10 -- order changed}write(‘Bisect‘);b := lbound + 0.5*(ubound-lbound) {bisection of interval}endelsebegin {STEP 2}write(‘FalseP‘);b := (lbound*fup-ubound*flow)/(fupflow);{to compute false position b}end;{STEP 3 -- begin convergence tests}if b<=lbound thenbeginb := lbound; {to bring argument within interval again}ubound := lbound; {to force convergence, since the functionargument b cannot be outside the interval [lbound, ubound]in exact arithmetic by either false position or bisection}end;if b>=ubound then {STEP 4}beginb := ubound; lbound := ubound {as in STEP 3}end;ifn := ifn+1; {to count function evaluations} {STEP 5}fb := fn1d(b, notcomp);if notcomp then halt; {safety check}write(ifn,’evalns: f(‘,b:16,‘)=‘,fb:10);writeln(‘width interval= ‘,(ubound-lbound): 10);if (ubound-lbound)>tol thenbegin (update of interval)if fb*flow<0.0 then {STEP 6}begin {STEP 7}fup := fb; ubound := b; {since root is in [lbound, b]}endelse {we could check the equal to zero case -- root found,but it will be picked up at STEPs 2, 3, or 4}beginflow := fb; lbound := b; {since root is in [b, ubound] }end; {else}end; {update of interval}end; {while loop}writeln(‘Converged to f(‘,b,‘)=‘,fb);writeln(‘Final interval width =‘,ubound-lbound);end; {alg18.pas = root1d}The algorithm above is able to halt if a point is encountered where the function is not computable.In later minimisation codes we will be able to continue the search for a minimum by presumingsuch points are not local minima.
However, in the present context of one-dimensional root-164Compact numerical methods for computersfinding, we prefer to require the user to provide an interval in which at least one root exists andupon which the function is defined. The driver program DR1618.PAS on the software diskette isin tended to allow users to approximately localise roots of functions using grid search, ,followed bya call to algorithm 18 to refine the position of a suspected root.Example 13.2. A test of root-finding algorithmsIn order to test the algorithm above it is useful to construct a problem which canbe adapted to cause either the bisection or the False Position methods somerelative difficulty.
Therefore, consider finding the root off(b)=z*[tanh(y) +w ](13.28)wherey = s * (b - t )(13.29)b *=t+ln[(l-w )/(l+w )]/(2s ) .(13.30)which has a root atNote that z is used to provide a scale for the function, s scales the abscissa while ttranslates the root to right or left. The position of the root is determined by w towithin the scaling s. Suppose now that s is large, for example s=100.
Then thefunction f(b) will change very rapidly with b near the root, but otherwise will beapproximated very well by(13.31)In fact using the grid search procedure (algorithm 16) we find the values in table13.1 given t=0·5, z=100, s=100 and w=0·99. (The results in the table havebeen computed on a Data General NOVA having 23-bit binary arithmetic.)The consequences of such behaviour can be quite serious for the False Positionalgorithm. This is because the linear approximation used is not valid, and a typicalstep using u = 0, 2, v = 1 givesb=1·00001/200=5·00005E-3.Since the root is near 0·473533, the progress is painfully slow and the methodrequires 143 iterations to converge.
Bisection, on the other hand, converges in 24iterations (nbis=1 in the algorithm above). For nbis=2, 25 iterations arerequired, while for nbis=5, which is the suggested value, 41 iterations areneeded. This may indicate that bisection should be a permanent strategy. However, the function (13.28) can be smoothed considerably by setting w=0·2 ands=1, for which the root is found near 0·297268. In this case the number ofiterations needed is again 24 for nbis=1 (it is a function only of the number ofbits in the machine arithmetic), 6 for nbis=5 and also 6 if nbis is set to a largenumber so no bisections are performed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
