Date: Sun, 12 Feb 95 18:54:30 +0000
Here is a transcription of algorithm #4 from Collected Algorithms
from ACM. I had to modify some characters to make it fit into ASCII,
specially the multiplicative operator has become * and greek letters
are substituted by their equivalent names.
Jose R. Valverde
European Bioinformatics Institute
txomsy@ebi.ac.uk
-------------------------------------------------------------------------
4. BISECTION ROUTINE
S. Gorn
Univeristy of Pennsylvania Computer Center
Philadelphia, Pa.
comment This procedure evaluates a function at the end-points
of a real interval, switching to an error exit (fools
exit) FLSXT if there is no change of sign. Otherwise
it finds a root by iterated bisection and evaluation
at the midpoint, halting if either the value of the
function is less than the free variable epsilon, or two suc-
cessive approximations of the root differ by less
than epsilon1. Epsilon should be chosen of the order of error in
evaluating the function (otherwise time would be
wasted), and epsilon1 of the order of desired accuracy. Epsilon1
must not be less than two units in the last place
carried by the machine or else indefinite cycling will
occur due to roundoff on bisection. Although this
method is of 0 order, and therefore among the slow-
est, it is applicable to any continuous function. The
fact that no differentiability conditions have to be
checked makes it, therefore, an 'old work-horse'
among routines for finding real roots which have
already been isolated. The free varaibles y1 and y2
are (presumably) the end-points of an interval within
which there is an odd number of roots of the real
function F. Alpha is the temporary exit fot the evalua-
tion of F.;
procedure Bisec(y1, y2, epsilon, epsilon1, F(), FLSXT) =: (x)
begin
Bisec: i := 1 ; j := 1 ; k := 1 ; x := y2
alpha: f := F(x) ; if (abs(f) <= epsilon) ; return
go to gamma[1]
First val: i := 2 ; f1 := f ; x := y1 ; go to alpha
Succ val: if (sign(f) = sign(f1)) ; go to delta[j] ; goto etha[k]
Sec val: j := 2 ; k := 2
Midpoint: x := y1 / 2 + y2 / 2 ; go to alpha
Reg delta: y2 := x
Precision: if (abs(y1 - y2) >= epsilon1) ; go to Midpoint
return
Reg etha: y1 := x ; go to Precision
integer (i, j, k)
switch gamma := (First val, Succ val)
switch delta := (FLSXT, Reg delta)
switch etha := (Sec val, Reg etha)
end Bisec
--------------------------------------------------------
CERTIFICATION OF ALGORITHM 4
BISECTION ROUTINE (S. Gorn, _Comm._ACM_,
March 1960)
Patty Jane Rader,* Argonne National Laboratory,
Argonne, Illinois
Bisec was coded for the Royal Precision LGP-30 computer,
using an interpretive floating point system (24.2) with 28 bits of
significance.
The following minor correction was found necessary.
alpha: go to gamma[1] should be go to gamma[i]
After this correction was made, the program ran smoothly for
F(x) = cos x, using the following parameters:
y1 y2 Epsilon Epsilon1 Results
0 1 .001 .001 FLSXT
0 2 .001 .001 1.5703
1.5 2 .001 .001 1.5703
1.55 2 .1 .1 1.5500
1.5 2 .001 .1 1.5625
These combinations test all loops of the program.
-------
* Work supported by the U. S. Atomic Energy Commission.