Rigorous Uniform Approximation of D-finite Functions Using Chebyshev Expansions
A wide range of numerical methods exists for computing polynomial approximations of solutions of ordinary differential equations based on Chebyshev series expansions or Chebyshev interpolation polynomials. We consider the application of such methods in the context of rigorous computing (where we need guarantees on the accuracy of the result), and from the complexity point of view.
It is well-known that the order- truncation of the Chebyshev expansion of a function over a given interval is a near-best uniform polynomial approximation of the function on that interval. In the case of solutions of linear differential equations with polynomial coefficients, the coefficients of the expansions obey linear recurrence relations with polynomial coefficients. Unfortunately, these recurrences do not lend themselves to a direct recursive computation of the coefficients, owing among other things to a lack of initial conditions.
We show how they can nevertheless be used, as part of a validated process, to compute good uniform approximations of D-finite functions together with rigorous error bounds, and we study the complexity of the resulting algorithms. Our approach is based on a new view of a classical numerical method going back to Clenshaw, combined with a functional enclosure method.
Key words and phrases:Rigorous computing, computer algebra, complexity, D-finite functions, recurrence relation, Chebyshev series, Clenshaw method, Miller algorithm, asymptotics, functional enclosure
2000 Mathematics Subject Classification:68W30, 33C45, 65L05, 65G20,
- 1 Introduction
- 2 Chebyshev Expansions of D-finite Functions
- 3 Convergent and Divergent Solutions
- 4 Computing Approximation Polynomials
- 5 Chebyshev Expansions of Rational Functions
- 6 Validation
Many of the special functions commonly used in areas such as mathematical physics are so-called D-finite functions, that is, solutions of linear ordinary differential equations (LODE) with polynomial coefficients . This property allows for a uniform theoretic and algorithmic treatment of these functions, an idea that was recognized long ago in Numerical Analysis [34, p. 464], and more recently found many applications in the context of Symbolic Computation [67, 54, 32]. The present article is devoted to the following problem.
Let be a D-finite function specified by a linear differential equation with polynomial coefficients and initial conditions. Let . Given and , find the coefficients of a polynomial written on the Chebyshev basis , together with a “small” bound such that for all .
Approximations over other real or complex segments (written on the Chebyshev basis adapted to the segment) are reduced to approximations on by means of an affine change of variables, which preserves D-finiteness.
A first motivation for studying this problem comes from repeated evaluations. Computations with mathematical functions often require the ability to evaluate a given function at many points lying on an interval, usually with moderate precision. Examples include plotting, numerical integration, and interpolation. A standard approach to address this need resorts to polynomial approximations of . We deem it useful to support working with arbitrary D-finite functions in a computer algebra system. Hence, it makes sense to ask for good uniform polynomial approximations of these functions on intervals. Rigorous error bounds are necessary in order for the whole computation to yield a rigorous result.
Besides easy numerical evaluation, polynomial approximations provide a convenient representation of continuous functions on which comprehensive arithmetics including addition, multiplication, composition and integration may be defined. Compared to the exact representation of D-finite functions by differential equations, the representation by polynomial approximations is only approximate, but it applies to a wider class of functions and operations on these functions. When we are working over an interval, it is natural for a variety of reasons to write the polynomials on the Chebyshev basis rather than the monomial basis. In particular, the truncations that occur during most arithmetic operations then maintain good uniform approximation properties. Trefethen et al.’s Chebfun [58, 17] is a popular numerical computation system based on this idea.
In a yet more general setting, Epstein, Miranker and Rivlin developed a so-called “ultra-arithmetic” for functions that parallels floating-point arithmetic for real numbers [20, 21, 30]. Various generalized Fourier series, including Chebyshev series, play the role of floating-point numbers. Ultra-arithmetic also comprises a function space counterpart of interval arithmetic, based on truncated series with interval coefficients and rigorous remainder bounds. This line of approach was revived with the introduction of “ChebModels” in recent work by Brisebarre and Joldeș . Part of the motivation for Problem 1.1 is to allow one to use arbitrary D-finite functions as “base functions” at the leaves of expression trees to be evaluated using ChebModels.
Finally, perhaps the main appeal of ultra-arithmetic and related techniques is the ability to solve functional equations rigorously using enclosure methods [44, 30, 37, 45, 60]. LODE with polynomial coefficients are among the simplest equations to which these tools apply. A third goal of this article is to begin the study of the complexity of validated enclosure methods, from a computer algebra point of view, using this simple family of problems as a prototype.
To specify the D-finite function , we fix a linear homogeneous differential equation of order with polynomial coefficients
We also write . Up to a change of variable, we assume that we are seeking a polynomial approximation of a solution of (1.1) over the interval . The uniform norm on this interval is denoted by . We also assume that for , so that (by Cauchy’s existence theorem for complex LODE) all solutions of (1.1) are analytic on . Besides the operator , we are given initial values
Many of the results actually extend to the case of boundary conditions, since we can compute a whole basis of solutions of (1.1) and reduce boundary value problems to initial value problems by linear algebra. Also note that the case of initial values given outside the domain of expansion may be reduced to our setting using numerical analytic continuation.
|coefficients of the operator ;||p. 1.1|
|coefficients of the operator||p. ii|
|sequences in with exponential decrease||p. 2.1|
|differentiation operator,||p. 1.1|
|differential operator,||p. 1.1|
|initial values,||p. 1.2|
|Chebyshev recurrence operator||p. 2.1, p. 2.5|
|degree- minimax polynomial approximation of||p. 2.5|
|truncated Chebyshev expansion operator||p. 2.3|
|order of||p. 1.1|
|shift operator,||p. 2.2|
|(usually) half-order of||p. ii, p. 3.4|
|singularities of , shifted by||p. 4.1|
|Chebyshev polynomials of the first kind||p. 2.1|
|(usually) unknown function,||p. 1.1, p. 5.1|
|approximation of computed by Algorithm 4.2||p. 4.2|
|inverse Joukowski transform of a function||p. 2.4|
|integer interval,||p. 1|
Table 1 summarizes for quick reference the notation used throughout this article. Notations related to Chebyshev expansions are detailed in Section 2.1 below. Notations from Theorem 2.1 are also repeatedly used in the subsequent discussion. Double brackets denote integer intervals .
Unless otherwise noted, we assume for simplicity that all computations are carried out in exact (rational) arithmetic. The rigor of the computation is unaffected if exact arithmetic is replaced by floating-point arithmetic in Algorithm 4.2 and by interval arithmetic in Algorithm 6.5. (In the case of Algorithm 5.6, switching to interval arithmetic requires some adjustments.) However, we do not analyze the effect of rounding errors on the quality of the approximation polynomial and error bound from Problem 1.1 when the computations are done in floating-point arithmetic. In simple cases at least, we expect that Algorithm 4.2 exhibits comparable stability to similar methods based on backward recurrence . Our experiments show a satisfactory numerical behaviour.
To account for this variability in the underlying arithmetic, we assess the complexity of the algorithms in the arithmetic model. In other words, we only count field operations in , while neglecting both the size of their operands and the cost of accessory control operations. The choice of the arithmetic complexity model for an algorithm involving multiple precision numerical computations may come as surprising. Observe however that all arithmetic operations on both rational and floating-point numbers of bit size bounded by may be done in time (see for instance Brent and Zimmermann’s book  for detail). In general, the maximum bit size of the numbers we manipulate is roughly the same as that of the coefficients of , which may be checked to be when represented as rational numbers, so that the bit complexity of the algorithm is actually almost linear in the total bit size of the output.
1.3. Summary of Results
As we will see in the next section, truncated Chebyshev expansions of analytic functions provide very good approximations of these functions over straight line segments. In the case of D-finite functions, their coefficients are known to satisfy linear recurrences. But computing Chebyshev series based on these recurrences is not entirely straightforward.
Roughly speaking, the conclusion of the present article is that these recurrences can nevertheless be used to solve Problem 1.1 efficiently for arbitrary D-finite functions. The techniques we use (backward recurrence and enclosure of solutions of fixed-points equations in function spaces) date back to the 1950s–1960s. The originality of this work is that we insist on providing algorithms that apply to a well-defined class of functions (as opposed to methods to be adapted to each specific example), and focus on controlling the computational complexity of these algorithms.
Our algorithm proceeds in two stages. We first compute a candidate approximation polynomial, based on the Chebyshev expansion of the function . No attempt is made to control the errors rigorously at this point. We then validate the output using an enclosure method.
The main results of this article are Theorems 4.4 (p. 4.4) and 6.6 (p. 6.6), stating respectively that each of these two steps can be performed in linear arithmetic complexity with respect to natural parameters, and estimating the quality of the results they return. Theorem 4.4 is based on a description of the solution space of the recurrence on Chebyshev coefficients that is more complete than what we could find in the literature and may be of independent interest.
This article is organized as follows. In Section 2, we review properties of Chebyshev series of D-finite functions and then study the recurrence relations satisfied by the coefficients of these series, whose use is key to the linear time complexity. Section 3 provides results on the asymptotics of solutions of these recurrences that will be critical for the computation of the coefficients. The actual algorithm for this task, described in Section 4, reminds of Fox and Parker’s variant [23, Chap. 5] of Clenshaw’s algorithm . A short description of a prototype implementation and several examples follow.
The part dedicated to the validation step starts in Section 5 with a study of Chebyshev series expansions of rational functions. Most importantly, we state remainder bounds that are then used in Section 6, along with an enclosure method for differential equations, to validate the output of the first stage and obtain the bound . We conclude with examples of error bounds obtained using our implementation of the validation algorithm, and some open questions.
2. Chebyshev Expansions of D-finite Functions
2.1. Chebyshev Series
Recall that the Chebyshev polynomials of the first kind are polynomials defined for all by the relation . They satisfy for all . The family is a sequence of orthogonal polynomials over with respect to the weight function , and hence a Hilbert basis of the space . (We refer the reader to books such as Rivlin’s  or Mason and Handscomb’s  for proofs of the results collected in this section.)
Expansions of functions on this basis are known as Chebyshev series. Instead of the more common
we write Chebyshev series as
This choice makes the link between Chebyshev and Laurent expansions as well as the action of recurrence operators on the (both discussed below) more transparent. The Chebyshev coefficients of the expansion of a function are given by
Now assume that is a solution of Equation (1.1). As such, it may be analytically continued to any domain that does not contain any singular point of the equation. Let
be the largest elliptic domain with foci in with this property. Since the singular points are in finite number, we have . The coefficients then satisfy for all ; and the Chebyshev expansion (2.2) of converges uniformly to on [39, Theorem 5.16]. Letting and , it is not hard to see that the are also the coefficients of the (doubly infinite) Laurent expansion of the function around the unit circle. The transformation sending to is known as the inverse Joukowski transform. It maps the elliptic disk to the annulus
The formula translates into . The coefficients are also related to those of the Fourier cosine expansion of .
Let be the vector space of doubly infinite sequences such that
The sequence of Chebyshev coefficients of a function that is analytic on some complex neighborhood of belongs to . Conversely, for all , the function series converges uniformly on (some neighborhood of) to an analytic function .
Truncated Chebyshev series are near-minimax approximations: indeed, they satisfy [59, Theorem 16.1]
where is the polynomial of degree at most that minimizes .
Even though itself can be computed to arbitrary precision using the Remez algorithm [12, Chap. 3], Equation (2.5) shows that we do not lose much by replacing it by . Moreover, tighter approximations are typically hard to validate without resorting to intermediate approximations of higher degree . The need for such intermediate approximations is actually part of the motivation that led to the present work. There exist a variety of other near-minimax approximations with nice analytical properties, e.g., Chebyshev interpolation polynomials. Our choice of truncated Chebyshev expansions is based primarily on the existence of a recurrence relation on the coefficients when is a D-finite function.
2.2. The Chebyshev Recurrence Relation
The polynomials satisfy the three-term recurrence
as well as the mixed differential-difference relation
which translates into the integration formula where . From these equalities follows the key ingredient of the approach developed in this article, namely that the Chebyshev coefficients of a D-finite function obey a linear recurrence with polynomial coefficients. This fact was observed by Fox and Parker [22, 23] in special cases and later proved in general by Paszkowski . Properties of this recurrence and generalizations to other orthogonal polynomial bases were explored in a series of papers by Lewanowicz starting 1976 (see in particular [35, 36]). The automatic determination of this recurrence in a symbolic computation system was first studied by Geddes .
The following theorem summarizes results regarding this recurrence, extracted from existing work [49, 35, 36, 52, 4] and extended to fit our purposes. Here and in the sequel, we denote by the skew Laurent polynomial ring over in the indeterminate , subject to the commutation rules
Likewise, is the subring of noncommutative Laurent polynomials in themselves with polynomial coefficients. The elements of identify naturally with linear recurrence operators through the left action of on defined by and . Recall that denotes the differential operator appearing in Equation (1.1).
There exist difference operators with the following properties.
The differential equation is satisfied if and only if
The left-hand side operator is of the form where and for all .
we have (this expression is to be interpreted as a polynomial identity in ). In particular, depends only on and satisfies the same symmetry property as .
Note that , as defined in Eq. (2.10), may be interpreted as an operator from the symmetric sequences to the sequences defined only for nonzero . A sloppy but perhaps more intuitive statement of the main point of Theorem 2.1 would be: “ if and only if , up to some integration constants”.
Assume . Benoit and Salvy [4, Theorem 1] give a simple proof that (2.9) holds for some . The fact that and can actually be taken to have polynomial coefficients and satisfy the properties listed in the last two items then follows from the explicit construction discussed in Section 4.1 of their article, based on Paszkowski’s algorithm [49, 35]. More precisely, multiplying both members of [4, Eq. (17)] by yields a recurrence of the prescribed form. The recurrence has polynomial coefficients since . Rebillard’s thesis [52, Section 4.1] contains detailed proofs of this last observation and of all assertions of Item ii. Note that, although Rebillard’s and Benoit and Salvy’s works are closest to the formalism we use, several of these results actually go back to [49, 35, 36].
There remains to prove the “if” direction. Consider sequences such that , and let be the Chebyshev coefficient sequence of the (analytic) function . We then have by the previous argument. This implies , whence finally by Lemma 2.2 below. ∎
The restriction to of the operator from Theorem 2.1 is injective.
With the notation of Theorem 2.1, we show by induction on that
First, we have since any sequence belonging to converges to zero as . Now assume that (2.11) holds, and let be such that for . Let . Observe that is stable under the action of , so . Since , we have
Hence, for , it holds that
Unless is ultimately zero, this implies that as , which contradicts the fact that . It follows that for large enough, and, using (2.12) again, that as soon as . Applying the hypothesis (2.11) concludes the induction. ∎
An easy-to-explain way of computing a recurrence of the form (2.9) is as follows. We first perform the change of variable in the differential equation (1.1). Then, we compute a recurrence on the Laurent coefficients of by the classical (Frobenius) method.
The function satisfies the homogeneous equation . The substitutions
yield (after clearing common factors and denominators)
We then set and extract the coefficient of (which amounts to replacing by and by ) to get the recurrence
Benoit and Salvy  give a unified presentation of several alternative algorithms, including Paszkowski’s, by interpreting them as various ways to perform the substitution , in a suitable non-commutative division algebra. In our setting where the operator is nonsingular over , they prove that all these algorithms compute the same operator .
As applied in Example 2.3, the method based on setting in the differential equation does not always yield the same operator as Paszkowski’s algorithm. It can be modified to do so as follows: instead of clearing the denominator of the differential equation in given by the rational substitution, move this denominator to the right-hand side, translate both members into recurrences, and then remove a possible common left divisor of the resulting operators .
that is, . In particular, if is a solution of a homogeneous Chebyshev recurrence, then so is , and is a symmetric solution. Not all solutions are symmetric. For instance, the differential equation corresponds to the recurrence which allows for .
2.3. Solutions of the Chebyshev Recurrence
Several difficulties arise when trying to use the Chebyshev recurrence to compute the Chebyshev coefficients.
A first issue is related to initial conditions. Here it may be worth contrasting the situation with the more familiar case of the solution of differential equations in power series. Unlike the first few Taylor coefficients of , the Chebyshev coefficients that could serve as initial conditions for the recurrence are not related in any direct way to initial or boundary conditions of the differential equation. In particular, as can be seen from Theorem 2.1 above, the order of the recurrence is larger than that of the differential equation except for degenerate cases. Hence we need to somehow “obtain more initial values for the recurrence than we naturally have at hand” 111Nevertheless, the recurrence (2.9) shows that the Chebyshev coefficients of a D-finite function are rational linear combinations of a finite number of integrals of the form (2.3). Computing these coefficients efficiently with high accuracy is an interesting problem to which we hope to come back in future work. See Benoit  for some results..
Next, also in contrast to the case of power series, the leading and trailing coefficients of the recurrence (2.9) may vanish for arbitrarily large values of even though the differential equation (1.1) is nonsingular. The zeroes of are called the leading singularities of (2.9), those of , its trailing singularities. In the case of Chebyshev recurrences, leading and trailing singularity sets are opposite of each other.
One reason for the presence of (trailing) singularities is clear: if a polynomial of degree is a solution of , then necessarily . However, even differential equations without polynomial solutions can have arbitrarily large leading and trailing singularities, as shown by the following example.
For all , the Chebyshev recurrence relation associated to the differential equation , namely
admits the leading singularity . For , the differential equation has no polynomial solution.
We do however have some control over the singularities.
With the notations of Theorem 2.1, the coefficients of the Chebyshev recurrence satisfy the relations
with for . In particular, is zero for all .
We proceed by induction on . When , assertion (2.13) reduces to , which follows from the second item of Theorem 2.1. In particular, this proves the result for . Now let and assume that the proposition holds when has order . Write where and is a differential operator of order at most . Letting be the Chebyshev recurrence operator associated to , we then have 
where the last term denotes the evaluation of at . Since
The case having already been dealt with, assume . Since and is a polynomial, it follows by extracting the coefficient of in the last equality and evaluating at that
Now for by the induction hypothesis, and the term involving vanishes for and . In each case, we obtain . ∎
Let be the Chebyshev recurrence operator associated to . The image by of a symmetric sequence satisfies for .
it follows from Proposition 2.8 with and that
that is, . ∎
Last but not least, Chebyshev recurrences always admit divergent solution sequences. Divergent solutions do not correspond to the expansions of solutions of the differential equation the recurrence comes from.
The Chebyshev recurrence associated to the equation is
In terms of the modified Bessel functions and , a basis of solutions of the recurrence is given by the sequences and . The former is the coefficient sequence of the Chebyshev expansion of the exponential function and decreases as . The later satisfies .
3. Convergent and Divergent Solutions
3.1. Elements of Birkhoff-Trjitzinsky Theory
Before studying in more detail the convergent and divergent solutions of the Chebyshev recurrence relation, we recall some elements of the asymptotic theory of linear difference equations. Much of the presentation is based on Wimp’s book [65, Appendix B], to which we refer the reader for more information.
For all , , , we call the formal expansion
a formal asymptotic series (FAS). The set of all FAS is denoted by .
Formal asymptotic series are to be interpreted as asymptotic expansions of sequences as . The product of two FAS is defined in the obvious way and is again an FAS. The same goes for the substitution for fixed , using identities such as . The sum of two FAS is not always an FAS, but that of two FAS sharing the same parameters is. Thus, it makes sense to say that an FAS satisfies a recurrence
with formal series coefficients of the form
Also, given FAS , the Casoratian
belongs to as well.
Following Wimp, we say that are formally linearly independent when their Casoratian is nonzero. Note that the elements of any subset of are then formally linearly independent as well. Indeed, it can be checked by induction on that FAS are formally linearly dependent if and only if there exists a relation of the form where the are FAS such that222Like Wimp, but unlike most authors, we consider recurrences rather than difference equations. Accordingly, we forbid factors of the form with in (3.1), so that the are actually constants in our setting. .
The FAS (3.1) is said to be an asymptotic expansion of a sequence , and we write , when for any truncation order , the relation
holds as .
The following fundamental result is known as the Birkhoff-Trjitzinsky theorem, or “main asymptotic existence theorem” for linear recurrences. It will be the starting point of our analysis of the computation of “convergent” solutions of the Chebyshev recurrence by backward recurrence.
the (formal) recurrence (3.2) possesses a system of formally linearly independent FAS solutions;
We note that many expositions of the Birkhoff-Trjitzinsky theorem warn about possible major gaps in its original proof. However, the consensus among specialists now appears to be that these issues have been resolved in modern proofs [27, 62]. Besides, under mild additional assumptions on the Chebyshev recurrence, all the information needed in our analysis is already provided by the more elementary Perron-Kreuser theorem333The Perron-Kreuser theorem yields the existence of a basis of solutions such that , under the assumption that . It does not require that the coefficients of (3.4) admit full asymptotic expansions, which makes it stronger than Theorem 3.3 in some respects. (cf. [26, 41, 43]) or its extensions by Schäfke . See also Immink  and the references therein for an alternative approach in the case of recurrences with polynomial coefficients, originating in unpublished work by Ramis.
Also observe that for any subfamily of the sequences from Theorem 3.3, the matrix is nonsingular for large . In particular, the can vanish only for finitely many . The more precise statement below will be useful in the sequel.
Assume that the sequences admit formally linearly independent asymptotic expansions of the form (3.1), with , . Then the Casorati determinant
for some , , and .
Write . The formal linear independence hypothesis means that , and hence , admit nonzero FAS as asymptotic expansions. Additionally,
has at most polynomial growth, so that the leading term of its asymptotic expansion must be of the form . ∎
3.2. Newton Polygon of a Chebyshev Recurrence
The formal solutions described in Theorem 3.3 may be constructed algorithmically using methods going back to Poincaré  and developed by many authors. See in particular Adams  and Birkhoff  for early history, Tournier  for a comparison of several methods from a Computer Algebra perspective, and Balser and Bothner  for a modern algorithm as well as more references.
Here, we are mostly interested in the parameters and that control the “exponential” growth rate of the solutions. We briefly recall how the possible values of these parameters are read off the recurrence using the method of Newton polygons. Consider again the Chebyshev recurrence operator
where denotes the leading coefficient of . Note that the degrees of the sum to . Let
be the sequence of all roots of the polynomials , with multiplicities, the roots corresponding to distinct edges being written in the order of increasing and the roots of each in that of increasing modulus. For all , let be the slope of the edge associated to . (Thus, each is repeated a number of times equal to the horizontal length of the corresponding edge, and we have .)
How does this relate to the asymptotics of Chebyshev series? Assume that is the leading factor of some FAS solution of . It is not too hard to see that, in order for asymptotically dominant terms of to cancel out, must be among the slopes of the Newton polygon of , and must be a root of the characteristic equation of the corresponding edge. This gives all possible values of and . Conversely, the construction behind Theorem 3.3 (i) yields a number of linearly independent FAS with given and equal to the multiplicity of as a root of the characteristic equation of the edge of slope . In the case of Chebyshev recurrences, the Newton polygon has the following symmetry property.
The slopes of the Newton polygon of and the roots of its characteristic equations satisfy and for all . In addition, none of the roots associated to the horizontal edge (if there is one) has modulus .
By Theorem 2.1, the coefficients of are related by . Hence, the Newton polygon is symmetric with respect to the vertical axis, and for all . Now fix , and let be the edge of slope . The characteristic equation of reads
where denotes the leading coefficient of a polynomial . Using the relation for lying on , we get
and hence .
There remains to prove that implies . Under the change of variable , the leading term with respect to of is . (The leading term is well-defined because the commutation relation between and preserves degrees.) Therefore, the characteristic equation associated to the slope (when there is one) of the recurrence operator obtained by changing into and into is
where is the leading coefficient of (1.1). Since is a right factor of , the characteristic polynomial associated to in the Newton polygon of <