2.4 Numerical Calculation of Padé Approximants

Quite probably, the first problem to be tackled by anyone interested in Padé approximants involves the calculation of values of a set of diagonal approximants using the coefficients of the Maclaurin series as data. Explicit calculation is only feasible by hand for the lowest-order approximants, and so a good numerical algorithm becomes an early requirement. In this section we discuss the methods available and the qualities which are required of a "good" method.


Here 6y=:0 as

Thus ôfi/fi is of order -q, and 6/1 /A, which by (2.7)

We start by considering the direct method of solution (see Section 2.1), and defer the detailed discussion of the various distinct objectives and the other methods available until this discussion can be put in context. Numeri­cal solution of the linear equations are real or complex numbers known to a certain is the core of the problem. We assume that the data coefficients in this way and completes the calculation of an

accuracy. If /<0 in (1), we take c\ =0. The simplest approach to (1) which can be recommended is to solve the linear system using Gaussian elimination a (real) point is given in the appendix. The function pade(- ■) is written for and may readily use with real, single-precision coefficients be adapted to suit other cases.

An obvious question to put is whether Gaussian elimination with full pivoting is the best method of solving (4.1). If we decide to use the direct method of calculating Padé approximants, most of the work is entailed in ror our the actual solution of (4.1) for the coefficients urposes, we need a method which will recognize whether the system of

Unless the data equations is degenerate, which occurs if coefficients are integers (or rational fractions with an exact binary representation) for which symbolic methods are suitable [Creddes, 19791, it is likely that rounding errors will prevent our computer program calculation which we obtain are generated by rounding error and are likely to be misleading. We stress that this sort of case does arise: for example, if f(z) is a geometric series, all its [L/2] Padé approximants are degenerate. Our numerical algorithm must be able to recognize similar cases in which a lower order of Padé approximation is quite possibly what is required. In short, one of our requirements is an accurate method of solving the system (4.1); a method which is more accurate than one might at first sight expect is needed to discriminate against exactly degenerate equations. The accuracy requirement alone is achieved using the Gaussian elimination method with double-precision arithmetic, as in the program in the appendix of Part I, or even higher precision if necessary, but this head-on attack is not quite the best approach.

We wish to ascertain not only whether the data coefficients correspond to an exactly degenerate problem, but also whether the given data constitute an ill-posed problem or correspond within error to a degenerate Padé approximation problem. We wish to know if the system can become degenerate when the coefficients are varied within their given accuracy limits. The method consisting of equilibration, partial pivoting for genera­tion of an LU decomposition and subsequent iterative refinement of the solution is especially useful in near-degenerate cases. We refer to Wilkinson [1963,1965] for details of this numerical method of solving general linear equations and its error analysis. Basically, this sophisticated method consists of finding the LU decomposition of the coefficient matrix in (4.1), which is used to generate a sequence of approximate solution vectors. The iterative refinement converges geometrically for stable well-posed problems, and also for the ill-posed problems which occur so frequently in Padé approximation. If the iteration does not converge properly, we conclude that the solution is generated by rounding error in the data coefficients, which almost certainly corresponds to a degenerate approximant. The numerical convergence char­acteristics of the iterative refinement permit better discrimination between the ill-conditioned but nondegenerate coefficient matrices and singular coefficient matrices corresponding to degenerate approximants. This method is implemented in the (Mark 6) N.A.G. library algorithm.

Whatever method of solution of (4.1) is chosen, the complete algorithm must contain an error exit for use when the requested approximant is degenerate according to agreed criteria. If the algorithm includes pointwise evaluation of an approximant, the algorithm should contain an error exit in case the approximant is being evaluated at a pole.

We have emphasized that a good numerical algorithm should contain a fairly sophisticated degeneracy test. If the Padé approximant corresponding to given data does not exist, we expect good algorithms to recognize this situation and to detect the degeneracy. We call such algorithms reliable. We might also expect that small variations of the data coefficients within their accuracy limits should not significantly affect the computed value of the Padé approximants. However, we expect to face ill-posed problems in Padé approximation, and so, more realistically, we expect the vanishingly small variation of the data coefficients to correspond to vanishingly small varia­tion in the computed numerator and denominator coefficients. This reflects the idea of mathematical continuity, and methods having this property are called stable. With ill-posed problems, we recognize that there will be magnification of rounding errors in the solution. Beyond these considera­tions, we expect our algorithm to be efficient, but efficiency is not as important as reliability and stability. Experienced Padé approximators are wary of degenerate and near-degenerate cases; they prefer to use more sophisticated approximation methods (see Part II, Chapter 1) rather than to resort to high-order approximants. Furthermore, the actual calculation of the Padé approximants is usually the least time consuming part of a complete calculation, in which it is the computation of the coefficients to high accuracy which is expensive in computer time. We note that solution of (4.1) by elimination requires 0{\M3) operations, unless iterative refinement is necessary. We consider methods of lower order, which are more efficient when M is large, later in the section and discuss the penalties incurred by increased efficiency.

Our concern for accuracy of the numerical algorithm stems primarily from experience. In the previous section, we mentioned the empirical rule of keeping M extra decimal places of precision for computation of an [L/M] approximant. Of course, one cannot say a priori that M, M— 1, 2M, or any other number of extra decimal places are generally required. We simply sound a clear warning about the minimum number of guarding figures likely to be required, and suggest that any calculations of [L/M] approximantsperformed with less than M guarding figures be looked at critically. From a numerical point of view, the Pade approximant derives its capacity to extrapolate certain power series beyond their circle of convergence from using the information contained in the tails of the decimal expansion of the data: naturally, accurate data are of paramount importance. We understand some of these ideas in terms of the following example. We consider a Stieltjes series (see Chapter 5; this and other material of the section may be more suitable for a second reading of this work) in which the coefficients are given by

c = f"uJd<j>(u), j=0,1,2,...       (4.2)


where «>0 and <p(u) is bounded and nondecreasing. (Our sign convention here differs from (5.1.2) in that/(— z) is a Stieltjes function.) A special case of (4.2) is given by



where «>0 and <p(u) is bounded and nondecreasing. (Our sign convention here differs from (5.1.2) in that/(— z) is a Stieltjes function.) A special case of (4.2) is given by



CJ=J*"JdU=jTT- V = 0J,2,...  (4.3)


For Stieltjes series, various sequences of Pade approximants are known to converge systematically, and Stieltjes functions form the major class of functions for which we have a complete convergence theory. We see that the equations (4.1) with coefficients given by (4.3) are



for L>M — 1. The coefficient matrix of (4.4) is a Hilbert segment which is notorious for its ill-conditioning. In fact, we have encountered ill- conditioning in an ideal case which is in no way degenerate. In this context, it is interesting to consider the condition number of the Hankel matrix


in which the elements of HM are defined by (4.2). The condition number of HM is defined by



A matrix with a large condition number corresponds to an ill-posed problem, and very crudely, one may expect that \ogwn(HM) is the number of significant figures lost in the solution of (1) with rounded data. In the general case defined by (4.2) and (4.5), Taylor [1978] has shown that provided defined by (4.3) and (4.5).


These results strongly support our rule of thumb that approximately M decimal places of accuracy are lost in the calculation of an [L/M] ap­proximant. In reality, this means that approximately M extra decimal places any further significant loss of accuracy to occur in the calculation of the numerator coefficients using (2.1b).

Hitherto in this section, we have tacitly assumed that our problem is the calculation of the coefficients of a particular Padé approximant with a view to calculating the approximant at a prespecified value of z. In practice, it is much more likely that our task consists of the tabulation of a particular [L/M] approximant at a number of values of z, or else that it consists of pointwise evaluation of a whole sequence of Padé approximants at a prespecified value of z. Which is the best algorithm to choose for Padé approximation depends on the problem at hand. It is normal to distinguish between the coefficient problem and the value problem. If tabulation of values of a particular approximant is required, it is probably best to calculate the coefficients first, and then to evaluate the approximants. If pointwise evaluation of a whole sequence is required, it may be better to consider methods such as using the e-algorithm or the Q.D. algorithm in which a whole sequence is calculated recursively.


It is unusual tor of precision are required of the data coefficient than is expected of the solution coefficients


The direct method of calculating Padé approximants is stable and reli­able, but may not be the most efficient such method. There are several "0(aM2)" methods, where a is typically 4 or 6, which have computational efficiency gained at the expense of reliability. This means that they are suitable for use in contents where the existence and nondegeneracy of the

approximants is not in question. Possibly, any of them may be developed into a reliable method at some future time. Our list of 0(aM2) methods, each with its own individual merits, for the coefficient problem is

Kronecker's algorithm for the numerator and denominator coeffi­cients.

The Q.D. algorithm for the corresponding continued-fraction repre­sentation.

Baker's algorithm for the numerator and denominator coefficients.

Watson's algorithm for the numerator and denominator coefficients.

Toeplitz and Hankel matrix methods for the denominator coeffi­cients.

Kronecker's algorithm [1881]. This is discussed in Part II, Section 1.1. The application to Pade approximation is the case in which all the interpo­lation points are confluent:



The coefficients of the Newton interpolating polynomial in this case are simply the Maclaurin series coefficients Cq, C|,..., m- Kronecker's algo­rithm is essentially based on exploitation of the (***) identity (3.5.12) to develop an antidiagonal sequence in the Padé table. Consequently, the algorithm is well defined only if all elements of the antidiagonal sequence exist and are nondegenerate. However, in every degenerate case, Kronecker's algorithm can be extended using the Euclidean algorithm [Claessens (thesis), 1976; McEliece and Shearer, 1978; Cordellier, 1979b; Graves-Morris, 1979].

The Q. D. algorithm. This algorithm develops a continued-fraction repre­sentation for a sequence of Padé approximants, and is fully discussed in Section 4.4. We require that all the elements of the corresponding fraction be finite and nonzero, which is equivalent to requiring that the convergents of the continued fraction be nondegenerate.


Baker's algorithm [1973a], We construct the sequence of numerators {•>),■(z), / = 0,1,2  } and denominators {#,(z), / = 0,1,2,...} of Padé ap­proximants in a staircase sequence given by




The coefficients of tj,(z), <9,(z) are related recursively using formulas based on the (*£) and (J*) identities (3.5.11). Essentially, this algorithm is a fast method for calculating an antidiagonal staircase sequence, and we refer to EPA, p. 77, and Pindor [1976] for details.Watson's algorithm [19731. We construct the sequence of numerators proximants in a staircase sequence given by



The coefficients of f)j(z),6i(z) are calculated recursively using formulae based on the (**) and (*|) identities (3.5.11) and the accuracy-through-order condition. This is a fast method for calculating staircase sequences parallel to the diagonal. We refer to EPA p. 79, Watson [1973], and Pade [1894]. Gragg's algorithm [1974] for the denominator polynomials is somewhat similar, but does not exploit the accuracy-through-order condition.



Toeplitz and Hankel matrix methods. The equations (4.1) have a Hankel coefficient matrix. The structure of an NXN Hankel matrix H is such that its (/, j) element is given by








The coefficient matrix of (4.6) is a Toeplitz matrix T. The structure of an NXN Toeplitz matrix T is such that its (i, j) element is given by



Special methods exist which exploit the structure of Toeplitz and Hankel matrices to enable them to be solved in nondegenerate cases by an 0(aN2) method. Such methods do not use pivoting, and the calculation is an iterative scheme based on an embedding principle. The Padé numerator coefficients are calculated by back substitution in (2.1b). We refer to Bareiss [1969], Trench [1964,1965], Rissanen [1973], Zohar [1974], Bultheel [1979], Graves-Morris [1979], Brent et al. [1980], and Bultheel and Wuytack [1980] for details.

We refer to Brezinski [1976] and Bultheel [1979, 1980a] for details of reliable calculation of a paradiagonal sequence of Padé approximants. Werefer to Luke [1980] for a note on the accuracy of the value problem for small values of |z|, and to Pindor [1979a] for some interesting conjectures on the value problem. An up-to-date bibliography is given by Wuytack [1979].

Exercises 1. Find the [3/3] Padé approximant of Euler's series (5.5.9) using each of the methods of this section.

Exercise 2. If multiplication and division are regarded as the significant arithmetic operations, show that Baker's algorithm is an 0(4M2) method, and that Kronecker's algorithm is an 0(4M2) method for MCL.