Quasi-analytical root-finding for non-polynomial functions

A method is presented for the calculation of roots of non-polynomial functions, motivated by the requirement to generate quadrature rules based on non-polynomial orthogonal functions. The approach uses a combination of local Taylor expansions and Sturm’s theorem for roots of a polynomial which together give a means of efficiently generating estimates of zeros which can be polished using Newton’s method. The technique is tested on a number of realistic problems including some chosen to be highly oscillatory and to have large variations in amplitude, both of which features pose particular challenges to root–finding methods.


Introduction
Root-finding is one of the classical problems in numerical analysis, and continues to stimulate research. Our motivation in this paper is to develop methods for the roots of functions with distinct real roots. Our particular application comes from the study of orthogonal functions, both classical and of types with more general orthogonality properties [3, for example], though we emphasize that our method is not limited to these functions and applies to any function which can be represented locally by a Taylor series.
A concrete motive for finding the roots of orthogonal functions is the derivation of quadrature rules, including quadrature rules on the unit circle. The classical orthogonal polynomials-Legendre, Jacobi, Hermite-have well-known quadratures associated with them and a correspondingly extensive body of work on efficient methods for their roots. There are also quadrature rules based on the orthogonality properties of special functions, such as the Bessel-function rules used for the evaluation of Hankel transforms.
While various methods exist which exploit particular special properties of the classical orthogonal functions, the methods of choice for general functions are probably those based on the use of ordinary differential equations for the function. The original form of this algorithm is that of Glaser, Liu, and Rokhlin [8], whose technique uses the second-order ordinary differential equation for the function to generate high-order approximations for stepping between successive roots. Their method uses a Prüfer transform to step from one root to an approximation of the succeeding root, which is then polished using Newton's method applied to a high (thirty or more) order Taylor expansion of the function. The coefficients of the Taylor expansion are efficiently generated using a recursion derived from the differential equation, which allows the function to be evaluated from its Taylor expansion to near machine precision.
The o.d.e. techniques have since been extended using the theory of non-oscillatory phase functions [6] with a recent paper by Bremer [5] presenting a root-finding algorithm which yields roots of the function in O(1) operations per root after a set-up operation which requires very little computation. The algorithm also allows parallel computation of roots, unlike the previous version [8] which required that roots be found sequentially.
Another approach to root-finding for general functions [2] is to represent the functions as combinations of Chebyshev polynomials on sub-intervals of the appropriate domain. Roots, if any, can then be found on each sub-interval as the eigenvalues of the companion matrix of the Chebyshev approximation. This gives good approximations to the roots but does require that the interpolation be generated in advance.
While o.d.e. techniques are probably most convenient for many functions, they do require that a differential equation be available for the stepping between roots, and for generation of the Taylor expansion. When the differential equation is not available, we are forced back on more classical root-finding methods with the corresponding difficulties of root isolation and accurate evaluation. In this paper, we make use of the Taylor series to step from root to root, and as a means of isolating approximations to roots which can then be polished using Newton's method. The approach is not as powerful as the o.d.e.-based algorithm, but is more flexible in that it can be used when a differential equation is either not available or is computationally difficult to use for expansions.
Our root-finding method is based on a combination of two ideas. The first is the use of truncated Taylor series to approximate a function near some point. This is valid over a finite interval whose size depends on the order of terms retained in the Taylor series. Over that interval, the expansion is an accurate representation of the underlying function, to within some user-specified tolerance. The second idea is the use of Sturm sequences to find a root of the Taylor series in the interval. Our argument is that if the Taylor polynomial faithfully represents the function over an interval, then its roots in that interval are good approximations to the roots of the function proper.

Algorithm
The underlying principle of the root-finding algorithm is to locally replace the function f (x) whose distinct real roots are to be found with a polynomial approximation of controlled accuracy. In principle, if the Taylor polynomial (truncated Taylor series) is accurate enough over some range, roots of the Taylor polynomial will be close approximations of roots of the underlying function f (x). Our method is thus composed of two basic tools: the Taylor polynomial as a local representation of the function, and the Sturm sequence which can be used to isolate roots of the Taylor polynomial as initial guesses for the roots of the function proper.

Taylor polynomials
The Taylor polynomial T N ( x, x 0 ) of order N is an approximation of a function f (x) in the vicinity of a point x 0 , with a corresponding estimate of the remainder, We will use this error estimate to fix the range of validity of the Taylor polynomial used in approximating the function whose roots are to be found.
Our method does not depend on any particular means of generating derivatives for the coefficients of the Taylor polynomial, and to some degree a choice of method will depend on the application. In the approach closest to ours [5,8], an o.d.e. is used to generate the required information about the function: our aim here is to provide a method which can be used when an o.d.e. is not readily available. We note that in the general case, high-order derivatives of an analytic function can be stably and accurately found using numerical evaluation of Cauchy integrals [1] though in practice we believe that in most applications there is usually enough local information about the function to provide efficient methods for the computation of derivatives without recourse to Cauchy integrals.
In the applications which motivate this work, the function f (x) is often defined by a sequence of functions {f n (x)} which satisfies a three-term recursion of the form so that one means of generating the derivatives required for the Taylor polynomial is repeated differentiation of the recursion, Although this is the form of recursion used for orthogonal polynomials, for which b n−1 is independent of x, there is no requirement that a n or b n−1 be polynomial, and we impose no such restriction. For example, the Bessel function recursion makes use of rational functions [9, 8.471 and in our numerical tests we will include functions defined by recursions based on algebraic functions.

Sturm sequences
The Sturm sequence is a method for the isolation of roots of a polynomial on the real line, dating to the nineteenth century [13] and which has previously been implemented for use in a polynomial root-finding algorithm based on bisection [10]. Sturm's theorem states that for any polynomial p(x), with simple real roots, there exists a sequence of polynomials p i (x), i = 1, 2, . . . , m. At any point x, the function s(x) can be calculated as the number of sign changes in the sequence p i (x). The number of roots lying between x = a and x = b, with a < b, is then the difference s(a) − s(b). This gives us a means of determining the number of roots of our Taylor polynomial as a first step towards finding roots of f (x). The Sturm sequence for any polynomial p(x) with simple real roots is generated by initializing the sequence with p 1 (x) = p(x), p 2 (x) = p (x). Subsequent terms p k (x) are generated as the remainder of the polynomial division p k−1 (x)/p k (x) with its sign reversed. An implementation of this algorithm in C has been published [10] and consists of the following steps for a polynomial p(x) = N n=0 a n x n : 1. Set p 1 (x) = p(x) and p 2 = p (x); 2. Scale p 2 (x) on the absolute value of its leading coefficient; The output of this algorithm is a sequence of m polynomials which can be evaluated at a point x to check for sign changes as follows: Given these algorithms for evaluation of the Sturm sequence and for counting sign-changes, we can check for and, if necessary, bracket one root in an interval 0 < x < h, by the following bisection method, where α is a desired interval size for the bracketing: The algorithm isolates the root in [0, h] with smallest absolute value, i.e., the 'next' root in our stepping procedure. Replacing h with −h in the first two steps gives a procedure which steps in the negative direction, which allows our root-finding algorithm to be applied starting from an arbitrary point. For example, this allows the method to work forwards and backwards from some starting point to find all roots of a function in a given interval, even when there is no symmetry or other information available about the distribution of roots. This ability to start from an arbitrary point means that the method can also be implemented in parallel to seek multiple roots independently.

Summary
To summarize, we give an algorithm for stepping from an initial point x 0 to the next largest root of f (x), f (x 0 + x) = 0. This step can be repeated to find subsequent roots, and reversed to find roots x 0 − x. The user supplies a starting point, which may be a previously-found root, and an error tolerance ε which is used to control the precision of the Taylor expansion by setting the range over which a root is sought.
To step from point x 0 to the next largest root, using expansions of order N: 1. Localize region containing roots: 2. Isolate one root: (a) When interval containing roots is found, use bisection method of Section 2.2 to isolate one root x; (b) Refine root using Newton's method; (c) Set x 0 to x 0 + x + δ and repeat.
The tolerance ε is user-defined and will depend to some degree on the problem being solved. The small increment δ is required to shift the starting point off the root which has just been found in order to take that root out of the interval covered by the Sturm sequence, to machine precision.

Some notes on implementation
In our numerical examples, we consider functions which have amplitudes comparable to machine precision, a stringent test of a numerical root-finding method. In our experience, the choice of N, ε, and δ affect the performance of the algorithm and we outline how they can be set.
Firstly, the order N of the method interacts with the tolerance ε in setting the step size h. At high order, small tolerance, and large |t N |, h cannot be reliably evaluated, limiting the maximum usable order. We find that priority should be given to ε and the order adjusted rather than vice versa. Given that h N = ε/|t N | a usable criterion is to fix a minimum value of h N , one or two orders of magnitude greater than machine precision say, and use this to fix a minimum value for ε/|t N | which will allow accurate evaluation of h. In our second test case, ε varies over many orders of magnitude at various values of N and the method still gives accurate results.
Secondly, too small a value of δ risks making the method fail on functions of very small amplitude. In our implementation, we shift x 0 by repeatedly incrementing it by δ until |f (x 0 )| > f min , with f min some small tolerance set by the user. Values of these settings are given in the principal tests which we present.
Finally, in principle, there might be functions for which the method of Section 2.3 fails to terminate. We have not encountered this problem in our tests, but for completeness the algorithm should stop and report failure if the number of steps exceeds some user-defined limit.

Numerical examples
We present a number of numerical examples intended to illustrate the behaviour of our algorithm and to show its performance on the kind of problem for which it is intended. The algorithm of Section 2.3 is applied with the Newton refinement step used to evaluate roots to near machine precision in all cases.

Trigonometric plus Gaussian function
Our first test case is used to graphically illustrate the behaviour of the method on a simple function which has closely-spaced roots or near-zeros which are not roots. The function is a combination of a Gaussian and a trigonometric function, which can have a pair of zeros near x = 0 as well as the usual zeros of cos x. The parameter σ controls the spacing of these zeros while the value of a determines whether or not they exist. For a < 0, there are no zeros near the origin, though the function may come very close to zero, while for a > 0 there are two zeros which can be made arbitrarily close. Figure 1 shows the behaviour of the function for different values of σ when there are two roots near the origin. Figure 2 shows the behaviour of the algorithm on the test function for three different values of σ when there are two well-separated roots near zero (a = 1/4) and when the function approaches, but does not touch, the axis (a = −1/128). In each case, the root-finding was initialized at x = −2. Four roots were sought for a = 1/4, and two for a = −1/128, using eighth-order Taylor polynomials with ε = 10 −10 and δ = 10 −12 . The roots are indicated by solid boxes, while intermediate values of the function at the points x 0 used to generate Taylor expansions are shown as circles. The first obvious point to note is that the method does indeed locate the roots accurately. The second point of interest is that the adaptivity is clearly shown by the varying density of the circles on the curve, especially near maxima and minima. This is especially pronounced for σ = 1/4, where the circles are quite dense near the peaks of the function, and separate as the curve descends through the two central roots before bunching again around the minimum.

Highly-oscillatory trigonometric function
The next test case is f (x) = sin(1/x) which is chosen for being highly oscillatory making root-finding a challenging task as x → 0. Roots of the function for x > 0 are x n = 1/(nπ ), n ∈ N, so that the spacing between successive roots |x n − x n+1 | ≈ 1/(πn 2 ) = πx 2 n , x → 0. As a test of the root-finding method, we seek roots starting from points increasingly close to the origin. A basic requirement of the method is that it should work until the distance between successive roots is comparable to machine precision, at which point roots become indistinguishable.
The test involved seeking ten roots moving toward the origin from a starting point x 0 = 10 −m , m = 0, 1, . . . , 7. Roots found were checked to ensure that they corresponded to successive values of n in x n = 1/(nπ ). To allow for the increasing frequency of oscillation as x 0 → 0, the tolerance for setting the step size h was set to ε = x 2 0 /10, while the small increment required to shift the starting point away from a just-found root was set to δ = x 2 0 /1000. The measure of computational effort, shown in Fig. 3, is the total number of function evaluations in finding the ten roots, as a function of m and of the order of the Taylor polynomial. As expected, increasing the order of the method reduces the computational effort considerably-by more than two orders of magnitude for the cases considered here-and starting nearer the origin leads to an increase in effort as the roots become more closely spaced and the interval of validity of the Taylor polynomial becomes smaller so that more function evaluations are needed. In this test case, the tolerance is varied over many orders of magnitude 10 −15 ≤ ε ≤ 10 −3 in order to adapt to the increasing frequency of oscillation of the function, with 4 ≤ N ≤ 10. In all cases, the method continued to give accurate answers with no difficulties arising in selecting a step size h.

Orthogonal SR-functions and para-orthogonal polynomials on the unit circle
Our final set of examples is made up of functions given in terms of a recursion relation. These functions are not polynomials so that many classical root-finding methods are ruled out, nor do we have enough information for use in an o.d.e. method [5,8].
The functions in question have some orthogonality properties defined using a positive measure dφ(x) on the interval (−1, 1). A sequence of such functions {W n (x)} ∞ n=0 has been studied previously [3] using the orthogonality properties for n, m = 0, 1, 2, . . ., with ρ 2m ρ 2m+1 = 0, δ n,m = 0,if n = m and δ n,m = 1,if n = m. It has been proved that the sequence of orthogonal functions {W n (x)} ∞ n=0 satisfies the following three-term recurrence relation We refer to the functions W n (x) which satisfy (7) as SR-functions (square-root), to emphasize that they are not polynomial, though they do reduce to polynomials in the special case c m+1 = 0.
From the recurrence relation (7), one can see that an SR-function can be written as where A n (x) ∈ P n and B n−1 (x) ∈ P n−1 satisfy A n (−x) = (−1) n A n (x) and B n−1 (−x) = (−1) n−1 B n−1 (x),and P n denotes the space of the polynomials of degree at most n.
Our interest in these functions comes from the fact that in the interval (−1, 1) an SR-function, W n (x), which satisfies the orthogonality properties (6) has exactly n simple zeros [3].
Furthermore, the orthogonal SR-functions have an interesting connection with orthogonal polynomials on the unit circle. Let dμ(z) be a positive measure on the unit circle T = {z = e iθ : 0 ≤ θ ≤ 2π }. Then, we call the sequence of polynomials {S n (z)} ∞ n=0 such that It is well known that [12] all zeros of these polynomials lie inside the unit disk and that they satisfy the relations S n (z) = zS n−1 (z) − α n−1 S * n−1 (z), n ≥ 1, where α n−1 = −S n (0) and the reciprocal polynomial of S n (z) is S * n (z) = z n S n (1/z). Para-orthogonal polynomials on the unit circle (POPUC) associated with the sequence of OPUC {S n (z)} ∞ n=0 can be given by zS n−1 (z) − τ n S * n−1 (z), n ≥ 1, where {τ n } ∞ n=1 is any sequence of complex numbers such that |τ n | = 1. Unlike the OPUC, the POPUC has all its n zeros simple and lying on the unit circle T.
We must emphasize that, from relation (9), the zeros of the POPUC (z − 1)R n (z) which lie on the unit circle and are used in quadrature rules on the unit circle [11] can easily be computed from z j = e iθ j , θ j = 2 arccos(x j ), where x j , j = 1, 2, · · · , n,are the zeros of the orthogonal SR-function W n (x).
Note that when η = 0, the orthogonal SR-functions W n (x) reduce to the monic ultraspherical polynomials C (λ−1/2) n (x). Varying λ and η gives quite different behaviour to the family of functions and can lead to quite testing conditions for a root-finding algorithm. Figure 4 shows the roots of W 24 (x) with λ = 0.25, η = 0.9. Roots were found working downwards from x 0 = 1 − 10 −7 , with δ = 10 −7 . Other than the square-root behaviour at the end points, this is a reasonably well-behaved function without a large variation in amplitude. The roots are well-separated as is clear from Fig. 4. Figure 5 is another matter: the function is of quite large amplitude for x −1/2 and is then exponentially small, with oscillations of increasing frequency as x → 1, where the function has a square root singularity. Nonetheless, our method finds the roots in −1 ≤ x ≤ 1 accurately, as can be seen from the lower plot of Fig. 5 which gives the detail of the function and roots for 0.9 ≤ x ≤ 1. The roots are correctly identified despite the small amplitude of the function (note the scaling on the vertical axis) and the high frequency oscillations.
Finally, Fig. 6 shows the computational effort required to find the roots in Fig. 5, as a function of the order N of the method. The reduction in effort, shown as the number of function evaluations required, is quite marked, falling by more than two orders of magnitude as N increases from 4 to 9, though the improvement is slower as N increases beyond this point.
Furthermore, the orthogonal SR-functions W n (x) satisfy the recurrence relation (7) with the same coefficients {c n } ∞ n=1 and {d n+1 } ∞ n=1 given in (12). Figure 7 shows results of root-finding for n = 29, κ = 0.8. Roots were found using eighth-order Taylor polynomials, starting from x 0 = −0.9999, with δ = 10 −7 and ε = 10 −12 . From the upper plot, it is clear that the roots have been found successfully. The lower plot shows that, importantly, the algorithm has correctly detected and separated two roots which are close but distinct, near x = √ 2/2. This is an essential feature of any root-finding method, especially one used for orthogonal or para-orthogonal functions which form the basis of a quadrature rule.
Finally, to complete the calculation, Fig. 8 shows the roots of the POPUC (z − 1)R n (z) for n = 29, associated with the measure dμ (κ) (z) given by (11) with κ = 0.8, computed from the real roots of W 29 (x), with the addition of the root z = 1.

Conclusions
We have presented a method which is well-suited to the problem of finding simple real roots of non-polynomial functions when a differential equation is not available for use in alternative algorithms. We have demonstrated that the method works on some realistic problems, including some deliberately chosen to pose a challenge to root-finding methods. The technique is accurate and efficient, and can be implemented using standard coding tools. We note finally that the approach is easily modified to extract roots of first or other derivatives, such as function extrema, which are important in many applications, and that it is easily implemented on parallel systems.