Rigorous Uniform Approximation of D-finite Functions Using Chebyshev Expansions

# Rigorous Uniform Approximation of D-finite Functions Using Chebyshev Expansions

Alexandre Benoit Alexandre Benoit; Éducation nationale; France http://alexandre.benoit.83.free.fr/ Mioara Joldes Mioara Joldes; CNRS, LAAS, 7 Avenue du Colonel Roche, 31077 Toulouse, Cedex 4, France http://homepages.laas.fr/mmjoldes/  and  Marc Mezzarobba Marc Mezzarobba; CNRS, UMR 7606, LIP6, F-75005, Paris, France; Sorbonne Universités, UPMC Univ Paris 06, UMR 7606, LIP6, F-75005, Paris, France http://marc.mezzarobba.net/
###### Abstract.

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,
This research was partly supported by the Austria Science Fund (FWF) grants P22748-N18 and Y464-N18, and before that by the MSR-Inria joint research center.
This work is in the public domain. As such, it is not subject to copyright. Where it is not legally possible to consider this work as released into the public domain, any entity is granted the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

## 1. Introduction

### 1.1. Background

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 [56]. 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.

###### Problem 1.1.

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ș [10]. 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.

### 1.2. Setting

To specify the D-finite function , we fix a linear homogeneous differential equation of order  with polynomial coefficients

 (1.1)

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

 (1.2) y(i)(0)=ℓi,0≤i≤r−1.

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.

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 [65]. 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 [9] 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.

Note that earlier versions of the present work appeared as part of the authors’ PhD theses [3, 29, 42].

### 1.4. Outline

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 [16]. 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 [53] or Mason and Handscomb’s [39] 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

 (2.1) ∑n′unTn=u02T0+u1T1+u2T2+⋯,

we write Chebyshev series as

 (2.2) ∞∑n=−∞cnTn(x),c−n=cn.

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

 (2.3) cn=1π∫1−1f(x)Tn(x)√1−x2dx.

for all . The series (2.2) converges to  in the sense for all (see for example [39, Chap. 5.3.1]). We denote by the associated orthogonal projection on the subspace of polynomials of degree at most .

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

 (2.4) Er={x∈C:|x+√x2−1|

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

 Ar={z∈C:r−1<|z|

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

 (∀n∈N)(cn=c−n)and(∃α<1)(cn=On→∞(αn)).

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]

 (2.5) ∥f−πd(f)∥∞≤(4π2ln(d+1)+4)∥f−p∗d∥∞

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 [13]. 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

 (2.6) 2xTn(x)=Tn−1(x)+Tn+1(x),

as well as the mixed differential-difference relation

 (2.7) 2(1−x2)T′n(x)=n(Tn−1(x)−Tn+1(x))

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 [49]. 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 [24].

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

 (2.8) Sλ=λS(λ∈Q),Sn=(n+1)S.

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).

###### Theorem 2.1.

[49, 35, 36, 52, 4] Let be analytic functions on some complex neighborhood of the segment , with Chebyshev expansions

 u(x)=∞∑n=−∞unTn(x),v(x)=∞∑n=−∞vnTn(x).

There exist difference operators  with the following properties.

1. The differential equation is satisfied if and only if

 (2.9) P⋅(un)=Q⋅(vn).
2. The left-hand side operator  is of the form where and for all .

3. Letting

 (2.10) δr(n)=2rr−1∏i=−r+1(n−i),I=12n(S−1−S),

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”.

###### Proof.

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. ∎

###### Lemma 2.2.

The restriction to  of the operator  from Theorem 2.1 is injective.

###### Proof.

With the notation of Theorem 2.1, we show by induction on  that

 (2.11) (v∈C)∧(|n|≥r⟹(Qr⋅v)n=0)⟹v=0.

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

 2nQr+1 =δr+1(n)(S−1−S)Ir =2((n+r)(n+r−1)S−1δr(n)−(n−r)(n−r+1)Sδr(n))Ir =2((n+r)(n+r−1)S−1−(n−r)(n−r+1)S)Qr.

Hence, for , it holds that

 (2.12) (n+r)(n+r−1)wn−1=(n−r)(n−r+1)wn+1.

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.

###### Example 2.3.

The function satisfies the homogeneous equation . The substitutions

 x=z+z−12ddx=2zz−z−1ddz

yield (after clearing common factors and denominators)

 (z+1)(z−1)(z4+18z2+1)^y′′(z)+2(z4−2z2−19)z^y′(z)=0.

We then set and extract the coefficient of  (which amounts to replacing  by  and  by ) to get the recurrence

 (n−2)(n−3)cn−3+(n−1)(17n−38)cn−1−(n+1)(17n+38)cn+1−(n+2)(n+3)cn+3=0.

Benoit and Salvy [4] 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 .

###### Remark 2.4.

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 .

###### Definition 2.5.

Following Rebillard, we call the recurrence relation (2.9) computed by Paszkowski’s algorithm (or any equivalent method) the Chebyshev recurrence associated to the differential equation (1.1).

###### Remark 2.6.

By Theorem 2.1(ii) and with its notation, for any sequence , we have the equalities

 ∀n,∑kbk(n)un+k=−∑kb−k(−n)un+k=−∑kbk(−n)u−n−k,

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 [3] 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.

###### Example 2.7.

For all , the Chebyshev recurrence relation associated to the differential equation , namely

 (n+1)(n−k−3)cn−3+(n−1)(5n+k+7)cn−1+8n(n+1)(n−1)cn−(n+1)(5n−k−7)cn+1−(n−1)(n+k+3)cn+3=0,

admits the leading singularity . For , the differential equation has no polynomial solution.

We do however have some control over the singularities.

###### Proposition 2.8.

With the notations of Theorem 2.1, the coefficients of the Chebyshev recurrence satisfy the relations

 (2.13) bj−i(−j)=−bj+i(−j),|j|≤r−1,i∈N,

with  for . In particular, is zero for all .

###### Proof.

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 [4]

 (2.14) δr(n)−1P=Iδr−1(n)−1P♭+pr(12(S+S−1))

where the last term denotes the evaluation of  at . Since

 Iδr−1(n)−1=(nδr(n))−1((n−r+2)(n−r+1)S−1−(n+r−2)(n+r−1)S)

by the commutation rule (2.8), relation (2.14) rewrites as

 P =1n∑k((n−r+2)(n−r+1)b♭k+1(n−1) −(n+r−2)(n+r−1)b♭k−1(n+1))Sk +δr(n)pr(12(S+S−1)).

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

 (2.15) −jbk(−j)=(j+r−2)(j+r−1)b♭k+1(−j−1)−(j−r+2)(j−r+1)b♭k−1(−j+1).

Now  for by the induction hypothesis, and the term involving  vanishes for and . In each case, we obtain . ∎

###### Corollary 2.9.

Let  be the Chebyshev recurrence operator associated to . The image by  of a symmetric sequence  satisfies for .

###### Proof.

Since

 (P⋅u)n=∑k∈Zbk(n)un+k=∑i∈Zbi−n(n)ui,

it follows from Proposition 2.8 with  and that

 ∑i∈Zbi−n(n)ui=−∑i∈Zb−i−n(n)ui=−∑i∈Zbi−n(n)ui,

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.

###### Example 2.10.

The Chebyshev recurrence associated to the equation  is

 (P⋅u)n=u(n+1)+2nu(n)−u(n−1)=0.

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.

###### Definition 3.1.

For all , , , we call the formal expansion

 (3.1) ¯u(n)=n!καneπ(n)J∑j=0(lnn)j∞∑i=0βj,inθ−i/ρ

where

 π(n)=π1n1/ρ+⋯+πρ−1n(ρ−1)/ρ

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

 (3.2) ¯bs(n)¯u(n+s)+⋯+¯b0(n)¯u(n)=0

with formal series coefficients of the form

 (3.3) ¯bk(n)=nτk/ω(βk,0+βk,1n−1/ω+βk,2n−2/ω+⋯)∈C((n−1/ω)).

Also, given  FAS , the Casoratian

 C(n)=det(¯uj(n+i))0≤i,j

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. .

###### Definition 3.2.

The FAS (3.1) is said to be an asymptotic expansion of a sequence , and we write , when for any truncation order , the relation

 un=n!καneπ(n)J∑j=0(lnn)j(I−1∑i=0βj,inθ−i/ρ+O(nθ−I/ρ))

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.

###### Theorem 3.3.

[6, 7, 61, 27] Consider a linear recurrence

 (3.4) bs(n)un+s+⋯+b0(n)un=0

whose coefficients admit asymptotic expansions (in the sense of Definition 3.2) of the form (3.3) for some integer . Then,

1. the (formal) recurrence (3.2) possesses a system of  formally linearly independent FAS solutions;

2. for any formally linearly independent solutions of (3.2), there exists complex sequences defined in some neighborhood of infinity, with the property that for all , and such that is a basis of the solution space of (3.4) for .

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 [55]. See also Immink [28] 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.

###### Lemma 3.4.

Assume that the sequences admit formally linearly independent asymptotic expansions of the form (3.1), with , . Then the Casorati determinant

 C(n)=∣∣ ∣ ∣ ∣ ∣∣e0,ne1,n⋯es−1,ne0,n+1es−1,n+1⋮⋮e0,n+s−1e1,n+s−1⋯es−1,n+s−1∣∣ ∣ ∣ ∣ ∣∣

satisfies

 C(n)=βe0,ne1,n⋯es−1,nnθ((lnn)λ+O((lnn)λ−1)),n→∞,

for some , , and .

###### Proof.

Write . The formal linear independence hypothesis means that , and hence , admit nonzero FAS as asymptotic expansions. Additionally,

 C′(n)=det(ej,n+iej,n)0≤i,j

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é [50] and developed by many authors. See in particular Adams [1] and Birkhoff [6] for early history, Tournier [57] for a comparison of several methods from a Computer Algebra perspective, and Balser and Bothner [2] 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

 P=b−s(n)S−s+⋯+b0(n)+⋯+bs(n)Ss

from Section 2.2. The Newton polygon of  is defined as the lower convex hull of the points (see Figure 1). To each edge of the polygon is attached a characteristic equation

 χi(α)=∑k:pk∈[pi,pj]lc(bk)αk−i,

where denotes the leading coefficient of . Note that the degrees of the sum to . Let

 αs,αs−1,…,α1,α−1,…,α−s+1,α−s

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.

###### Proposition 3.5.

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 .

###### Proof.

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

 χi(α)=∑k:pk∈ϵilc(bk)αk−ℓ(i)=∑k:pk∈ϵi(−1)1+degbklc(b−k)αk−ℓ(i),

where denotes the leading coefficient of a polynomial . Using the relation for lying on , we get

 χi(α) =±∑k:pk∈ϵ−i(−1)κi(−k+ℓ(−i))lc(bk)α−k+ℓ(−i) =±αℓ(i)−ℓ(−i)χ−i((−1)−κiα−1),

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

 χhoriz(α):=∑k:degpk=maxidegpilc(bk)αk−i=ar(α+α−12),

where  is the leading coefficient of (1.1). Since  is a right factor of , the characteristic polynomial associated to  in the Newton polygon of <