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.