Date: Thu, 02 Feb 95 15:04:49 +0000 Here's a transcription of algorithm #2 from Collected Algorithms from ACM. I have had to change a few details (like substituting the multiply operator by '*', the minor-or-equal operator by '<=', the epsilon symbol by 'Epsilon'). Jose R. Valverde European Bioinformatics Institute Jose.Valverde@ebi.ac.uk (currently) JRValverde@enlil.iib.uk (oldies, still valid) JRamon@mvax.fmed.uam.es --------------------------------------------------------------------------- 2. Rootfinder J. Wegstein National Bureau of Standards, Washington 25, D. C. comment This procedure computes a value of g=x satis- fying the equation x = f(x). The procedure calling statement gives the function, an initial approxi- mation a=/= 0 to the root, and a tolerance parameter Epsilon for determining the number of sig- nificant figures in the solution. This accelerated iteration or secant method is described by the author in _Communications_, June, 1958.; procedure Root(f(), a, Epsilon) =: (g) begin Root: b := a ; c := f(b) ; g := c if (c = a) ; return d := a ; b := c ; e := c Hob: c := f(b) g := (d * x - b * e) / (c - e - b + d) if (abs((g - b) / g) <= Epsilon) ; return e := c ; d := b ; b := g ; go to Hob end ------ CERTIFICATION 2. Rootfinder, J. Wegstein, _Communications_ACM_, February, 1960 Henry C. Thacher, Jr., *Argonne National Labora- tory, Argonne Illinois Rootfinder was coded for the Royal-Precision LGP-30 Com- puter, using an interpretive floating point system with 28 bits of significance. The translation from *Algol* was made by hand. Provision was made to terminate the iteration after ten cycles if convergence had not been secured. The program was tested against the following functions: 1/3 (1) f(x) = (x + 1) (Root = 1.324718) (2) f(x) = tan x -1 (3.a) f(x) = 2*PI*alpha + tan x (alpha = 1, 2, 3, 4) (4.a) f(x) = sinh alpha*x (alpha = -1.2, -0.5, 0.5, 1.2) Selected results were as follows: f(x) alpha epsilon x_k-1 x_k 1 1.3 10^-7, 10^-6 1.3247233 1.3258637 (1) --- ----- 1 1.3 10^-5 1.3247165 (1) -- 2 5 10^-3 -.4674691 -.36021288 (1,2) -------- 2 4 10^-3 +.84880381 +.69496143 (1, 2) -------- 3.1 1 10^-5 7.7252531 -- 3.2 2 10^-5 14.066155 -- 3.3 3 10^-5 20.371026 --- 3.4 4 10^-5 26.665767 ---- (1) No convergence after 10 iterations. Underlined figures are in- correct. (2) For this function f'(0) = 1; so convergence is not to be ex- pected at this root. However, the algorithm did not find any other root. It should be noted that the convergence criterion used fails for a zero root. The provision to terminate after a given number of cycles is therefore essential. Also, double precision is desirable. ------------ * Work supported by the U.S. Atomic Energy Commission. ----------------------- REMARK ON ALGORITHM 2 ROOTFINDER (J. Wegstein, _Communications_ACM_, February, 1960). Henry C. Thacher, Jr., *Argonne National Laboratory Argonne Illinois This algorithm has the convergence factor y - Y (y - Y) f'' k k-2 2 --------- = ------------------------------ + O(y - Y) y - Y 2(f' - 1) + (y - y ) f'' k-1 k-1 k-1 k-2 where Y is the desired root, and the derivatives f' and f'' are evaluated there. Convergence is thus second order, provided that |f''| |y -Y| < 2 |f' - 1|. k-1 The algorithm is, however, somewhat unstable numerically be- cause of the factor f(y ) - f(y ) - y + y in the de- k-1 k-2 k-1 k-2 nominator. Experience has shown that the minimum for epsilon is about one half the precision being used. Provisions to indicate when round- off errors are causing random oscillations of g would be a desirable addition. The criterion used for terminating the iteration renders the algorithm unsuitable for a zero root. A preliminary test for a zero root would be desirable. In addition, the algorithm should include provision for exit after a stated number of iterations. Algorithm 15 appears to offer advantages along these lines. ------------ * Work supported by the U.S. Atomic Energy Commission. REMARKS ON ALGORITHMS 2 and 3 (_Comm._ _ACM_, February 1960), ALGORITHM 15 (_Comm._ _ACM_, August 1960) AND ALGORITHMS 25 AND 26 (_Comm._ACM__, November 1960). J. H. Wilkinson National Physical Laboratory, Teddington. Algorithms 2, 15, 25 and 26 were all concerned with the cal- culation of zeros of arbitrary functions by successive linear or quadratic interpolation. The main limiting factor on the accuracy attainable with such procedures is the condition of the _method_ of evaluating the function in the neighborhood of the zeros. It is this condition which should determine the tolerance which is allowed for the realtive error. With a well-conditioned method of evaluation quite a strict convergence criterion will be met, even when the function has multiple roots. For example, a real quadratic root solver (of a type similar to Algorithm 25) has been used on ACE to find the zeros of triple- diagonal matrices T having t[ij] = a[i], t[i+1,i] = b[i+1], t[i,i+1] = c[i+1]. As an extreme case I took a[1] = a[2] = ... = a5 = 0, a[6] = a[7] = ... = a[10] = 1, a[11] = 2, b[i] = 1, c[i] = 0 so that the func- tions which was being evaluated was x^5 (x-1)^5 (x-2). In spite of the multiplicity of the roots, the answers obtained using float- ing point arithmetic with a 46-bit mantissa had errors no greater than 2^-44. Results of similar accuracy have been obtained for the same problem using linear interpolation in place of the quadratic. This is beacuse the method of evaluation which was used, the two- term recurrence relation for the leading principal minors, _is_a_ _very_well-conditioned_method_of_evaluation_. Knowing this, I was able to set a tolerance of 2^-42 with confidence. If the _same_function_ had been evaluated from its explicit polynomial expansion, then a tolerance of about 2^-7 would have been necessary and the mul- tiple roots would have obtained with very low accuracy. To find the zero roots it is necessary to have an absolute toler- ance for |x[r+1] - x[r]| as well as the relative tolerance condition. It is undesirable that the preliminary detection of a zero root should be necessary. The great power of root finders of this type is that, since we are not saddled with the problem of calculating the derivative, we have great freedom of choice in evaluating the rootfinder so that it finds the zeros of x = f(x) since the true func- tion x - f(x) is arbitrarily separated into two parts. The formal advantage of using this formulation is very slight. Thus, in Certi- fication 2 (June 1960), the calculation of the zeros of x = tan x was attempted. If the function (-x + tan x) were used with a general zero finder then, provided the method of evaluation was, for example x = nPI + y 3 5 y y --- - ----- - ... 3 30 tan x - x = -nPI + ------------------- , cos y the multiple zeros at x = 0 could be found as accurately as any of the others. With a slight modification of common sine and co- sine routines, this could be evaluated as (sin y - y) - y (cos y - 1) -nPI + ----------------------------- 1 + (cos y - 1) and the evaluation is then well-conditioned in the neighborhood of x = 0. As regards the number of iterations needed, the re- striction to 10 (Certification 2) is rather unreasonably small. For example, the direct evaluation of x^60 - 1 is well conditioned, but starting with the values x = 2 and x = 1.5 a considerable number of iterations are needed for Newton's method, starting with x = 2. If the time for evaluating the derivative is about the same as that for evaluating the function (often is much longer), then linear interpolation is usually faster, and quadratic inter- polation much faster, than Newton. In all of the algorithms, including that for Bairstow, it is use- ful to have some criterion which limits the permissible change from one value of the independent variable to the next [1]. This condition is met to some extent in Algorithm 25 by the condition S4, that abs(fprt) < abs(x2 * 10), but there the limitation is placed on the permissible increase in the value of the function from one step to the next. Algorithms 3 and 25 have tolerances on the size of the function and on the size of the remainders r1 and r0 respectively. They are very difficult tolerances to assign since these quantitites may take very small values without our wishing to accept the value of x as a root. In Algorithm 3 (Comm. ACM June 1960) it is useful to return to the original polynomial and to iterate with each of the computed factors. This elimininates the loss of accuracy which may occur if the factors are not found in in- creasing order. This presumably was the case in Certification 3 when the roots of x^5 + 7x^4 + 5x^3 + 6x^2 + 3s + 2 = 0 were attempted. On ACE, however, all roots of thispolynomial were found very accurately and convergence was very fast unsing single- precision, but the roots emerged in increasing order. The reference to _slow_ convergence is puzzling. On ACE, convergence was fast for all the initial approximations to p and q which were tried. When the initial approximations used were such that the real root x = -6.35099 36103 and the spurious zero were found first, the remaining two quadratic factors were of lower accuracy, though this was, of course, rectified by iteration in the original polynomial. When either of the other two factors was found first, then all factors were fully accurate even without iteration in the original polynomial [1]. REFERENCE [1] J. H. Wilkinson. The evaluation of the zeros of ill-conditioned polynomials Parts I and II. _Num. Math._ 1 (1959), 150-180. ---------------------------------------------------------------------- I guess that this is an accurate translation into C: ----- /* * Rootfinder: * J. Wegstein * National Bureau of Standards, Washington 25, D. C. * * This procedure computes a value of g=x satis- * fying the equation x = f(x). The procedure calling * statement gives the function, an initial approxi- * mation a=/= 0 to the root, and a tolerance * parameter Epsilon for determining the number of sig- * nificant figures in the solution. This accelerated * iteration or secant method is described by the * author in _Communications_, June, 1958.; * * * Implemented by: * Jose R. Valverde * European Bioinformatics Institute. * Jose.Valverde@ebi.ac.uk (currently) * * JRValverde@enlil.iib.uk (oldies, still valid) * JRamon@mvax.fmed.uam.es * * Disclaimer: * I made this out of my wits. The European Bioinformatics Institute * does not guarantee or otherwise support or accept any liability * in any manner whatsoever (nor explicit or implicit or ever) with * respect to this software. * * Jose R. Valverde. 2 February 1995 */ double Root(f, a, epsilon) double (*f)(double); double a; double epsilon; { /* Root */ double b, c, d, e,g; b = a; c = (*f)(b); g = c; if (c == a) return (g); d = a; b = c; e = c; for (;;) { /* Hob: */ c = (*f)(b); g = (d * c - b * e) / (c - e - b + d); if (fabs((g - b)/g) <= epsilon) return (g); e = c; d = b; b = g; } }