4.5.3 Viskovatov's Algorithm for the General Case

In order to generalize Viskovatov's algorithm to degenerate cases, one must make an inspired choice for the representation of the particular \L/M] Pade aDDroximant as a continued fraction. We consider

 

 

which should be compared with the general C-fraction (5.12) and Magnus P-fraction (5.16). In (5.17), each p,(z) is a polynomial with ^(0)^=0 and each a; is a positive integer: each p^z) and a, are fully defined in (5.21) and (5.22). As usual, we suppose that we are given the series expansion

 

 

and that we are requested to find the [L/M] Pade approximant of (5.18) with L>M. If L<M, we should remove the factor zk from (5.18), where k is the highest power such that zk factors (5.18), and start with its reciprocal series. We choose to consider an infinite series in (5.18) for convenience of representation; in fact only L + M+ 1 terms of (5.18) are required, and questions of convergence do not arise. Following Werner [1979, 1980], we define the generalized Viskovatov reduction:

We may then define the polynomial Pj(z) of degree strictly less than a by

Initialization. Let

 

 

and

 

 

Iteration. For j=0,1,2,..., proceed recursively. Define a] to be the least integer for which

 

 

 

 

 

Our aim is to find the [/-/w -] Pade approximant of the left-hand side of (5.22). If at >lj + m , then let n=j. In this case, (5.17) is the Pade approxi­mant required, as is proved in Theorem 4.5.1. If

 

 

the [L/M] Pade approximant to the given series (5.18) does not exist, as is proved in Theorem 4.5.1.

Next, we consider the reciprocal series defined by

 

 

so that, formally at least,

 

 

In order to find the [lj/m ] Pade approximant of the left-hand side of (5.25), we seek the [my//y—ay] Pade approximant of the left-hand side of (5.24V Hence we define

 

and

 

 

 

 

 

The requirement that

 

is guaranteed by (5.21), and also

Pade approximant

in view of (5.23V Hence we need to find the

 

of the series

defined by (5.24) and (5.26), and the process

 becomes recursive.

Having defined the iterative process for the construction of (5.17), we must verify that it fulfills its specification. From (5.26), we note that

 

 

and therefore

 

 

where deg{/?(z)} is the degree of p(z). We see that the fraction (5.17) has at most L + M divisors, which occurs if a,— 1 for / = 0,1,2,..., n, and hence the process terminates.

Theorem 4.5.1. The generalized Viskovatov reduction defined in (5.19)— (5.26) is reliable: if the [L/M] Pade approximant exists, the process finds if, if not, (5.23) is satisfied and indicates this fact.

Proof Suppose that the [L/M] approximant to (5.20) exists, yet the algorithm fails. Let the approximant have the Maclaurin expansion

 

 

for

 

 

 

 

 

and

derived by the reduction is identical when applied to

 

itself. Initially, with

we find

 

so that

Since the reduction process (5.25) for calculating the [L/M] approximant is based on the series-expansion coefficients only, the representation (5.17)

 

This division of

 

 

yields 0n(z) as the quotient anc

 

as the remainder. Unless on(z) itself is the [L/ Wm] Padé approximant of

 

/(z), (5.21) ensures that

 

and so we may use the Maclaurin

expansion

 

 

 

 

We sunnose that

 

in (5.31) is reduced sequentially by the

process indicated by (5.19)-(5.26) as far as possible. Let

 

 

fory=0,1k. Hence we reduce the Padé approximant sequentially by

 

 

where we define

 

 

and

 

 

Suppose that the process fails at stage k, namely because (5.23) is satisfied with j=k. From (5.33),

 

 

where lk >ak >lk +mk. Because the left-hand side of (5.37) has degree at most lk, then either ak *Zlk or else T(k\z) = 0. This is incompatible with a failure at stage k. Consequently the reduction process cannot fail if the [L/M] Padé approximant exists.

The proof of the converse result, that R(0)(z) constructed by (5.17)-(5.26 > is in fact the [L/M] Padé approximant, uses the definitions of these equations. From (5.22) and (5.28), we see that

 

 

where

 we define

 

 

and for

 

We define

 

 

 We see directly from (5.40) that

 

 

 

Let us suppose that

 

 

 

It then follows from (5.25) and (5.41) that (5.42) is true with k^k~\. Using (5.39), we have proved that

 

 

by induction. This completes the proof of the theorem.

 

variable

 

is also a factor ol

 

but it is obvious that each common factor of

Hence we see

that we may view the reduction process and Theorem 4.5.1 as a proof, using

In the course of the proof we can see the connection between (5.31). (5.331 and the Euclidean algorithm. The connection is direct using the

the Euclidean algorithm, that P(0\z) and (?<0|(z ) have no common factors. Of course, the proof in Theorem 1.4.3 is more direct.

 

There is also a connection between (5.33) and the recurrence relations (2.14). From (5.33), we find that

 

 

 

 

for

 

We may suppose that the relative normalization

 

is fixed by defining

 

 

 

 

From (5.35), (5.36), and (5.44), it follows that

 

 

and

 

 

which are remarkably similar to (5.2.14) when viewed as a recurrence for j=n—\,n — 2,...,l. In view of (5.34), (5.41), and (5.45), the recurrence defines the numerator and denominator of

 

 

However, the recurrence

(5.48)

 

 

initialized by (5.45), is probably the more convenient

form.

It is now clear how to define a modified Viskovatov algorithm, based on (5.17)-(5.26) and incorporating a generalization of (3.1) for efficiency. In fact, we may as well consider the problem of finding an [L/M] Padé approximant, with L>M, for

 

In (5.49), the coefficients gj°\ dj0) for / = 0.1,2           L + M are given, and

^ 7^0. This problem is solved using bigradients in Section 1.6. We outline a generalized Viskovatov algorithm, emphasizing only the alterations to (5.17)—(5.26) [Bultheel, 1980b],

Initialization. Define l0—L, m0=M.

 

 

 

 

Recurrence. The process is iterative for7 = 0,1,2,.... We construct the polynomial pj{z) of degree at most lJ — mJ such that

(5.50)

with a ■ > lj —mr In fact, we need to construct the coefficients of

(5.511

 

coefficients of the series

 

(5.521

The construction of the coefficients of (5.51) and (5.52) is performed iteratively, according to the algebra implied by the following analysis. Let

 

 

For k = 0,1,...,/—m define

 

 

and

 

 

HencePj(z) is specified by (5.51), and for the coefficients of (5.52) we have

(5.53)

We define a; as the least integer such that (5.53) holds with h\tn 7^0, and hence the coefficients of (5.52) are fully specified.

 

 

If

 

Otherwise, let

or if (5.53) vanishes exactly, the iteration is terminated

the Padé approximant required is nonexistent.

 

and

 

 

 

 

for

 

In this case,

 

 

the /th stage of the reduction

process (5.54) may be iterated. The algorithm is now fully specified.

The generalized Viskovatov algorithm in (5.49)-(5.54) is designed to yield an [L/M] Padé approximant reliably. It follows a path on the Padé table which avoids nonexistent entries in blocks in the table. We conclude that the path that the convergents of the fraction (5.17) follows in the Padé table depends on the values of the coefficients in (5.49) as well as the values of L and M. As a computational algorithm, it can be recommended wherever exact arithmetic is available and the values of a0, a,, ...,«„_, are unambigu­ously defined. The techniques of the generalized Viskovatov algorithm can be adapted to rational interpolation using a generalized Thiele fraction (see Part II, Section 1.1, and Werner [ 1979a, b]) except that the algorithm must incorporate a facility for reordering the interpolation points if necessary [Thacher and Tukey, 1960; Graves-Morris and Hopkins, 1978; Graves- Morris, 1980a], Similar techniques based on the Euclidean algorithm allow Kronecker's algorithm to be generalized to avoid inessential degeneracies (see Section 2.4; Graves-Morris [1980a]).