Construction of rational approximations by means of REDUCE

Abstract

1. In recent years the rational approximations have been widely used to solve physical and computational problems /1,2/. When a real function f(x) is repeatedly calculated on a ≤ × ≤ b, it is reasonable to replace it by a polynomial or rational approximation on [a,b]. For example, if f(x) is a composite combination of elementary and special functions any of which can be calculated by means of the corresponding standard program, the f(x) values are obtainable by these programs. This method, however, involves unjustified losses in the computing time and often provides a too high accuracy inadequate to the problem in question. In this case it is more convenient to use the corresponding approximation. There exist iteration algorithms which ensure the best (in the sense of absolute or relative error) rational approximations based on P.L. Chebyshev theory /2,3/. Unfortunately, these algorithms are cumbersome and do not guarantee the convergence if the choice of the initial approximation is unsuccessful, see ref./2/. The present paper treats simple algorithms (Padé-Chebyshev approximation /1/ and Paszkowski algorithm /4/) providing approximations similar to the best ones with a relatively moderate computer resources required. In this case the calculation of the approximation coefficients reduces to the solution of a system (generally speaking, ill-conditioned) of linear algebraic equations. The errors of the Padé-Chebyshev approximations and the corresponding best approximations are compared in the paper /5/ where one of the methods of computation of the Padé-Chebyshev approximations is described. 2. The Analytic Computations System Reduce is a rather convenient tool of realization of algorithms of the rational approximation construction. This system saves one the trouble of inventing an effective algorithm of approximated-function computation if this function can be given in an analytic form or if the terms in the Taylor series expansion are known or determined analytically by the differential equation. The possibility of using the rational arithmetic (without round-off errors) is essential because the coefficients of rational approximations are not stable with respect to the perturbations of initial data and to the round-off errors. Specifically, the error is minimized which arises in solving the ill-conditioned systems of linear equations and when converting a power series into a series of the Chebyshev polynomials and vice versa. The ALGOL - like input language and the convenient tools of solving the problems of linear algebra ensure the simplicity and compactness of programs. For example, the program of computation of the Padé-Chebyshev coefficients occupies sixty two cards. 3. We compute the Padé-Chebyshev approximations by the standard “cross multiplied scheme”/1/. By means of the change of variable x → [(b-a)x+a+b]/2 the approximation on an arbitrary finite range [a,b] is reduced to the approximation on [-1, 1]. We shall, therefore, restrict ourselves to the case when the function f(x) is to be approximated on the [-1, 1] range by the expression of the form R(x) = a0+a1x+…+an/b0+ b1x+…+bmxm (1) where m, n are given non-negative integer numbers, a0, a1,…, an, b0, b1,…, bm are numerical coefficients to be determined. If the function is specified by a power series ƒ(x)= @@@@ CkXk (2) the corresponding finite sum (the number of terms is determined by the user on the basis of the accuracy required) is converted, using the well-known economizing procedure, to the polynomial @@@@(x) = @@@@ γk Tk (3) where Tk are the k-power Chebyshev polynomials. Solving the system of linear equations 1/2 @@@@ βj (γi+j γli-jl) = 0, i=n+1,…, n+m 1/2 @@@@ βj (γi+j + γli-jl) = αi, i=0,1, …, n (4) we determine the coefficients αi, βi of the rational approximation R(x) = α0+α1T1 + … + αnTn/β0 + β1T1 + … + βm Tm (5) As usual /1,4/ @@@@ dj denotes that the first term in the sum is to be replaced by d0/2. The system (4) is homogeneous and the solution is determined to within the non-zero factor; it is quite natural since the fraction will not change if the numerator and denominator are simultaneously multiplied or divided by a non-zero values. Therefore, the system (4) is completed with the normalization condition, for example β m= 1. Finally, the standard transformation reduces (5) to the form (1). The absolute error of the approximation (1) is of the form δ(x) = &PHgr;(x) / @@@@ bjxj (6) where &PHgr;(x) = @@@@ bjxj ƒ(x) - @@@@ ajxj The above-described algorithm is equivalent to the following procedure: the numerator &Pgr;(x) in (6) is expanded in a series of the Chebyshev polynomials and the first m+n+1 terms are equated to zero. The Paszkowski algorithm (4) leads to the rational approximation R(x) of the form (1) such that the first m+n+1 terms in the expansions of f(x) and R(x) coincide (it will be noted that this approximation does not always exist). The program of this algorithm is similar to that of Padé-Chebyshev algorithm. For even functions the approximation may be sought for in the form R(x) = a0 + a1x2 + … + an(x2)n/b0 + b1x2 + … + bm(x2)m (7) For this case the algorithms involved permit a convenient modification. If f(x) is an odd function it is reasonable to find an approximation of the form (7) to the even function f(x)/x and then to multiply the result by x such that the approximation be of the form R(x) = x a0 + a1x2 + … + an(x2)n/b0 + b1x2 + … + bm(x2)m (8) A large relative error at x=0 is thereby avoided. 4. To estimate the quality of the rational approximation R(x) to the function f(x), the error functions δ(x) = ƒ(x)-R(x) and δ(x)=Δ(x)/ƒ(x) are calculated and the error curves (i.e. the plots of these functions) are built. The absolute error of the approximation coincides with the value of Δ = max \Δ(x)\ and the relative error, with δ = max\δ (x)\. Actually, the error functions are calculated in the finite number of the check points that are uniformly distributed over the range where the function is approximated. To find it out to what extent the approximation R(x) to f(x) differs from the best (in the sense of the absolute and relative error) approximation to f(x) of the same form, one can use the generalized de la Valléé-Poussin theorem /3/. Let us, for example, estimate the similarity of the approximation (1) to the approximation of the same form with the best absolute error. For simplicity, we assume that for the best approximation an≠0, bm≠0 (this condition is usually always fulfilled). Then if a given approximation R(x) is close enough to the best one, the function Δ (x) takes on non- zero values λ1,-λ2,…, (-1)n+mλn+m+2 with alternating signs at successive points x1 < x2 < … < xn+m+2 in the x range by virtue of the generalized Chebyshev theorem /3/. Assume that λ = min {|λ1|,|λ2, …, |λm+n+2|} Then, according to the generalized de la Valléé-Poussin theorem /3/, the following inequality is valid: λ ≤ Δmin ≤ Δ (9) where Δ min is the best absolute error of the approximations. It is clear that Δ coincides with the largest (in the absolute value) extremum of Δ(x) and the least (in the absolute value) extremum of this function (up to a sign) can be used λ. For example, for the Padé-Chebyshev approximation (7) to the function cos(π/4 x) on [-1,1] for m=n=2 the absolute error Δ = =0.68 5 .10-10 and λ = 0.663 .10-10) (2400 check points and the condition of the de la Valléé- Poussin theorem is fulfilled). It is clear that this approximation is rather close to the best one. The errors of the Padé-Chebyshev approximants and the approximants obtained with the help of the Paszkowski algorithm are usually of the same order of magnitude. Consider, as an example, the approximations (1) to the function ex on the range [-1,1] for m=n=3. In this case the exponent is replaced with the truncated Taylor series (2) up to x10/10! inclusive. For the Padé-Chebyshev approximation Δ = 0.33 .10-6, δ = 0.20 . 10-6, and Paszkowski algorithm gives Δ = 0.25 .10-6, δ = 0.26.10-6. The above-described algorithms lead to much lower maximum errors as compared to