Page 173 - Compact Numerical Methods For Computers

162 Compact numerical methods for computers
generate the straight line. Thus
b = u-f (u)/f '(u) (13.26)
which gives the zero of the line, is suggested as the next point approximating the
root of f and defines an iteration if b replaces u. The iteration converges very
rapidly except in a variety of delightful cases which have occupied many authors
(see, for instance, Acton 1970, chap 2, Henrici 1964, chap 4) and countless
careless programmers for many hours. The difficulty occurs principally when f´(u)
becomes small. The bisection/False Position combination requires no derivatives
and is thoroughly reliable if carefully implemented. The only situation which may
upset 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 approximation
to Newton’s formula (13.26) by approximating
f '(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 a
pair which straddle at least one root is called the secant algorithm.

Algorithm 18. Root-finding by bisection and False Position
procedure root1d(var 1bound, ubound: real; {range in which
root 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 be
located. Zero is an acceptable value.}
var noroot: boolean {to indicate that interval
on entry may not contain a root since both
function values have the same sign});
{alg18.pas == a root of a function of one variable
Copyright 1988 J.C.Nash
}
var
nbis: integer;
b, fb, flow, fup : real;
notcomp: boolean;
begin
writeln(‘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 displayed
by 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);