3.6 The Q.D. Algorithm and the Root Problem

Then /(z) has the representation

If we are given the formal expansion of /(z),



and we know that f(z) is a meromorphic function which is analytic at the origin, it is natural to wonder if the Pade method is useful for locating the poles and zeros of /(z). Indeed, a vast theory of root solving exists; we will give a simple method which works when the poles of f(z) have distinct moduli, and so also do the zeros. In this case, we may order the poles as with



where h(z) is analytic and zero-free in |z|<min(|uL+1 \,\uM+] |). The special cases where /(z) has either no poles or no zeros cause no problems in principle. However, the condition that the poles and zeros respectively have distinct moduli is important both in principle and in practice. To find the poles of/(z), which means the numerical values of uf, we form [L/M] Pade approximants to /(z). Keeping M fixed, and with increasing order of approximation, we expect that







and, corresponding to (6.8), we define

In fact, this result is a consequence of de Montessus's theorem. If we write







and use the explicit determinantal formula (5.1) for

follows from the usual expression for the product of the roots that


Similarly, by considering M — 1 roots,



We define



For the particular case of M= 1, the appropriate definition is



Then, from (6.6), (6.7), and (6.8),



Recalling the duality theorem (Theorem 1.5.1), the reciprocal of [L/M] is the [M/L] Pade approximant of {/(z)}"1. Let C\M/L) denote the Hankel determinant of the [M/L] Pade approximant of g(z)= {/(z)}"1. With this notation, we find, as in (6.6), that

(6.11)where (6.13) follows from (6.12) by Hadamard's formula (1.6.9). For the particular case of L = 1, the appropriate definition is







By the construction (6.12), the Lth zero of g(z) is

given by



To determine the numerical values of um, vh we need some computational rules.

Product rule.



Proof. From the definitions (6.8), (6.13),



Addition rule.




Proof. Using the definitions (6.8), (6.13) and Sylvester's identity,


Similarly, we find



yields the same answer.


Having established these identities, we display them pictorially in Figure


Figure 1. Pictorial representation of the rhombus rules: (a) product rule, (b) addition rule.




can be con­

structed with the aid of these identities (6.16), (6.17), the initializing values


Jiven by (6.9), and the artificial initializing values



All the quantities




for the first column and



for the first row. We verify that (6.18) gives the correct initializing values by using (6.17) with M=0 together with (6.8) and (6.13). Likewise, (6.19) is verified by using (6.17) with L = 0 together with (6.8) and (6.13). The quantities u(L/M),v{L/M) are usually displayed in the u-v table [Gragg, 1972], shown in Table 1.



Summary. The first two columns and the first row initialize the construc­tion of the u-v table, according to (6.9), (6.18), and (6.19). The other entries are constructed by the Q.D. algorithm, expressed by (6.16) and (6.17). The poles m,,m2,m3,... of /(z) are shown as the limiting values of the "m- columns", and the zeros u,, v2, v3,... of/(z) are shown as limiting values of the "u-rows".

In principle, the rate of convergence of the Q.D. algorithm for poles and zeros is geometric. To see this, we refer to analysis of Section 6.2. From (6.2.11), we find that



and from (6.2.22) we note that the dominant part of (6.20) is given by



where A" is a nonzero constant (independent of L). An inspection of the equations (6.2.11), (6.2.14)-(6.2.22), and (6.2.33) shows that


Equation (6.22) is Hadamard's formula. It holds under the conditions that f(z) is analytic in the disc          except for precisely M poles, counting

multiplicity, in the annulus 0<|r|<R.

Let us assume that the poles at uM_x, uM, and uM+l have distinct moduli:




Substitute (6.22) into (6.8) to show that







The hypothesis (6.23) ensures that O<0<1, and (6.24) shows that the «-columns of the u-v table converge geometrically, in principle.

In practice, the Q.D. algorithm, as expressed in the summary preceding, is unstable. Rounding error accumulates in the «-columns whose limits are

theoretically u2,u^,u4  according to (6.10) and (6.24). To demonstrate


this, we consider





second column

which is obtained from (6.13). Hence, from (6.10),


and we crudely estimate v{L/M) by the formula



We deduce that calculation of the

necessarily involves substantial

rounding error (low relative precision), introduced by cancellation of com­parable quantities in the first M-column. This low relative precision is directly transmitted by the multiplication rule (6.16) to each element of the column (m(L/2), L = 0,1,2,...}. Hence calculation of u2 by the basic Q.D. algorithm is necessarily unstable, and a similar argument extends to the poles m3,m4,.... To remedy this instability, Henrici [1958] proposed the progressive form of the Q.D. algorithm. The order of calculation is changed so as to avoid the necessarily inaccurate arithmetic computations.

Progressive Q.D. algorithm. This algorithm may be used for the computation of the poles of /(z), provided the moduli of these poles are distinct. It is initialized by construction of the coefficients of the power series





are constructed iteratively


from the coefficients using the identity



Hence the quantitities


are constructed, according to (6.14). The other initializing equations (6.9), (6.18), and (6.19) are retained. The progressive form of the Q.D. algorithm requires that the order of calculation using (6.16), (6.17) be as indicated in Table 2. The entries in this schematic section of the u-v table are calculated in numerical order 1,2,3,..., 11,12,13,...,21,22,... as indicated there.

Table 2. Elements of the u-v table showing schematically the direction of calculation with the progressive form of the Q.D. algorithm.



Notice that the progressive form of the Q.D. algorithm may be naturally regarded as a practical method for finding the zeros of a meromorphic




defined by its power series coefficients gh pro-

vided that these zeros of g(z) are known to have distinct moduli. For a treatment of the difficult case where the zeros of g(z) have equal modulus, we refer to Henrici [1974, p. 642], or to Rutishauser [1954, p. 35].

Exercise 1. If f(z) is a polynomial, the sequence (6.9) is ill defined. Explain how the Q.D. algorithm may continue to be used in this case. Exercise 2. Is it a good idea to extrapolate the entries of a M-column of the u-v table using the e-algorithm?