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

 

 

 

 

it

 

and use the explicit determinantal formula (5.1) for

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

(6.6)

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

 

 

 

 

where

 

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

(6.22)

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:

 

 

(6.23)

Substitute (6.22) into (6.8) to show that

 

 

(6.24)

where

 

 

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

 

 

Consequently,

 

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

 

function

 

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?