Alexandre Benoit

Mioara Joldes
email: joldes@laas.fr

Marc Mezzarobba

RIGOROUS UNIFORM APPROXIMATION OF D-FINITE FUNCTIONS USING CHEBYSHEV EXPANSIONS

Keywords: 2000 Mathematics Subject Classification. 68W30, 33C45, 65L05, 65G20, Rigorous computing, computer algebra, complexity, D-finite functions, recurrence relation, Chebyshev series, Clenshaw method, Miller algorithm, asymptotics, functional enclosure

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

as mathematical physics are so-called D-finite functions, that is, solutions of linear ordinary differential equations (LODE) with polynomial coefficients [START_REF] Stanley | Differentiably finite power series[END_REF]. 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 [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF][START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF][START_REF] Kauers | The holonomic toolkit[END_REF]. The present article is devoted to the following problem. Approximations over other real or complex segments (written on the Chebyshev basis adapted to the segment) are reduced to approximations on [-1, 1] 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 y 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 y. 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 [START_REF] Trefethen | Computing numerically with functions instead of numbers[END_REF][START_REF] Driscoll | The chebop system for automatic solution of differential equations[END_REF] is a popular numerical computation system based on this idea.

In a yet more general setting, Epstein, Miranker and Rivlin developed a socalled "ultra-arithmetic" for functions that parallels floating-point arithmetic for real numbers [START_REF] Epstein | Ultra-arithmetic I: function data types[END_REF][START_REF] Epstein | Ultra-arithmetic II: intervals of polynomials[END_REF][START_REF] Kaucher | Self-validating numerics for function space problems[END_REF]. 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 Joldes , [START_REF] Brisebarre | Chebyshev interpolation polynomial-based tools for rigorous computing[END_REF]. Part of the motivation for Problem 1.1 is to allow one to use arbitrary Dfinite 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 [START_REF] Moore | Methods and applications of interval analysis[END_REF][START_REF] Kaucher | Self-validating numerics for function space problems[END_REF][START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF][START_REF] Neumaier | Taylor forms -Use and limits[END_REF][START_REF] Tucker | Validated Numerics: A Short Introduction to Rigorous Computations[END_REF]. 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 y, we fix a linear homogeneous differential equation of order r with polynomial coefficients (1.1) L • y = a r y (r) + a r-1 y (r-1) + • • • + a 0 y = 0, a i ∈ É[x].

We also write L = a r ∂ r + • • • + a 1 ∂ + a 0 . Up to a change of variable, we assume that we are seeking a polynomial approximation of a solution y of (1.1) over the interval [-1, 1]. The uniform norm on this interval is denoted by • ∞ . We also assume that a r (x) = 0 for x ∈ [-1, 1], so that (by Cauchy's existence theorem for complex LODE) all solutions of (1.1) are analytic on [-1, 1]. Besides the operator L, we are given r 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 i, j = {i, i + 1, . . . , j}.

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 a i coefficients of the operator L; a r (0) = 0 p. 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 p and error bound B 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 [START_REF] Wimp | Computation with Recurrence Relations[END_REF]. 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 n may be done in time O(n(ln n) O( 1) ) (see for instance Brent and Zimmermann's book [START_REF] Brent | Modern Computer Arithmetic[END_REF] for detail). In general, the maximum bit size of the numbers we manipulate is roughly the same as that of the coefficients of p, which may be checked to be O(d ln d) 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 y. 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. 18) and 6.6 (p. 33), 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 [START_REF] Benoit | Algorithmique semi-numérique rapide des séries de Tchebychev[END_REF][START_REF] Joldes | Rigorous polynomial approximations and applications[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF]. 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 [START_REF] Fox | Chebyshev polynomials in numerical analysis[END_REF]Chap. 5] of Clenshaw's algorithm [START_REF] Clenshaw | The numerical solution of linear differential equations in Chebyshev series[END_REF]. 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 B. We conclude with examples of error bounds obtained using our implementation of the validation algorithm, and some open questions.

Chebyshev Expansions of D-finite Functions

(x) = 1/ √ 1 -x 2
, and hence a Hilbert basis of the space L 2 (w). (We refer the reader to books such as Rivlin's [START_REF] Rivlin | The Chebyshev polynomials[END_REF] or Mason and Handscomb's [START_REF] Mason | Chebyshev polynomials[END_REF] for proofs of the results collected in this section.)

Expansions of functions f ∈ L 2 (w) on this basis are known as Chebyshev series. Instead of the more common

(2.1) n ′ u n T n = u 0 2 T 0 + u 1 T 1 + u 2 T 2 + • • • ,
we write Chebyshev series as

(2.2) ∞ n=-∞ c n T n (x), c -n = c n .
This choice makes the link between Chebyshev and Laurent expansions as well as the action of recurrence operators on the c n (both discussed below) more transparent.

The Chebyshev coefficients c n = 1 2 u n of the expansion of a function f are given by

(2.3) c n = 1 π 1 -1 f (x)T n (x) √ 1 -x 2 dx.
for all n ∈ . The series (2.2) converges to f in the L 2 (w) sense for all f ∈ L 2 (w) (see for example [START_REF] Mason | Chebyshev polynomials[END_REF]Chap. 5.3.1]). We denote by π d : f → d n=-d c n T n the associated orthogonal projection on the subspace of polynomials of degree at most d. Now assume that f is a solution of Equation (1.1). As such, it may be analytically continued to any domain U ⊂ that does not contain any singular point of the equation. Let

2

) around the unit circle. The transformation x = z+z -1 2 sending f (x) to f (z) is known as the inverse Joukowski transform. It maps the elliptic disk E r to the annulus

A r = {z ∈ : r -1 < |z| < r}. The formula T n (cos θ) = cos(nθ) translates into T n ( z+z -1 2 ) = z n +z -n

2

. The coefficients c n are also related to those of the Fourier cosine expansion of θ → f (cos θ).

Let C ⊂ be the vector space of doubly infinite sequences (c n ) n∈ such that

(∀n ∈ AE)(c n = c -n ) and (∃α < 1)(c n = O n→∞ (α n )).
The sequence of Chebyshev coefficients of a function f that is analytic on some complex neighborhood of [-1, 1] belongs to C. Conversely, for all c ∈ C, the function series ∞ n=-∞ c n T n (x) converges uniformly on (some neighborhood of) [-1, 1] to an analytic function f (x).

Truncated Chebyshev series are near-minimax approximations: indeed, they satisfy [START_REF] Trefethen | Approximation theory and approximation practice[END_REF]Theorem 16.1]

(2.5) f -π d (f ) ∞ ≤ 4 π 2 ln(d + 1) + 4 f -p * d ∞
where p * d is the polynomial of degree at most d that minimizes f -p ∞ . Even though p * d itself can be computed to arbitrary precision using the Remez algorithm [START_REF] Cheney | Introduction to approximation theory[END_REF]Chap. 3], Equation (2.5) shows that we do not lose much by replacing it by π d (f ). Moreover, tighter approximations are typically hard to validate without resorting to intermediate approximations of higher degree [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. 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 (c n ) when f is a D-finite function.

2.2. The Chebyshev Recurrence Relation. The polynomials T n satisfy the three-term recurrence (2.6) 2xT

n (x) = T n-1 (x) + T n+1 (x),
as well as the mixed differential-difference relation

(2.7) 2(1 -x 2 )T ′ n (x) = n(T n-1 (x) -T n+1 (x)) which translates into the integration formula 2nc n = c ′ n-1 -c ′ n+1 where c ′ n T n = ( c n T n ) ′ .
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 [START_REF] Fox | Chebyshev methods for ordinary differential equations[END_REF][START_REF] Fox | Chebyshev polynomials in numerical analysis[END_REF] in special cases and later proved in general by Paszkowski [START_REF] Paszkowski | Zastosowania numeryczne wielomianow i szeregow Czebyszewa[END_REF]. 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 [START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF][START_REF] Lewanowicz | A new approach to the problem of constructing recurrence relations for the Jacobi coefficients[END_REF]). The automatic determination of this recurrence in a symbolic computation system was first studied by Geddes [START_REF] Geddes | Symbolic computation of recurrence equations for the Chebyshev series solution of linear ODE's[END_REF].

The following theorem summarizes results regarding this recurrence, extracted from existing work [START_REF] Paszkowski | Zastosowania numeryczne wielomianow i szeregow Czebyszewa[END_REF][START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF][START_REF] Lewanowicz | A new approach to the problem of constructing recurrence relations for the Jacobi coefficients[END_REF][START_REF] Rebillard | Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF][START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF] and extended to fit our purposes. Here and in the sequel, we denote by É(n) S, S -1 the skew Laurent polynomial ring over É(n) in the indeterminate S, subject to the commutation rules Theorem 2.1. [START_REF] Paszkowski | Zastosowania numeryczne wielomianow i szeregow Czebyszewa[END_REF][START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF][START_REF] Lewanowicz | A new approach to the problem of constructing recurrence relations for the Jacobi coefficients[END_REF][START_REF] Rebillard | Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF][START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF] Let u, v be analytic functions on some complex neighborhood of the segment [-1, 1], with Chebyshev expansions

u(x) = ∞ n=-∞ u n T n (x), v(x) = ∞ n=-∞ v n T n (x).
There exist difference operators P, Q ∈ É[n] S, S -1 with the following properties.

(i) The differential equation L • u(x) = v(x) is satisfied if and only if (2.9) P • (u n ) = Q • (v n ).
(ii) The left-hand side operator P is of the form

P = s k=-s b k (n)S k where s = r + max i (deg a i ) and b -k (-n) = -b k (n) for all k. (iii) Letting (2.10) δ r (n) = 2 r r-1 i=-r+1 (n -i), I = 1 2n (S -1 -S),
we have Q = Q r = δ r (n)I r (this expression is to be interpreted as a polynomial identity in É(n) S, S -1 ). In particular, Q depends only on r and satisfies the same symmetry property as P .

Note that I, as defined in Eq. (2.10), may be interpreted as an operator from the symmetric sequences (u |n| ) n∈ to the sequences (u n ) n∈ \{0} defined only for nonzero n. A sloppy but perhaps more intuitive statement of the main point of Theorem 2.1 would be: "( ) r L • u = w if and only if δ r (n)P • u = w, up to some integration constants".

Proof. Assume L • u = v. Benoit and Salvy [4, Theorem 1] give a simple proof that (2.9) holds for some P, Q ∈ É(n) S, S -1 . The fact that P and Q 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 [START_REF] Paszkowski | Zastosowania numeryczne wielomianow i szeregow Czebyszewa[END_REF][START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF]. More precisely, multiplying both members of [START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF]Eq. (17)] by δ r (n) yields a recurrence of the prescribed form. The recurrence has polynomial coefficients since δ r (n)I r ∈ Q S, S -1 . Rebillard's thesis [START_REF] Rebillard | Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF]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 [START_REF] Paszkowski | Zastosowania numeryczne wielomianow i szeregow Czebyszewa[END_REF][START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF][START_REF] Lewanowicz | A new approach to the problem of constructing recurrence relations for the Jacobi coefficients[END_REF].

There remains to prove the "if" direction. Consider sequences u, v ∈ C such that Proof. With the notation of Theorem 2.1, we show by induction on r ≥ 1 that

P • u = Q • v,
(2.11) (v ∈ C) ∧ |n| ≥ r =⇒ (Q r • v) n = 0 =⇒ v = 0.
First, we have (ker Q 1 ) ∩ C = {0} since any sequence belonging to C converges to zero as n → ±∞. Now assume that (2.11) holds, and let v ∈ C be such that

(Q r+1 • v) n = 0 for |n| ≥ r + 1. Let w = Q r • v. Observe that C is stable under the action of Q(n) S, S -1 , so w ∈ C. Since r ≥ 1, we have 2n Q r+1 = δ r+1 (n)(S -1 -S)I r = 2 ((n + r)(n + r -1)S -1 δ r (n) -(n -r)(n -r + 1)Sδ r (n))I r = 2 ((n + r)(n + r -1)S -1 -(n -r)(n -r + 1)S)Q r .
Hence, for |n| ≥ r + 1, it holds that (2.12)

(n + r)(n + r -1)w n-1 = (n -r)(n -r + 1)w n+1 .
Unless w n is ultimately zero, this implies that w n+1 /w n-1 → 1 as n → ∞, which contradicts the fact that w ∈ C. It follows that w n = 0 for |n| large enough, and, using (2.12) again, that w n = 0 as soon as |n| ≥ r. 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 x = 1 2 (z + z -1 ) in the differential equation (1.1). Then, we compute a recurrence on the Laurent coefficients of û(z) = u(x) by the classical (Frobenius) method.

Example 2.3. The function y(x) = arctan(x/2) satisfies the homogeneous equation (x 2 + 4) y ′′ (x) + 2x y ′ (x) = 0. The substitutions

x = z + z -1 2 d dx = 2z z -z -1
d dz yield (after clearing common factors and denominators)

(z + 1)(z -1)(z 4 + 18z 2 + 1) ŷ′′ (z) + 2(z 4 -2z 2 -19)z ŷ′ (z) = 0.
We then set ŷ(z) = ∞ n=-∞ c n z n and extract the coefficient of z n (which amounts to replacing z by S -1 and z d dz by n) to get the recurrence

(n -2)(n -3)c n-3 + (n -1)(17n -38)c n-1 -(n + 1)(17n + 38)c n+1 -(n + 2)(n + 3)c n+3 = 0.
Benoit and Salvy [START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF] give a unified presentation of several alternative algorithms, including Paszkowski's, by interpreting them as various ways to perform the substitution x → 1 2 (S + S -1 ), d dx → (S -S -1 ) -1 (2n) in a suitable non-commutative division algebra. In our setting where the operator L is nonsingular over [-1, 1], they prove that all these algorithms compute the same operator P . Remark 2.4. As applied in Example 2.3, the method based on setting x = 1 2 (z +z -1 ) 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 z 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 P, Q ∈ É(n) S, S -1 . 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 (u n ) n∈ , we have the equalities

∀n, k b k (n)u n+k = - k b -k (-n)u n+k = - k b k (-n)u -n-k , that is, P • (u n ) n∈ = -P • (u -n ) n∈ . In particular, if (u n ) n∈ is a solution of a homogeneous Chebyshev recurrence, then so is (u -n ) n∈ , and (u n + u -n ) is a symmetric solution. Not all solutions are symmetric. For instance, the differential equation y ′ (x) = xy(x) corresponds to the recurrence -c n-2 + 4n c n + c n+2 = 0 which allows for u -2 = 3, u -1 = 12, u 0 = 1, u 1 = 2, u 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 y, the Chebyshev coefficients c 0 , c 1 , . . . 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 2s 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" 1 .

Next, also in contrast to the case of power series, the leading and trailing coefficients b ±s of the recurrence (2.9) may vanish for arbitrarily large values of n even though the differential equation (1.1) is nonsingular. The zeroes of b s (n -s) are called the leading singularities of (2.9), those of b -s (n + s), 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 y = y |n| T n of degree d is a solution of L • y = 0, then necessarily b -s (d + s) = 0. 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 k ∈ , the Chebyshev recurrence relation associated to the differential equation

y ′′ (x) + (x 2 + 1) y ′ (x) -k x y(x) = 0, namely (n + 1)(n -k -3) c n-3 + (n -1)(5n + k + 7) c n-1 + 8n(n + 1)(n -1) c n -(n + 1)(5n -k -7) c n+1 -(n -1)(n + k + 3) c n+3 = 0,
admits the leading singularity n = k. For k = 1, 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) b j-i (-j) = -b j+i (-j), |j| ≤ r -1, i ∈ AE, with b k = 0 for |k| > s. In particular, b s (n) is zero for all n ∈ 1, r -1 .
Proof. We proceed by induction on r. When j = 0, assertion (2.13) reduces to b -i (0) = -b i (0), which follows from the second item of Theorem 2.1. In particular, this proves the result for r = 1. Now let r ≥ 2 and assume that the proposition holds when L has order r -

1. Write L = L ♭ + ∂ r p r (x) where p r ∈ É[x] and L ♭ is a differential operator of order at most r -1. Letting P ♭ = k∈ b ♭ k (n)S k be the Chebyshev recurrence operator associated to L ♭ , we then have [4] (2.14) δ r (n) -1 P = Iδ r-1 (n) -1 P ♭ + p r ( 1 2 (S + S -1 ))
where the last term denotes the evaluation of p r at x = 1 2 (S + S -1 ). 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)
1 Nevertheless, 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 [START_REF] Benoit | Algorithmique semi-numérique rapide des séries de Tchebychev[END_REF] for some results.

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

P = 1 n k (n -r + 2)(n -r + 1)b ♭ k+1 (n -1) -(n + r -2)(n + r -1)b ♭ k-1 (n + 1) S k + δ r (n)p r ( 1 2 (S + S -1 )).
The case j = 0 having already been dealt with, assume 0 < |j| < r. Since δ r (-j) = 0 and p r is a polynomial, it follows by extracting the coefficient of S k in the last equality and evaluating at n = -j that j) for |j| < r -1 by the induction hypothesis, and the term involving b ♭ k±1 vanishes for j = ∓(r -1) and j = ∓(r -2). In each case, we obtain b j-i (-j) = -b j+i (-j).

(2.15) -jb k (-j) = (j + r -2)(j + r -1)b ♭ k+1 (-j -1) -(j -r + 2)(j -r + 1)b ♭ k-1 (-j + 1). Now b ♭ j-i (-j) = -b ♭ j+i (-
Corollary 2.9. Let P be the Chebyshev recurrence operator associated to L. The image by P of a symmetric sequence

(u |n| ) n∈ satisfies (P • u) n = 0 for |n| < r. Proof. Since (P • u) n = k∈ b k (n)u n+k = i∈ b i-n (n)u i ,
it follows from Proposition 2.8 with j = -n and |n| < r that

i∈ b i-n (n)u i = - i∈ b -i-n (n)u i = - i∈ b i-n (n) u i , that is, (P • u) n = -(P • u) n .
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 y ′ = y is

(P • u) n = u(n + 1) + 2n u(n) -u(n -1) = 0.
In terms of the modified Bessel functions I ν and K ν , a basis of solutions of the recurrence is given by the sequences (I ν (1)) ν∈ and (K ν (1)) ν∈ . The former is the coefficient sequence of the Chebyshev expansion of the exponential function and decreases as Θ(2 -ν ν! -1 ). The later satisfies K ν (1) = Θ(2 ν (ν -1)!).

Convergent and Divergent Solutions

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 [START_REF] Wimp | Computation with Recurrence Relations[END_REF]Appendix B], to which we refer the reader for more information. Definition 3.1. For all ρ ∈ AE\{0}, J ∈ AE, κ ∈ É, α, π j , θ, β j,i ∈ , we call the formal expansion

(3.1) ū(n) = n! κ α n e π(n) J j=0 (ln n) j ∞ i=0 β j,i n θ-i/ρ where π(n) = π 1 n 1/ρ + • • • + π ρ-1 n (ρ-1
)/ρ a formal asymptotic series (FAS). The set of all FAS is denoted by B.

Formal asymptotic series are to be interpreted as asymptotic expansions of sequences as n → ∞. The product of two FAS is defined in the obvious way and is again an FAS. The same goes for the substitution n → n + k for fixed k ∈ , using identities such as

(n + k) θ = n θ (1 + kθn -1 + • • • ).
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 ū ∈ B satisfies a recurrence

(3.2) bs (n)ū(n + s) + • • • + b0 (n)ū(n) = 0
with formal series coefficients of the form

(3.3) bk (n) = n τ k /ω (β k,0 + β k,1 n -1/ω + β k,2 n -2/ω + • • • ) ∈ ((n -1/ω )).
Also, given s FAS ū0 , . . . , ūs-1 ∈ B, the Casoratian

C(n) = det(ū j (n + i)) 0≤i,j<s
belongs to B as well. Following Wimp, we say that ū1 , . . . , ūs ∈ B are formally linearly independent when their Casoratian is nonzero. Note that the elements of any subset of {ū 1 , . . . , ūs } are then formally linearly independent as well. Indeed, it can be checked by induction on s that s FAS ū1 , . . . , ūs are formally linearly dependent if and only if there exists a relation of the form μ1 (n)ū 1 (n) + • • • + μn (s)ū s (n) = 0 where the μk are FAS such that 2 μk (n + 1) = μk (n). Definition 3.2. The FAS (3.1) is said to be an asymptotic expansion of a sequence (u n ) ∈ AE , and we write u n ∼ ū(n), when for any truncation order I, the relation

u n = n! κ α n e π(n) J j=0 (ln n) j I-1 i=0 β j,i n θ-i/ρ + O(n θ-I/ρ ) holds as n → ∞.
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. [START_REF] Birkhoff | Formal theory of irregular linear difference equations[END_REF][START_REF] Birkhoff | Analytic theory of singular difference equations[END_REF][START_REF] Turrittin | The formal theory of systems of irregular homogeneous linear difference and differential equations[END_REF][START_REF] Immink | Reduction to canonical forms and the Stokes phenomenon in theory of linear difference equations[END_REF] Consider a linear recurrence

(3.4) b s (n)u n+s + • • • + b 0 (n)u n = 0
whose coefficients b 0 , . . . , b s admit asymptotic expansions (in the sense of Definition 3.2) b0 , . . . , bs of the form (3.3) for some integer ω ≥ 1. Then, 2 Like Wimp, but unlike most authors, we consider recurrences rather than difference equations.

Accordingly, we forbid factors of the form e πρn with | Im πρ| > π in (3.1), so that the µ k (n) are actually constants in our setting.

(i) the (formal) recurrence (3.2) possesses a system of s formally linearly independent FAS solutions; (ii) for any s formally linearly independent solutions ē1 , . . . , ēs ∈ B of (3.2), there exists complex sequences e 1 = (e 1,n ) n≥N , . . . , e s = (e s,n ) n≥N defined in some neighborhood of infinity, with the property that e k ∼ ēk for all k, and such that (e 1 , . . . , e s ) is a basis of the solution space of (3.4) for n ≥ N .

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 [START_REF] Immink | Reduction to canonical forms and the Stokes phenomenon in theory of linear difference equations[END_REF][START_REF] Van Der Put | Galois theory of difference equations[END_REF]. 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 theorem 3 (cf. [START_REF] Guelfond | Calcul des différences finies[END_REF][START_REF] Meschkowski | Differenzengleichungen[END_REF][START_REF] Milne-Thomson | The calculus of finite differences[END_REF]) or its extensions by Schäfke [START_REF] Schäfke | Lösungstypen von Differenzengleichungen und Summengleichungen in normierten abelschen Gruppen[END_REF]. See also Immink [START_REF] Immink | On the relation between global properties of linear difference and differential equations with polynomial coefficients[END_REF] 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 (f 1 , . . . , f s ′ ) of the sequences e i from Theorem 3.3, the matrix (f j,n+i ) 1≤i,j≤s ′ is nonsingular for large n. In particular, the e i,n can vanish only for finitely many n. The more precise statement below will be useful in the sequel. Lemma 3.4. Assume that the sequences (e 0,n ) n , . . . , (e s-1,n ) n admit formally linearly independent asymptotic expansions of the form (3.1), with α i ∈ \ {0},

κ i ∈ É. Then the Casorati determinant C(n) = e 0,n e 1,n
C(n) = βe 0,n e 1,n • • • e s-1,n n θ ((ln n) λ + O((ln n) λ-1 )), n → ∞,
for some β ∈ \{0}, θ ∈ , and λ ∈ AE.

Proof. Write C(n) = e 0,n e 1,n • • • e s-1,n C ′ (n).
The formal linear independence hypothesis means that C(n), and hence C ′ (n), admit nonzero FAS as asymptotic expansions. Additionally,

C ′ (n) = det e j,n+i e j,n 0≤i,j<s
has at most polynomial growth, so that the leading term of its asymptotic expansion must be of the form n θ (ln n) λ .

Newton Polygon of a Chebyshev

Recurrence. The formal solutions described in Theorem 3.3 may be constructed algorithmically using methods going back to Poincaré [START_REF] Poincaré | Sur les équations linéaires aux différentielles ordinaires et aux différences finies[END_REF] and developed by many authors. See in particular Adams [START_REF] Adams | On the irregular cases of the linear ordinary difference equation[END_REF] and Birkhoff [START_REF] Birkhoff | Formal theory of irregular linear difference equations[END_REF] for early history, Tournier [START_REF] Tournier | Solutions formelles d'équations différentielles[END_REF] for a comparison of several methods 3 The Perron-Kreuser theorem yields the existence of a basis of solutions such that e i,n+1 /e i,n ∼ αn κ i , under the assumption that

κ i = κ j ⇒ |α i | = |α j |.
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.

S n

α 3 κ 3 α 2 α 1 α -1 α -2 κ 2 = • • • = κ -2 α -3 κ -3 Figure 1.
The Newton polygon of a Chebyshev recurrence.

from a Computer Algebra perspective, and Balser and Bothner [START_REF] Balser | Computation of formal solutions of systems of linear difference equations[END_REF] 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 + • • • + b 0 (n) + • • • + b s (n)S s
from Section 2.2. The Newton polygon of P is defined as the lower convex hull of the points 1). To each edge [p i , p j ] (i < j) of the polygon is attached a characteristic equation

p k = (k, -deg b k ) ∈ Ê 2 (see Figure
χ i (α) = k:p k ∈[pi,pj] lc(b k )α k-i ,
where lc(p) denotes the leading coefficient of p. Note that the degrees of the χ i sum to 2s. Let α s , α s-1 , . . . , α 1 , α -1 , . . . , α -s+1 , α -s be the sequence of all roots of the polynomials χ i , with multiplicities, the roots corresponding to distinct edges being written in the order of increasing i and the roots of each χ i in that of increasing modulus. For all k, let κ k be the slope of the edge associated to α k . (Thus, each κ k is repeated a number of times equal to the horizontal length of the corresponding edge, and we have

κ s ≤ κ s-1 ≤ • • • ≤ κ -s .)
How does this relate to the asymptotics of Chebyshev series? Assume that n! κ α n is the leading factor of some FAS solution ū of P • ū = 0. It is not too hard to see that, in order for asymptotically dominant terms of P • ū to cancel out, κ must be among the slopes of the Newton polygon of P , 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 κ i of the Newton polygon of P and the roots α i of its characteristic equations satisfy κ -i = -κ i and |α -i | = |α i | -1 for all i. In addition, none of the roots associated to the horizontal edge (if there is one) has modulus 1.

Proof. By Theorem 2.1, the coefficients b k of P are related by b -k (n) = -b k (-n).
Hence, the Newton polygon is symmetric with respect to the vertical axis, and κ -i = -κ i for all i. Now fix i, and let ǫ i = [p ℓ(i) , p r(i) ] be the edge of slope κ i . The characteristic equation of ǫ i reads

χ i (α) = k:p k ∈ǫi lc(b k )α k-ℓ(i) = k:p k ∈ǫi (-1) 1+deg b k lc(b -k )α k-ℓ(i) ,
where lc(b) denotes the leading coefficient of a polynomial b. Using the relation

deg b k -deg b ℓ(i) = κ i (k -ℓ(i)) for p k lying on ǫ i , we get χ i (α) = ± k:p k ∈ǫ-i (-1) κi(-k+ℓ(-i)) lc(b k )α -k+ℓ(-i) = ±α ℓ(i)-ℓ(-i) χ -i ((-1) -κi α -1 ), and hence |α -i | = |α i | -1 .
There remains to prove that

κ i = 0 implies |α i | = 1. Under the change of variable x = 1 2 (z + z -1
), the leading term with respect to θ = z d dz of ( d dx ) k is 2 k (z -z -1 ) -k . (The leading term is well-defined because the commutation relation between z and θ preserves degrees.) Therefore, the characteristic equation associated to the slope κ = 0 (when there is one) of the recurrence operator P 1 obtained by changing z into S -1 and θ into n is

χ horiz (α) := k:deg p k =maxi deg pi lc(b k )α k-i = a r α + α -1 2 ,
where a r is the leading coefficient of (1.1). Since P is a right factor of P 1 , the characteristic polynomial associated to κ = 0 in the Newton polygon of P divides χ horiz . But, due to the assumptions stated in Section 1.2, the polynomial a r (x) does not vanish for x ∈ [-1, 1], hence χ horiz (α) = 0 for |α| = 1.

Summing up, the asymptotic structure of the solutions of the Chebyshev recurrence may be described as follows. Similar observations were already made by Rebillard [START_REF] Rebillard | Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF]Chap. 5].

Corollary 3.6. For large enough N , the space of sequences (u n ) n≥N satisfying (P •u) n = 0 ("germs of solution at infinity of the Chebyshev recurrence") has a basis comprising s convergent sequences e 1 , . . . , e s and s divergent sequences e -1 , . . . , e -s , all with formally linearly independent FAS expansions, such that

e i,n = n! κi α n i e o(n) , n → ∞.
In particular, we have

ln |e i,n e -i,n | = o(n) for all i.
Proof. This follows from Theorem 3.3, Proposition 3.5, and the description of a basis of formal solutions at infinity using the Newton polygon.

Computing Approximation Polynomials

4.1. Clenshaw's Algorithm Revisited. At this point, we know that Chebyshev expansions of D-finite functions correspond to the symmetric, convergent solutions of the Chebyshev recurrences introduced in Section 2.2. The question we now face is to compute these solutions efficiently in spite of the various difficulties discussed above. Our algorithm for this task may be viewed as a systematized variant of a method originally due to Clenshaw [START_REF] Clenshaw | The numerical solution of linear differential equations in Chebyshev series[END_REF]. The link between Clenshaw's method 4 and the Chebyshev recurrence was observed long ago by Fox and Parker [START_REF] Fox | Chebyshev methods for ordinary differential equations[END_REF][START_REF] Fox | Chebyshev polynomials in numerical analysis[END_REF] and further discussed by Rebillard [START_REF] Rebillard | Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF]Section 4.1.3]. Based on the properties of the recurrence established in the last two sections, we can turn Clenshaw's method into a true algorithm that applies uniformly to differential equations of arbitrary order and degree. Both Clenshaw's original method and our algorithm are related to Miller's wellknown backward recurrence technique [START_REF] Bickley | Bessel functions. Part II. Functions of positive integer order[END_REF][START_REF] Wimp | Computation with Recurrence Relations[END_REF] to compute minimal ("slowest increasing") solutions of three-term recurrences. Miller's idea is to compute the coefficients u N , u N -1 , . . . , u 0 of a linear recurrence sequence in the backward direction, starting form arbitrary "initial conditions" u N +1 and u N +2 . When N goes to infinity (u N +1 , u N +2 being chosen once and for all), the computed coefficients u 0 , . . . , u N get close to those of a minimal solution with large u 0 , u 1 , in accordance with the intuition that "minimal solutions are the dominant ones when going backwards". This method behaves much better numerically that the standard forward recurrence. But its key feature for our purposes is that it allows one to compute a minimal solution characterized by its minimality plus one normalizing condition instead of two initial values.

Roughly speaking, our method amounts to a "block Miller algorithm" tuned to the special case of Chebyshev recurrences. We use the idea of backward recurrence to approximate the whole subspace of convergent solutions instead of a single minimal one. There remains to take care of the constraints related to the singularities of the recurrence, the symmetry condition and the initial values of the differential equation, all of which is done using linear algebra.

Denote (i) The space S of symmetric sequences (u |n| ) n∈ such that P • u = 0 has dimension s + r. Among the elements of E, these sequences are characterized by the linear equations

S = {n ≥ s : b -s (n) = 0}. Let E = {(u |n| ) n∈ : n ∈ AE\S ⇒ (P • u) n =
(4.1) (P • u) n = 0, n ∈ r, s -1 ∪ S.
These equations are linearly independent. Let t = |S|, and let p ≤ s -r + t denote the rank of (4.1), considered as a system of linear forms over E. For large enough n 0 , sequences u ∈ E are in bijection with their values (u n ) n∈J , where J = (S -s) ∪ n 0 , n 0 + 2s -1 , hence E has dimension 2s + t. It follows that dim S = dim E -p ≥ s + r.

Similarly, an element of E ∩ C is characterized, for large n 0 , by a convergent sequence (u n ) n≥n0 ("a convergent germ of solution") and values u n for n ∈ S. It belongs to S when additionally (4.1) is satisfied. Thus, by Corollary 3.6, we have dim(E ∩ C) = s + t and dim(S ∩ C) = s + t -q, where q ≤ p is the rank of the system (4.1) restricted to E ∩ C. But we already know from Theorem 2.1 that S∩C ≃ ker L, and hence dim S∩C = r. Therefore, we have r = s+t-q ≥ s+t-p ≥ r, and hence p = q = s + t -r.

The important fact is that the equations (4.1) are independent. The other statements are there to complete the picture but are not really used in the sequel. Incidentally, Proposition 4.1 implies that divergent solutions of a Chebyshev recurrence have the same exponential growth near positive and negative infinity.

The full procedure is stated as Algorithm 4.2. As the handling of boundary conditions is naturally incorporated into the algorithm, here, we see the initial conditions y (i) (0) = ℓ i as a special case of boundary conditions

(4.2) λ i (y) = ℓ i , 1 ≤ i ≤ r,
each of the form λ i (y) = q j=1 µ j y (rj) (x j ) with x j ∈ [-1, 1] and r j ≤ r. In general, the boundary conditions are assumed to be chosen so that the function y of interest is the unique solution of (1.1) satisfying (4.2). They are independent in the sense that the linear forms λ i : ker L → are linearly independent.

Motivated by Proposition 4.1, we "unroll" s + t linearly independent test sequences f i ∈ E. We then solve the linear system (4.3), consisting of the constraints (4.1) and of approximations of the boundary conditions (4.2), to select a single linear combination of the f i as output. Algorithm 4.2 takes as input both a target degree d and a starting index N . We will see in the next section how the choice of N influences the quality of the output. In practice, taking N = d + s usually yields good results.

(n)S k associated to L 2 set S = {n ≥ s : b -s (n) = 0} and I = S ∪ N, N + s -1 3 for n from N + 2s -1 down to 1 4 for i ∈ I 5 if n = i then set f i,n-s = 1 6 else if n ∈ I or n ≥ N + s then set f i,n-s = 0 7 else compute f i,n-s using the relation (P • f ) n = 0 8 using indeterminates η i , i ∈ I, set ỹn = i∈I η i f i,|n| , |n| ≤ N ỹn = 0, |n| > N, and ỹ(x) = N n=-N ỹn T n (x)
9 solve for (η i ) i∈I the linear system

(4.3) λ k (ỹ) = ℓ k , 1 ≤ k ≤ r, b -s (n)ỹ n-s + • • • + b s (n)ỹ n+s = 0, n ∈ r, s -1 ∪ S 10 return d n=-d ỹn T n (x)
The complexity of Algorithm 4.2 is easy to estimate. Its correctness is less obvious. At first sight, there could conceivably exist differential equations for which Algorithm 4.2 always fails, no matter how large N is chosen. It could happen for instance that the kernel of (4.1) (a system we know to have full rank over E by Proposition 4.1) always has a nontrivial intersection with Span{t i : i ∈ I}. It is not clear either that, when the algorithm does return a polynomial p, this polynomial is close to y. We prove in the next section that these issues do not occur. But already at this point, we note that if the result happens to be satisfactory, it is already possible to validate it (that is, to get a rigorous good upper-bound on y -p ∞ ) using the methods of Section 6. 4.2. Convergence. We now prove that Algorithm 4.2 converges. The fact that it does not fail for large N will come as a byproduct of the convergence proof. The proof, inspired in part by the analysis of the generalized Miller algorithm [START_REF] Zahar | A mathematical analysis of Miller's algorithm[END_REF][START_REF] Wimp | Computation with Recurrence Relations[END_REF], is based on the asymptotic behaviour of the solutions of the Chebyshev recurrence predicted by Theorem 3.3. The approach of backward recurrence algorithms based on this theorem was pioneered by Wimp [START_REF] Wimp | On recursive computation[END_REF].

Retaining the notation from the previous subsection, assume that the operator L and the boundary conditions λ i (y) = ℓ i are fixed. Write the Chebyshev expansion of y as

y(x) = ∞ n=-∞ y n T n (x).

Let y (N ) n

= ỹn , |n| ≤ N , be the coefficients computed by Algorithm 4.2 (run in exact arithmetic) when called with the starting index N .

The central result of the analysis of Algorithm 4.2 is the following theorem. It implies that when d is fixed and N → ∞, the polynomial output by Algorithm 4.2 converges at least exponentially fast to the truncated Chebyshev series π d (y). The base of the exponential is related to the asymptotics of the "slowest decreasing" convergent solution of the Chebyshev recurrence, in turn related to the location of the singular points of the differential equation (1.1).

Theorem 4.4. Algorithm 4.2 fails for finitely many N only. As N → ∞, its output satisfies

N max n=-N |y (N ) n -y n | = O(N τ e 1,N )
for some τ independent of N .

We write f (N ) = O pol (g(N )) when there exists τ ≥ 0 such that f (N ) = O(N τ g(N )) as n → ∞.

Proof. The finite sequence (y s+t to write these last s -r + t equations in the same form as the first r, so that (4.3) becomes

(4.4) λ (N ) i (y (N ) ) = N n=-N λ i,n y (N ) n = ℓ i , 1 ≤ i ≤ s + t.
Now let (e 1 , . . . , e s , e -1 , . . . , e -s ) be a basis of the solutions of P • u = 0 in the neighborhood of +∞ of the form provided by Corollary 3.6. Extend each e i to an element of E, and then the tuple to a basis of E by setting e s+1 = f n1 , . . . , e s+t = f nt where n

1 < n 2 < • • • < n t are the elements of S. Let (4.5) ∆ (N ) = e 1,N • • • e s+t,N e -1,N • • • e -s,N . . . . . . . . . . . . e 1,N +s-1 • • • e s+t,N +s-1 e -1,N +s-1 • • • e -s,N +s-1 λ (N ) 1 (e 1 ) • • • λ (N )
1 (e s+t ) λ

y (N ) = s+t k=-s γ (N ) k e k , γ (N ) k = ∆ (N ) k ∆ (N ) .
Algorithm 4.2 fails if and only if ∆ (N ) = 0. The sequence of "exact" Chebyshev coefficients of the function y defined by the input is likewise given by (4.7)

y = s+t k=1 γ k e k , γ k = ∆ k ∆ ,
where ∆ = det(λ i (e j )) 1≤i,j≤s+t and ∆ j denotes the determinant ∆ with the j-th column replaced by (ℓ 1 , . . . , ℓ s+t ) T .

Our goal is now to prove that γ

(N ) k → γ k fast as N → ∞.
To do that, we study the asymptotic behaviours of the determinants ∆ (N ) and ∆ (N ) k . We decompose ∆ (N ) into the four blocks indicated by Eq. (4.5) as follows:

∆ (N ) = A B C D .
The corresponding modified blocks in ∆

(N ) k are denoted A k , B k , C k , D k .
(We drop the explicit index for readability, but notice that these matrices depend on N .) The blocks B and C are nonsingular for large N , the first one by Lemma 3.4 and the second one because det C → ∆ = 0 as N → ∞. The Schur complement formula implies

∆ (N ) = -det(B) det(C) det(I -C -1 DB -1 A).
Setting e j = (e j,N , . . . , e j,N +s-1 ) T , the entry at position (i, j) in the matrix B -1 A satisfies (B -1 A) i,j = 0 for large N if j > s, and otherwise (where the notation • indicates the omission of the corresponding term) as N → ∞ by Lemma 3.4.

(B -1 A) i,j = det(e -
In view of our assumptions on the boundary conditions (4.2), we have λ i,n = O n→±∞ (n r ) for all i ≤ σ, where σ is the maximum derivation order appearing in (4.2). Additionally, the sequences λ i,n with r + 1 ≤ i ≤ s + t are ultimately zero. Therefore the entries of D satisfy

D i,j = λ (N ) i (e -j ) = O pol (e -j,N ).

This yields the estimate (DB

-1 A) i,j = O pol (e j,N ) for the j-th column of DB -1 A. Since C i,j = λ (N ) i (e j ) = λ i (e j ) + O pol (e j,N ),
we get (C -1 DB -1 A) i,j = O pol (e j,N ) as well, and

∆ (N ) = -det(B) det(C)(1 -tr(C -1 DB -1 A) + O( C -1 DB -1 A 2 )) = -det(B)(∆ + O pol (e 1,N )).
In particular, ∆ (N ) = 0 for all large enough N , hence, for any fixed differential equation, the algorithm fails at most for finitely many N .

We turn to the modified determinants

∆ (N ) k = A k B k C k D k .
For k > 0, the same reasoning as above (except that C k may now be singular) leads to

∆ (N ) k = -det(B) det(C k -DB -1 A k ) = -det(B)(det(C k ) + O pol (e 1,N )) = -det(B)(∆ k + O pol (e 1,N )), hence (4.8) γ (N ) k = ∆ (N ) k ∆ (N ) = γ k + O pol (e 1,N ), k > 0.
In the case k < 0, write

∆ (N ) k = -det(C) det(B k -AC -1 D k ).
The natural entrywise bounds on A and D yield (C -1 D k ) i,j = O(N r e -j,N +s-1 ) and from there as N → ∞, uniformly in n.

(AC -1 D k ) i,j = O(N r e 1,N e -j,N +s-1 ) = o(e -j,N ), so that (B k -AC -1 D k ) i,j ∼ e -j,N +i-1 , j = -k. For j = -k however, the j-th column of B k is zero and that of D k is constant, hence (B k -AC -1 D k ) i,j = O(e 1,N ), j = -k. It follows that det(B k + AC -1 D k ) = O pol (e -1,N • • • e k,N • • • e -s,
Remark 4.5. In the special case where the solution y is a polynomial, it is computed exactly.

Theorem 4.4 implies that the polynomial y (N ) = ỹ returned by Algorithm 4.2 satisfies (4.10)

y (N ) -π d (y) ∞ ≤ φ(N )N ! κ1 α N 1 ,
where κ 1 and α 1 are the asymptotic growth parameters defined in Section 3.2, for some φ with ln φ(N ) = o(N ). Thus, given ǫ > 0, it suffices to take N = O(ln(ǫ -1 )) in order to obtain y (N ) -π d (y) ∞ ≤ ǫ. The constant hidden in the O(•) depends on y. The estimate goes down to O ln(ǫ -1 )/ ln ln(ǫ -1 ) when κ 1 < 0, that is (by a similar argument as in the proof of Proposition 3.5), when the leading coefficient a r of the differential equation is a constant.

Comparing with Equation (2.5), we can state the following "effective near-minimax approximation" property.

max k=0 |y n+k | ≥ n! κ1 |α 1 | n ψ(n), ψ(n) = e o(n) ,
we see by comparing with (4.10) that choosing

N = d + o(d) is enough to get y (N ) -π d (y) ∞ ≤ p * d -y ∞ .

This last inequality in turn implies

y (N ) -y ∞ ≤ 4π -2 ln(d + 1) + 5 p * d -y ∞ .
When, in (4.12), κ 1 and α 1 are replaced by κ i and α i for some i > 1, the estimate N = d + o(d) still holds in the case where

κ i = κ 1 < 0. It becomes N = O(d) when either κ i = κ 1 = 0 or κ i < κ 1 < 0, and N = O(d ln d) in general (that is, when κ i < κ 1 = 0). Remark 4.7.
Assuming only e i,n+1 /e i,n ∼ α i n κi e i,n instead of full asymptotic expansions in Lemma 3.4, but with the additional hypothesis κ i = κ j ⇒ α i = α j , one can prove that

C(n) ∼ n→∞ e 0,n e 1,n+1 • • • e s-1,n+s-1 i<j κi=κj α i α j -1 ,
if the e i are sorted so that

κ 0 ≤ κ 1 ≤ • • • ≤ κ s-1 .
This leads to a weaker variant of Theorem 4.4 that does not rely on the Birkhoff-Trjitzinsky theorem.

Examples.

We have developed a prototype implementation of Algorithm 4.2 in Maple [38]. Our implementation uses exact rational arithmetic for operations on the coefficients of approximation polynomials. The experimental source code can be downloaded from http://homepages.laas.fr/mmjoldes/Unifapprox Besides Algorithm 4.2, it includes a (not entirely rigorous with respect to several minor outwards rounding issues) proof-of-concept implementation of the validation algorithm of Section 6 further discussed in that section. The gfsRecurrence package [START_REF] Benoit | Algorithmique semi-numérique rapide des séries de Tchebychev[END_REF] available on the same web page provides tools to compute Chebyshev recurrences from linear differential equations as discussed in Section 2.2.

For each of the following examples, Figure 2 shows the graph of the difference between the polynomial approximation of a given degree computed by the implementation and the known exact solution, illustrating the quality of the approximations. (ii) Next, we consider the fourth order initial value problem (taken from Geddes [24, p. 31])

y (4) (x) -y(x) = 0, y(0) = -y ′′ (0) = 3 2 , -y ′ (0) = y ′′′ (0) = 1 2 ,
with the exact solution

y(x) = 3 2 cos(x) - 1 2 sin(x).
(iii) Finally, the second-order differential equation (2x 2 + 1)y ′′ (x) + 8xy ′ (x) + (2x 2 + 5)y(x), y(0) = 1, y ′ (0) = 0, has complex singular points at z = ±i/ √ 2, relatively close to the interval [-1, 1], and admits the exact solution

y(x) = cos(x) 2x 2 + 1 .
On our test system, using Maple 17, the total computation time for each example is of the order of 0.05 to 0.1 s.

According to a classical theorem of de la Vallée Poussin [12, Section 3.4], the nearuniform amplitude of the oscillation observed in the first two examples indicates an approximation error very close to that of the minimax approximation. Table 2 in Section 6.2 (p. [START_REF] Lewanowicz | Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series[END_REF] gives numerical values of p -y ∞ and p * -y ∞ in each case. We will later extend these examples to include in the comparison the bounds on p -y ∞ output by the validation algorithm.

4.4.

A Link with the Tau Method. Besides Clenshaw's, another popular method for the approximate computation of Chebyshev expansions is Lánczos' tau method [START_REF] Lanczos | Trigonometric interpolation of empirical and analytical functions[END_REF][START_REF] Lanczos | Applied analysis[END_REF]. It has been observed by Fox [START_REF] Fox | Chebyshev methods for ordinary differential equations[END_REF] and later in greater generality (and different language) by El Daou, Ortiz and Samara [START_REF] El-Daou | A unified approach to the tau method and Chebyshev series expansion techniques[END_REF] that both methods are in fact equivalent, in the sense that they may be cast into a common framework and tweaked to give exactly the same result. We now outline how the use of the Chebyshev recurrence fits into the picture. This sheds another light on Algorithm 4.2 and indicates how the Chebyshev recurrence may be used in the context of the tau method. As in the previous sections, consider a differential equation L • y = 0 of order r, with polynomial coefficients, to some solution of which a polynomial approximation of degree d is sought. Assume for simplicity that there are no nontrivial polynomial solutions, i.e., (ker L) ∩ [x] = {0}.

In a nutshell, the tau method works as follows. The first step is to compute L • p where p is a polynomial of degree d with indeterminate coefficients. Since (ker L) ∩

[x] = {0}, the result has degree greater than d. One then introduces additional unknowns τ d+1 , . . . , τ d+m in such number that the system (4.13)

L • p = τ d+1 T d+1 + • • • + τ d+m T d+m λ i (p) = ℓ i (1 ≤ i ≤ r)
• f i = 0) or |n| > N + s (since ỹn = 0 for |n| > N ), so that Supp(Q • v) ⊂ N -s, N + s -1 .
It can be checked that the Chebyshev recurrence associated to L = ( d dx ) r is P = δ r (n) : indeed, in the language of [START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF], it must be the first element of a pair (P 1 , Q 1 ) satisfying

Q -1 1 P 1 = I -r . Thus δ r (n) • u = Q • v is equivalent to u (r) = v, whence (4.14) v(x) = d r dx r |n|>r (P • ỹ) n δ r (n) T n (x) = N -s≤|n|<N +s (P • ỹ) n δ r (n) T (r) n (x).
We see that the output ỹ(x) of Algorithm 4.2 satisfies an inhomogeneous differential equation of the form

L • ỹ = τ N -s T (r) N -s (x) + • • • + τ N +s-1 T (r)
N +s-1 (x). (However, the support of the sequence (v n ) itself is not sparse in general.)

This point of view also leads us to the following observation. In comparison, the best known arithmetic complexity bound for the conversion of arbitrary polynomials of degree d from the Chebyshev basis to the monomial basis is O(M (d)), where M stands for the cost of polynomial multiplication [START_REF] Pan | New fast algorithms for polynomial interpolation and evaluation on the Chebyshev node set[END_REF][START_REF] Bostan | Power series composition and change of basis[END_REF].

Proof. As already mentioned, the Taylor series expansion of a function that satisfies an LODE with polynomial coefficients obeys a linear recurrence relation with polynomial coefficients. In the case of an inhomogeneous equation L • u = v, the recurrence operator does not depend on v, and the right-hand side of the recurrence is the coefficient sequence of v. Now ỹ satisfies L • ỹ = v where v is given by (4.14). The coefficients (P • ỹ) n /δ r (n) of (4.14) are easy to compute from the last few Chebyshev coefficients of ỹ. One deduces the coefficients v n in linear time by applying repeatedly the non-homogeneous recurrence relation

(4.15) T ′ n-1 (x) = -T ′ n+1 (x) + 2xT ′ n (x) + 2T n (x)
obtained by differentiation of the equation (2.6), and finally those of the expansion of ỹ on the monomial basis using the recurrence relation they satisfy.

Chebyshev Expansions of Rational Functions

This section is devoted to the same problems as the rest of the article, only restricted to the case where y(x) is a rational function. We are interested in computing a recurrence relation on the coefficients y n of the Chebyshev expansion of a function y, using this recurrence to obtain a good uniform polynomial approximation of y(x) on [-1, 1], and certifying the accuracy of this approximation. All this will be useful in the validation part of our main algorithm.

Our primary tool is the change of variable x = 1 2 (z + z -1 ) followed by partial fraction decomposition. Similar ideas have been used in the past with goals only slightly different from ours, like the computation of y n in closed form [START_REF] Einwohner | A MACSYMA package for the generation and manipulation of Chebyshev series[END_REF][START_REF] Mathar | Chebyshev series expansion of inverse polynomials[END_REF]. Indeed, the sequence (y n ) n∈AE turns out to obey a recurrence with constant coefficients. Finding this recurrence or a closed form of y n are essentially equivalent problems. However, we will use results regarding the cost of the algorithms that do not seem to appear in the literature. Our main concern in this respect is to avoid conversions of polynomial and series from the monomial to the Chebyshev basis and back. We also need simple error bounds on the approximation of a rational function by its Chebyshev expansion.

(S + S -1 )) • (y n ) = (a n ).
Proof. This is actually the limit case r = 0 of Theorem 2.1, but a direct proof is very easy: just write

deg b i=-deg b b i z i ∞ n=-∞ y n z n = ∞ n=-∞ ∞ i=-∞ b i y n-i z n = ∞ n=-∞ a n z n , x = z + z -1 2 ,
and identify the coefficients of like powers of z.

As in the general case (Section 2), this recurrence has spurious (divergent) solutions besides the ones we are interested in. However, we can explicitly separate the positive powers of z from the negative ones in the Laurent series expansion

(5.1) ŷ(z) = y z + z -1 2 = ∞ n=-∞ y n z n , ρ -1 < |z| < ρ,
using partial fraction decomposition. From the computational point of view, it is better to start with the squarefree factorization of the denominator of ŷ:

(5.2)

β(z) = z deg b b z + z -1 2 = β 1 (z)β 2 (z) 2 • • • β k (z) k
and write the full partial fraction decomposition of ŷ(z) in the form

(5.3) ŷ(z) = q(z) + k i=1 βi(ζ)=0 i j=1 h i,j (ζ) (ζ -z) j , q(z) = n q n z n ∈ É[z],
where h i,j ∈ É(z). The h i,j may be computed efficiently using the Bronstein-Salvy algorithm [START_REF] Bronstein | Full partial fraction decomposition of rational functions[END_REF] (see also [START_REF] Gourdon | Effective asymptotics of linear recurrences with rational coefficients[END_REF]).

We obtain an identity of the form (5.1) by expanding the partial fractions corresponding to poles ζ with |ζ| > 1 in power series about the origin, and those with |ζ| < 1 about infinity. The expansion at infinity of

h i,j (ζ) (ζ -z) j = (-1) j z -j h i,j (ζ) (1 -ζz -1 ) j
converges for |z| > |ζ| and does not contribute to the coefficients of z n , n ≥ 0 in the complete Laurent series. It follows from the uniqueness of the Laurent expansion of ŷ on the annulus ρ -1 < |z| < ρ that 5

(5.4)

∞ n=0 y n z n = q(z) + k i=1 βi(ζ) = 0 |ζ| > 1 i j=1 h i,j (ζ) (ζ -z) j .
We now extract the coefficient of z n in (5.4) and use the symmetry of (y n ) n∈ to get an explicit expression of y n in terms of the roots of b( 1 2 (z + z -1 )).

Proposition 5.2. The coefficients of the Chebyshev expansion y(x) = n y |n| T n (x) are given by (5.5)

y n = q n + k i=1 i j=1 βi(ζ) = 0 |ζ| > 1 n + j -1 j -1 h i,j (ζ)ζ -n-j (n ≥ 0)
where the q n ∈ É, β i ∈ É[z] and h i,j ∈ É(z) are defined in Equations ( 5.2) and (5.3).

Note that (5.4) also yields a recurrence of order deg b on (y n ) n∈AE , instead of 2 deg b for that from Proposition 5.1, but now with algebraic instead of rational coefficients in general.

Truncation Error.

We can now explicitly bound the error in truncating the Chebyshev expansion of y. Proposition 5.3. Let y ∈ É(x) have no pole within the elliptic disk E ρ (see (2.4)).

Assume again the notations from (5.2) and (5.3). For all d ≥ deg q, it holds that

n>d y n T n ∞ ≤ k i=1 i j=1 βi(ζ) = 0 |ζ| > 1 |h i,j (ζ)|(d + 2) j-1 (|ζ| -1) j |ζ| -d-1 = O(d deg b ρ -d ) as d → ∞. Proof. We have n>d y n T n ∞ ≤ n>d |y n | because T n ∞ ≤ 1 for all n. Using the inequality n>d n + j -1 j -1 t n+j ≤ (d + 2) j-1 t d+1 ∞ n=0 n + j -1 j -1 t n+j = (d + 2) j-1 t d+j+1 (1 -t) j 5
To prevent confusion, it may be worth pointing out that in the expression ŷ(z) = q(z)

+ k i=1 βi(ζ) = 0 |ζ| > 1 i j=1 ( h i,j (ζ) (ζ -z) j + h i,j (ζ -1 ) (ζ -1 -z) j )
the Laurent expansion of a single term of the form

h i,j (ζ) (ζ-z) j + h i,j (ζ -1 ) (ζ -1
-z) j is not symmetric for j > 1, even if q(z) = 0.

for t < 1, the explicit expression from Proposition 5.2 yields

n>d |y n | ≤ n>d k i=1 i j=1 βi(ζ)=0 |ζ|>1 n + j -1 j -1 |h i,j (ζ)||ζ| -n-j ≤ k i=1 i j=1 βi(ζ)=0 |ζ|>1 |h i,j (ζ)|(d + 2) j-1 (|ζ| -1) j |ζ| -d-1 . Since |ζ| > 1 actually implies |ζ| > ρ when b( 1 2 (ζ + ζ -1
)) = 0, the asymptotic estimate follows.

Computation.

There remains to check that the previous results really translate into a linear time algorithm. We first state two lemmas regarding basic operations with polynomials written on the Chebyshev basis. Proof. It suffices to loop over the indices (i, j) and, at each step, add to the coefficients of T |i±j| in the product the contribution coming from the coefficient of T i in a and that of T j in b, according to the formula 2T i T j = T i+j + T i-j .

The naïve Euclidean division algorithm [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Algorithm 2.5] runs in linear time with respect to the degree of the dividend when the divisor is fixed. Its input and output are usually represented by their coefficients in the monomial basis, but the algorithm is easily adapted to work in other polynomial bases. We end up with Algorithm 5.6. In view of future needs, it takes as input a polynomial of arbitrary degree already written in the Chebyshev basis in addition to the rational function (of bounded degree) y. The details of the algorithm are only intended to support the complexity estimates, and many improvements are possible in practice.

β i (ζ) = 0 ∧ |ζ| > 1 ⇒ 1 < ρ -≤ |ζ| ≤ ρ + 7 compute M ≥ k i=1 i j=1 j(deg β i ) sup ρ-≤|ζ|≤ρ+ |h ′ i,j (ζ)| + |ζ -1 h i,j (ζ)| ρ -j - 8 set ǫ ′ := min ρ --1, M -1 1 -ρ -1 - D+1 ǫ 4 , with D = deg b 9 compute approximations ζ ∈ É[i] of the roots ζ of β i such that ζ -ζ < ǫ ′ 10 for 0 ≤ n ≤ d ′ 11 set ỹn = q n + Re k i=1 i j=1 βi(ζ)=0 |ζ|>1 n + j -1 j -1 h i,j ( ζ) ζ-n-j 12 return ỹ(x) = d ′ n=-d ′ ỹn T n (x)
Proof. Firstly, we prove the error bound ỹ -y ∞ ≤ ǫ.

Let A = {ζ : ρ -≤ |ζ| ≤ ρ + } and M 0 = sup ζ∈A |h ′ i,j (ζ)|, M 1 = sup ζ∈A |ζ -1 h i,j (ζ)|.
On A, we have

d dζ h i,j (ζ) ζ -n-j ≤ (M 0 + (n + j)M 1 ) |ζ| -n-j ≤ (n + j)(M 0 + M 1 )ρ -n-j - .
By Proposition 5.2, observing that the condition |ζ -ζ| < ρ --1 from Step 8 implies [ζ, ζ] ⊂ A, we have the inequalities

|y n -ỹn | ≤ k i=1 i j=1 βi(ζ)=0 |ζ|>1 n + j -1 j -1 h i,j (ζ)ζ -n-j -h i,j ( ζ) ζ-n-j ≤ k i=1 i j=1 j (deg β i ) n + j j (M 0 + M 1 )ρ -n-j - ǫ ′ ≤ M n + D D ρ -n -ǫ ′ .
Therefore, the output of Algorithm 5.6 satisfies

y n -ỹn ∞ ≤ d ′ n=-d ′ |y n -ỹn | + 2 n>d ′ y n T n ≤ 2M ǫ ′ (1 -ρ -1 -) D+1 + 2 ǫ 4 ≤ ǫ.
By Proposition 5.3, for all ǫ > 0, there exists d ′ ≤ K ln(ǫ -1 ) with n>d y n T n ≤ ǫ. In addition, we have deg q ≤ d, hence the degree of the output, as computed in Step 5, satisfies d ′ ≤ max(d, K ln(ǫ -1 )).

Turning to the complexity analysis, Steps 1 and 4-8 have constant cost. So does each iteration of the final loop, assuming the powers of ζ-1 are computed incrementally. Steps 2 and 3 take O(d) operations by Lemmas 5.4 and 5.5, and do not depend on ǫ. Regarding Step 9, it is known [START_REF] Pan | Optimal and nearly optimal algorithms for approximating polynomial zeros[END_REF]Theorem 1.1(d)] that the roots of a polynomial with integer coefficients can be approximated with of absolute accuracy η in O(η -1 ) arithmetic operations. (In fact, the bit complexity of the algorithm beyond this statement is also softly linear in ln(ǫ -1 ).) Since M depends neither on ǫ nor on d, we have ǫ ′ = Ω(ǫ), and hence the cost of Step 9 is in O(ln(ǫ -1 )).

Validation

Assume that we have computed a polynomial of degree d

p(x) = d n=-d ỹn T n (x)
which presumably is a good approximation on [-1, 1] of the D-finite function y defined by (1.1,1.2). As stated in Problem 1.1, our goal is now to obtain a reasonably tight bound B such that y -p ∞ ≤ B.

6.1. Principle. The main idea to compute the bound is to convert the initial value problem defining y into a fixed-point equation T (y) = y, verify explicitly that T maps some neighborhood of p into itself, and conclude that p must be close to the "true" solution. This is one simple instance of the functional enclosure methods mentioned in Section 1.1 and widely used in interval analysis [START_REF] Moore | Methods and applications of interval analysis[END_REF]. In most cases (with the notable exception of [START_REF] Epstein | Ultra-arithmetic II: intervals of polynomials[END_REF][START_REF] Kaucher | Self-validating numerics for function space problems[END_REF][START_REF] Kaucher | Validating computation in a function space[END_REF]), these methods are based on interval Taylor series expansions such as so-called Taylor models [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF][START_REF] Neumaier | Taylor forms -Use and limits[END_REF]. Here, we describe an adaptation designed to work in linear time with polynomial approximations written on the Chebyshev basis.

The first step is to reduce the differential initial value problem defining y to a linear integral equation of the Volterra type and second kind. Let α 0 , α 1 , . . . , α r ∈ É[x] be such that

L = a r ∂ r + • • • + a 1 ∂ + a 0 = ∂ r α r + • • • + ∂α 1 + α 0 ,
and, for k ∈ 0, r , define

L k = ∂ r-k α r + • • • + ∂α k+1 + α k ∈ É[x] ∂ .
Observe that for all a ∈ É[x] and k ∈ AE, the operator ∂ k a is of the form

∂ k a = a∂ k + * ∂ k-1 + • • • , so that α r = a r . In particular, α r does not vanish on [-1, 1]. Lemma 6.1. The solution y ∈ X of L•y = 0 such that y (k) (0) = ℓ k for k ∈ 0, r-1 satisfies (6.1) α r (x)y(x) = g(x) + x 0 K(x, t)y(t)dt, where (6.2) K(x, t) = - r-1 k=0 (x -t) k k! α r-1-k (t) = r-1 j=0 β j (t)x j , β j (t) = r-1-j i=0 (-1) i+1 i!j! t i α r-1-j-i (t), and
(6.3) g(x) = r-1 k=0 (L r-k • y)(0) x k k! = r-1 k=0 k j=0 k i=j i j α (i-j) r-k+i (0)ℓ j x k k! .
Proof. We have

x 0 L k • y = (L k+1 • y)(x) -(L k+1 • y)(0) + x 0 α k y.
Applying this formula to integrate r times the equation L

• y = 0 yields (α r y)(x) + x 0 (α r-1 y)(x 1 )dx 1 + • • • + x 0 x1 0 • • • xr-1 0 (α 0 y)(x r )dx r • • •dx 2 dx 1 = (L 1 •y)(0) x r-1 (r -1)! + • • • + (L r-1 •y)(0)x + (L r •y)(0),
and (6.1) follows using the relation

x 0 x1 0 • • • x k-1 0 f (x k )dx k • • • dx 2 dx 1 = x 0 (x -t) k-1 (k -1)! f (t)dt
(which can be obtained by repeated integrations by parts).

Let X denote the Banach space of continuous functions from [0, 1] to Ê (or ), equipped with the uniform norm. With K and g as in Lemma 6.1, define an operator T : X → X by (6.4)

T (f )(x) = α r (x) -1 g(x) + x 0 K(x, t)f (t)
-y ∞ ≤ γ i T i (p) -p ∞ , γ i = ∞ j=0 A ij (ij)! , holds for all i ≥ 1. Proof. Let V : f → α -1 r
Kf denote the linear part of T . Equation (6.4) rewrites again as (Id -V ) • y = g. More generally, we have

(Id -V i ) • y = T i (y) -V i • y = T i (f ) -V i • f = (Id +V + • • • + V i-1 ) • g for all f ∈ X.
The operator V is continuous, and since

|V i • y(x)| ≤ x 0 A x1 0 A • • • xi-1 0 A y ∞ dx i • • • dx 2 dx 1 = A i i! y ∞ , the (subordinate) norm of V i satisfies V i ≤ A i i! .
Therefore, the series j V ij converges, (Id -V i ) is invertible, and (Id

-V i ) -1 ≤ γ i . Writing p -y = (Id -V i ) -1 • [(Id -V i ) • p -(T i (p) -V i • p)] = (Id -V i ) -1 • (p -T i (p))
yields the announced result.

The interesting fact about this bound is that we can effectively compute rigorous approximations of the right-hand side. Indeed, p and g are explicit polynomials, so that it is not too hard to compute T i (p) -p approximately while keeping track of the errors we commit, and deduce a bound on its norm.

Choosing i = 1 (and γ 1 = e A ) in Proposition 6.2 already yields a nontrivial estimate, but (as we shall see in more detail) larger values of i are useful since γ i → 1 as i → ∞. In particular, we have (6.5)

γ i ≤ 1 1 -A i /i!
for large i. This crude bound will be enough for our purposes. In practice, though, it is better to compute an approximation of γ i that is closer to reality (e.g., using the first few terms of the series and a bound on the tail) in order to reduce the required number of iterations of V .

f ∞ ≤ M ≤ √ d + 1 f ∞ in O(d) arithmetic operations. Proof. It suffices to take M = d n=-d |c n |. Indeed, we have f 2 = 1 π 1 -1 f (t) 2 √ 1 -t 2 dt 1/2 = |c 0 | 2 + 4 d n=1 |c n | 2 . It follows that f 2 ≤ f ∞ ≤ M ≤ √ d + 1 f 2 ,
where the first inequality results from the integral expression of f 2 , the second one, from the fact that T n ∞ ≤ 1 for all n, and the last one from the Cauchy-Schwarz inequality.

This results in Algorithm 6.5. Again here, the suggested bounds for α -1 r ∞ , p -p i ∞ and γ i are only intended to support the complexity estimate, and tighter choices are possible in practice. Algorithm 6.5. Input:

A differential operator L = a r ∂ r + • • •+ a 1 ∂ + a 0 ∈ É[x] ∂ of order r such that a r (x) = 0 for x ∈ [-1, 1]
, initial values ℓ 0 , ℓ 1 , . . . , ℓ r-1 . A degree-d polynomial p written on the Chebyshev basis. An accuracy parameter ǫ > 0. Output: A real number B > 0 such that y -p ∞ ≤ B, where y is the unique solution of L • y = 0 satisfying y(0) = ℓ 0 , y ′ (0) = ℓ 1 , . . . , y (r-1) (0) = ℓ r-1 .

using the commutation rule x∂

= ∂x -1, compute polynomials α 0 , α 1 , . . . , α r ∈ É[x] such that L = ∂ r α r + • • • + ∂α 1 + α 0 2 define K ∈ É[t, x], (β j ) r-1 j=0 ∈ É[x] r and g ∈ É[x]
as in Lemma 6.1

3 compute A ≥ max |α r (x) -1 K(x, t)| : 0 ≤ t ≤ x ≤ 1 or -1 ≤ x ≤ t ≤ 0
(e.g., using Algorithm 5.6 to expand α r (x) -1 in Chebyshev series), and define (γ i ) ∞ i=1 as in Proposition 6.2

4 compute the minimum i such that A i /i! ≤ 1/2 5 set p 0 = p 6 for k = 0, 1, . . . , i -1 7 compute q k+1 (x) = g(x) + r-1 j=0 x j x 0 β j (t)p k (t)dt ∈ É[x] 8 compute p k+1 ∈ É[x] such that p k+1 -α r (x) -1 q k+1 (x) ∞ ≤ ǫ using Algo- rithm 5.6 9 compute δ ≥ p -p i ∞ using Lemma 6.4 10 return B = γ i (δ + e A ǫ)
We now prove that Algorithm 6.5 works as stated, and estimate how tight the bound it returns is. Theorem 6.6. Algorithm 6.5 is correct: its output B is an upper bound for yp ∞ . For fixed L, as d → ∞ and ǫ → 0, the bound B satisfies

B = O ( y -p ∞ + ǫ)(d + ln(ǫ -1 )) 1/2
and the algorithm performs O(d + ln(ǫ -1 )) arithmetic operations.

Proof. Denote by V the linear part of T and recall from the proof of Proposition 6.2 that V i ≤ A i /i!. For all k, the polynomial p k+1 computed on line 8 satisfies

p k+1 -T (p k ) ∞ ≤ ǫ, hence we have p i -T i (p) ∞ ≤ p i -T (p i-1 ) ∞ + T (p i-1 ) -T 2 (p i-2 ) ∞ + • • • + T i-1 (p 1 ) -T i (p 0 ) ∞ ≤ p i -T (p i-1 ) ∞ + V p i-1 -T (p i-2 ) ∞ + • • • + V i-1 p 1 -T (p 0 ) ∞ ≤ e A ǫ.
By Proposition 6.2, it follows that

(6.6) p -y ∞ ≤ γ i p -T i (p) ∞ ≤ γ i ( p -p i ∞ + p i -T i (p) ∞ ) ≤ γ i (δ + e A ǫ).
This establishes the correctness of the algorithm.

We now turn to the tightness statement. Letting D = deg(p-p i )+ 1, Lemma 6.4 implies that (6.7)

δ ≤ √ D p -p i ∞ where (6.8) p -p i ∞ ≤ p -y ∞ + y -T i (p) ∞ + T i (p) -p i ∞ ≤ 1 + V i p -y ∞ + e A ǫ ≤ 3 2 p -y ∞ + e A ǫ.
Looking at the definition of p k+1 in step

) D ≤ max d, C 2 ln(ǫ -1 ) + C 1 i + 1. (6.9
Plugging (6.8) and (6.9) into (6.7) yields the estimate and the result then follows from the definition of δ since γ i ≤ 2. Finally, the only steps whose cost depends on d or ǫ are lines 7, 8, and 9 (and the number of loop iterations does not depend on these parameters either). By Proposition 5.7, the degrees of p k and q k are all in O(d + ln(ǫ -1 )). The cost of step 7 is linear in this quantity by Lemmas 5.4 and 6.3. The same goes for line 8 by Proposition 5.7, and for line 9 by Lemma 6.4.

Another way to put this is to say that Algorithm 6.5 can be modified to provide an enclosure of p -y ∞ . Indeed, Equations (6.7) and (6.8) imply .

In particular, if we restrict ourselves to polynomials p satisfying p -y ∞ ≤ e -Γd for some fixed Γ, and if ǫ is chosen such that (6.12) then D ≤ max(Γ -1 , 2EC 2 ) ln p -y -1 ∞ + O(1) as p -y ∞ tends to zero, so there exists K (computable as a function of y, Γ, and E) such that (6.13) B ≤ Kb ln(b -1 ) ≤ K p -y ∞ ln( p -y -1 ∞ ) for small p -y ∞ . Of course, since p -y ∞ is what we want to estimate, we do not know the "correct" choice of ǫ beforehand. But, assuming p -y ∞ is indeed small enough, we can search for a suitable ǫ iteratively, starting, say, with ǫ = 2 -d and checking whether (6.13) holds at each step. As our hypotheses imply d = O(ln p -y -1 ∞ ), the whole process requires at most O(ln p -y -1 ∞ ) operations. By combining these tightness guarantees with Corollary 4.6 and lower bounds on p * d -y ∞ such as (4.11), one can devise various strategies to obtain certified polynomial approximations of a given D-finite function y and relate the computed error bounds to p * d -y ∞ .

p -y 2E ∞ ≤ ǫ ≤ p -y E ∞ , E > 1,
Example 6.7. Table 2 gives validated error bounds obtained for the polynomials computed in Section 4.3 using the code presented there. In each case, a naïve implementation of Algorithm 6.5 was called on the polynomial p returned by Algorithm 4.2. (In the third example, the rough bound A suggested in Step 3 was manually replaced by a tighter one to keep the number of iterations small.) The remaining input parameter ǫ was manually set to approximately p -y 2 ∞ based on a heuristic estimate of p -y ∞ . In practice, this makes the term e A ǫ in the error bound (6.6) small, so that the main contribution to the error bound B in practice is p i -p ∞ .

Besides the upper bound B, the table gives a lower bound b ≤ p -y ∞ obtained as discussed above. For comparison, we include the "true" value of p -y ∞ , as well as the error y -p * ∞ corresponding to the minimax polynomial of degre d, computed using Sollya [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF].

The last four columns indicate the values of the parameters D and i and the running time of both algorithms. It can be observed that our choice of ǫ makes D grow significantly larger than d, and that a naïve implementation of Algorithm 6.5, despite its interesting theoretical complexity, is far from being efficient in practice. Nevertherless, for simple examples at least, the total running time remains reasonable. Note for comparison that plotting the error curves shown on Figure 2 is about 1 to 2 times slower than computing the error bounds.

Unfortunately, the above complexity results come short of providing what we may call "validated near-minimax approximations", at least in a straightforward way. More precisely, following Mason and Handscomb [39,Def. 3.2], call an approximation scheme mapping a function y to a polynomial p d of degree at most d near-minimax if it satisfies

p d -y ∞ ≤ Λ(d) p * d -y ∞
where Λ(d) does not depend on y. It is then natural to ask for polynomial approximations where p * d -y ∞ not only satisfies the above inequality, but also comes with an explicit upper bound satisfying a similar inequality, that is (6.14)

p d -y ∞ ≤ B ≤ Λ(d) p * d -y ∞ , Λ independent of y.
We thus leave open the following question. Another subject for future work is the following. In the timespan since we prepared the first draft of this work, an article by Olver and Townsend [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] has appeared that studies a similar question-how to obtain polynomial approximations of solutions of linear ODEs on the Chebyshev basis "in linear time"-from a Numerical Analysis perspective. On first sight at least, the motivations, language, and techniques look quite different from ours, and there appears to be little overlap between the actual results. Yet the methods have common ingredients. Roughly speaking, our algorithm may also be viewed as a coefficient spectral method in the terminology of Olver and Townsend. Their method is more general in the sense that it can deal with non-polynomial coefficients, which also means that they do not directly exploit the Chebyshev recurrence. Instead, the computation of the approximation polynomials (for which we use a block Miller algorithm) boils down to the fast and numerically stable solution of a linear system similar to (4.4). There is no validation of the solution. It is intriguing to understand these links in detail and determine if the best features of the two methods can somehow be combined.

Beyond non-polynomial coefficients, an interesting research direction concerns the case of nonlinear ODEs. We may expect algorithms of a different kind (probably based on Newton's method instead of recurrences) for the computation of polynomial approximations, but some of the ideas used in the present article may still apply. And, closer to what we do here, it is natural to ask for a generalization to other families of orthogonal polynomials, starting with the rest of the class of Gegenbauer polynomials.

Problem 1 . 1 .

 11 Let y : [-1, 1] → R be a D-finite function specified by a linear differential equation with polynomial coefficients and initial conditions. Let d ∈ AE. Given y and d, find the coefficients of a polynomial p(x) = d n=0 c n T n (x) written on the Chebyshev basis (T n ), together with a "small" bound B such that |y(x) -p(x)| ≤ B for all x ∈ [-1, 1].

2. 1 .

 1 Chebyshev Series. Recall that the Chebyshev polynomials of the first kind are polynomials T n (x) ∈ É[x] defined for all n ∈ by the relation T n (cos θ) = cos(nθ). They satisfy T -n = T n for all n. The family (T n ) n∈AE is a sequence of orthogonal polynomials over [-1, 1] with respect to the weight function w

(2. 4 )

 4 E r = {x ∈ : |x + x 2 -1| < r}be the largest elliptic domain with foci in ±1 with this property. Since the singular points are in finite number, we have 1 < r ≤ ∞. The coefficients c n then satisfy c n = O(α n ) for all α > r -1 ; and the Chebyshev expansion (2.2) of f converges uniformly to f on E r[START_REF] Mason | Chebyshev polynomials[END_REF] Theorem 5.16]. Letting x = cos θ and z = e iθ , it is not hard to see that the c n are also the coefficients of the (doubly infinite) Laurent expansion of the function f (z) = f ( z+z -1

(2. 8 )

 8 Sλ = λS (λ ∈ É), Sn = (n + 1)S. Likewise, É[n] S, S -1 ⊂ É(n) S, S -1 is the subring of noncommutative Laurent polynomials in S themselves with polynomial coefficients. The elements of É[n] S, S -1 identify naturally with linear recurrence operators through the left action of É[n] S, S -1 on defined by (n • u) n = nu n and (S • u) n = u n+1 . Recall that L denotes the differential operator appearing in Equation (1.1).

Lemma 2 . 2 .

 22 and let y ∈ C be the Chebyshev coefficient sequence of the (analytic) function L • u. We then have P • u = Q • y by the previous argument. This implies Q • v = Q • y, whence finally y = v by Lemma 2.2 below. The restriction to C of the operator Q from Theorem 2.1 is injective.

  0} be the space of symmetric sequences whose restriction to n ∈ AE satisfies the Chebyshev recurrence, except possibly when n ∈ S. Proposition 4.1.

  (ii) The space S ∩C of symmetric, convergent sequences (u |n| ) n∈ such that P • u = 0 has dimension r. Its elements are the elements of E satisfying (4.1), and the equations (4.1) are independent as linear forms on E ∩ C as well.Proof. First observe that a sequence u ∈ E automatically satisfies (P • u) n = 0 for -n ∈ AE\S too, by Remark 2.6. Then u belongs to S if and only if (P • u) n = 0 for n ∈ S and for |n| < s. But the equations (P • u) n = 0 with |n| < r are trivial by Corollary 2.9, thus u ∈ S if and only if u ∈ E and (4.1) is satisfied.

Algorithm 4 . 2 .

 42 Input: a linear differential operator L of order r, boundary conditions λ 1 (y) = ℓ 1 , . . . , λ r (y) = ℓ r as in (4.2), a target degree d > s, an integer N ≥ max(d, max{n : b -s (n) = 0}). Output: an approximation ỹ(x) = d n=-d ỹn T n (x) of the corresponding solution y of L • y = 0. 1 compute the Chebyshev recurrence operator P = s k=-s b k

Proposition 4 . 3 .

 43 For fixed L, λ i , and ℓ i , Algorithm 4.2 runs in O(N ) arithmetic operations.

n

  ) n∈ of E characterized by the conditions y (N ) N = • • • = y (N ) N +s-1 = 0 from Step 3, along with the linear system (4.3) solved in Step 9.By writing the linear forms λ 1 , . . . , λ r : C → that express the boundary conditions (4.2) as λ i (y) = ∞ n=-∞ λ i,n y n , we define "truncations" λ i,n y n that make sense even for divergent series. (Abusing notation slightly, we apply the λ i and λ (N ) i indifferently to functions, formal Chebyshev series or their coefficient sequences.) The system (4.3) consists of the equations λ (N ) i (y (N ) ) = ℓ i and of the symmetry and extension-through-singularities constraints (4.1). We introduce additional linear forms λ r+1 = λ (N ) r+1 , . . . , λ s+t = λ (N )

  determinant with the column involving e j replaced by (0, . . . , 0 s times , ℓ 1 , . . . , ℓ r , 0, . . . , 0 s -r + t times ) T . By Cramer's rule, provided ∆ (N ) = 0, the sequence (y (N ) n ) n decomposes on the basis (e j ) s+t j=-s of ker P ⊂ as (4.6)

∆

  N e 1,N ) = O pol det((N ) = -det(B) det(C)O pol (e 1,N /e k,N ) -det(B) det(C)(1 + O(e 1,N )) = O pol e 1,N e k,N , k < 0. Combining (4.6), (4.7) with (4.8), (4.9) finally yields y (N ) n = y n + O pol e 1,N s k=1 e k,n + e -k,n e -k,N

Corollary 4 . 6 .

 46 Let L and (ℓ k ) r k=1 be fixed. Given d ∈ N, there exists N such that Algorithm 4.2, called with parameters L, (ℓ k ), d, and N , computes a polynomial p d of degree at most d satisfying p d -y ∞ ≤ 4π -2 ln(d + 1) + 5 p * d -y ∞ in O(ln p d -y -1 ∞ ) arithmetic operations. There is a different way of looking at this, starting with the lower bound [12, Sec. 4.4, Theorem 5(i)] on the quality of the minimax polynomial approximation of degree d of a function y in terms of the Chebyshev coefficients y n of y. In the (typical) case where (4.12) r-1

Figure 2 .

 2 Figure 2. Plot of the error p d -y between a degree-d approximation p d computed by Algorithm 4.2 and the exact solutiony, for each of the problems listed in Section 4.3 and for d ∈ {30, 60, 90}. (i) The first example is adapted from Kaucher and Miranker [30, p. 222]. It concerns the hyperexponential function y(x) = e x/2 √ x + 16 , which can be defined by the differential equation 2 (x + 16)y ′ (x) -(x + 15)y(x) = 0, y(0) = 1 4 .

  has a (preferably unique) solution. The output is the value of p obtained by solving this system; it is an exact solution of the projection π d (L • y) = 0 of the original differential equation. Now let p = d n=-d p n T n and extend the sequence (τ n ) by putting τ n = 0 for |n| ∈ d + 1, d + m and τ -n = τ n . It follows from (4.13) that P • (p n ) = 1 2 Q • (τ n ) where P and Q are the recurrence operators given by Theorem 2.1. Denoting Supp u = {|n| : u n = 0}, we also see from the explicit expression of Q that Supp(Q • τ ) ⊂ d, d + m + 1 . Hence the coefficients p n of the result of the tau method are given by the Chebyshev recurrence, starting from a small number of initial conditions given near the index |n| = d. "Conversely," consider the polynomial ỹ computed in Algorithm 4.2 with N = d, and let v = n v n T n = L • ỹ. We have P • ỹ = Q • v by Theorem 2.1. But the definition of ỹ in the algorithm also implies that (P • ỹ) n = 0 when |n| ≤ N -s (since the ỹn , |n| ≤ N are linear combinations of sequences (f i,n ) |n|≤N recursively computed using the recurrence P

Proposition 4 . 8 .

 48 Assume that Equation 1.1 has no polynomial solution. The expression on the monomial basis of the polynomial ỹ(x) returned by Algorithm 4.2 with N = d can be computed in O(d) arithmetic operations, all other parameters being fixed.

5. 1 .Proposition 5 . 1 .

 151 Recurrence and Explicit Expression. Let y(x) = a(x)/b(x) ∈ É[x] be a rational function with no pole in [-1, 1]. As usual, we denote by (y n ) n∈ , (a n ) n∈ and (b n ) n∈ the symmetric Chebyshev coefficient sequences of y, a and b. The Chebyshev coefficient sequence (y n ) n∈ obeys the recurrence relation with constant coefficients b( 1 2

Lemma 5 . 4 .

 54 The product ab where the operands a, b ∈ É[x] and the result are written in the Chebyshev basis may be computed in O((deg a)(deg b)) operations.

Lemma 5 . 5 .

 55 The division with remainder a = bq + r (deg r < deg b) where a, b, q, r ∈ É[x] are represented in the Chebyshev basis may be performed in O(deg a) operations for fixed b. Proof. Assume n = deg a > deg b = m. The classical polynomial division algorithm mainly relies on the fact that deg(a -b -1 m a n x n-m b) < n where a = i a i x i and b = i b i x i . From the multiplication formula 2T n T m = T n+m + T n-m follows the analogous inequality deg(a -2b -1 m a n T n-m b) < n where a k , b k now denote the coefficients of a and b in the Chebyshev basis. Performing the whole computation in that basis amounts to replacing each of the assignments a ← a -b -1 m a n x n-m b repeatedly done by the classical algorithm by a ← a -2b -1 m a n T n-m b. Since the polynomial T n-m b has at most 2m nonzero coefficients, each of these steps takes constant time with respect to n. We do at most n -m such assignments, hence the overall complexity is O(n).

Proposition 5 . 7 .Algorithm 5 . 6 .

 5756 Algorithm 5.6 is correct. As d → ∞ and ǫ → 0 with all other parameters fixed, it runs in O(d + ln(ǫ -1 )) arithmetic operations and returns a polynomial of degree d ′ ≤ max(d, K ln(ǫ -1 )), where K depends on y, but not on f or d. Input: a rational fraction y(x) = a(x)/b(x), the Chebyshev coefficients of a polynomial f = d n=-d f n T n (x), an error bound ǫ. Output: the Chebyshev coefficients of an approximation ỹ(x) of f y such that ỹ -f y ∞ ≤ ǫ.

1 3 6

 3 convert a and b to Chebyshev basis 2 compute the polynomial g = af , working in the Chebyshev basis 3 compute the quotient q and the remainder r in the Euclidean division of g by b 4 compute the partial fraction decomposition of ŵ(z) = w(x) = r(x)/b(x), where x = 1 2 (z + z -1 ), using the Bronstein-Salvy algorithm5 find d ′ ≥ deg q such that n>d ′y n T n ≤ ǫ/4 using Proposition 5.compute ρ -and ρ + such that

6. 2 .Lemma 6 . 3 .Lemma 6 . 4 .

 26364 Algorithm. The actual computation of T i (p) -p ∞ relies on the following lemmas. One can compute an antiderivative of a polynomial 6 of degree at most d written on the Chebyshev basis in O(d) arithmetic operations.Proof. If f = n c n T n and f ′ = n c ′ n T n , then we have 2nc n = c ′ n-1 + c ′ n+1 according to Equation (2.7). Let f (x) = d n=-d c n T n (x) (with c n = c -n) be a polynomial of degree d, given on the Chebyshev basis. Then one can compute M ≥ 0 such that

  y ∞ + e A ǫ = O ( y -p + ǫ)(d + ln(ǫ -1 )) 1/2

  ǫ , and this b is a computable lower bound for p -y ∞ . Furthermore, using (6.6), we haveδ ≥ γ -1 i p -y -e A ǫ, and hence b ≥ 1 3 √ D p -y ∞ -( √ D + 2)e A ǫ .Comparing with the upper bound on B resulting from (6.10), we deduceB b ln(b -1 ) ≤ 9D p -y ∞ + 6(D + √ D)e A ǫ p -y ∞ -( √ D + 2)e A ǫ ln( p -y -1 ∞ )

enclosure of y -

  104 [2.3 • 10 -53 , 4.3 • 10 -52 ] 3.4 • 10 -52 3.4 • 10 -52 102 2 0.05 0.54 60 10 -194 [9.0 • 10 -99 , 2.4 • 10 -97 ] 2.0 • 10 -97 1.9 • 10 -97 192 2 0.05 1.07 90 10 -284 [4.6 • 10 -144 , 1.5 • 10 -142 ] 1.2 • 10 -142 1.1 • 10 -142 282 2 0.06 1.87 (ii) 30 10 -88 [6.0 • 10 -45 , 9.8 • 10 -44 ] 5.9 • 10 -44 5.6 • 10 -44 42 3 0.06 0.06 60 10 -206 [6.7 • 10 -104 , 1.5 • 10 -102 ] 8.8 • 10 -103 8.5 • 10 -103 72 3 0.07 0.10 90 10 -334 [2.0 • 10 -169 , 5.1 • 10 -168 ] 3.1 • 10 -168 3.0 • 10 -168 102 3 0.08 0.23 (iii) 30 10 -18 [1.2 • 10 -10 , 2.4 • 10 -9 ] 1.6 • 10 -9

Question 6 . 8 .

 68 Given a D-finite function y and a degree bound d, what is the complexity of computing a pair (p d , B) with deg p d ≤ d satisfying (6.14) for some Λ(d)? For instance, can it be done in O(d) + ln( p * d -y -1 ∞ ) arithmetic operations when y is fixed?

Table 1 .

 1 Notation.

	3

  The next proposition is an instance of a classical bound (compare, e.g.,Rall [51, Chap. 1]). As we recall below, when p is close to the fixed point y, its iterated image T i (p) remains close to y, which yields good upper bounds on the distance y -p ∞ .

dt , so that (6.1) becomes T (y) = y. Proposition 6.2. Assume that A is an upper bound on |α r (x) -1 K(x, t)| for x ∈ [-1, 1] and t between 0 and x. Then the bound p

  7, we see that deg q k+1 ≤ deg p k + C 1 for some C 1 > 0 depending on L only. Additionally, according to Proposition 5.7, there exists C 2 (again depending on L only) such that deg p k+1 ≤ max deg q k , C 2 ln(ǫ -1 ) . It follows by induction that deg p

k ≤ max d, C 2 ln(ǫ -1 ) + C 1 k for all k, whence

Table 2 .

 2 1.1 • 10 -9 79 3 0.05 0.74 60 10 -36 [2.2 • 10 -19 , 6.1 • 10 -18 ] 4.1 • 10 -18 3.0 • 10 -18 151 3 0.06 1.6 90 10 -54 [4.8 • 10 -28 , 1.7 • 10 -26 ] 1.1 • 10 -26 7.7 • 10 -27 223 3 0.10 2.7 Bounds, parameters appearing in Algorithm 6.5 and running time of Algorithms 4.2 and 6.5 for the examples of Section 4.3: (i) y(x) = e x/2 / √ x + 16, (ii) y(x) = 3 2 cos x + 1 2 sin x, and (iii) y(x) = (cos x)/(2x 2 + 1).

The Clenshaw method we are referring to in this text should not be confused with the Hornerlike scheme for Chebyshev polynomials known as Clenshaw's algorithm[START_REF] Clenshaw | A note on the summation of Chebyshev series[END_REF].

References 37

Acknowledgements

We thank Alin Bostan, Nicolas Brisebarre, Élie de Panafieu, Bruno Salvy and Anne Vaugon for useful discussions and/or comments on various drafts of this work, and Moulay Barkatou for pointing out Ramis' method to us.

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,


Roger Faye

Salam Sawadogo

Felix Mora-Camino

A fuzzy modeling approach for the long term management of water resource systems

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


A new neural network structure dedicated to the control of nonlinear systems

Tsurng-Jehng Shen, Felix Mora-Camino, Karim Achaibou

To cite this version:

Tsurng-Jehng Shen, Felix Mora-Camino, Karim Achaibou. A new neural network structure dedicated to the control of nonlinear systems. MWSCAS 1995, 38th Midwest Symposium on Circuits and Systems, Aug 1995, Rio de Janeiro, Brazil. pp 433-436, 10.1109433-436, 10. /MWSCAS.1995.504469. hal-01022426 .504469. hal-01022426


Felix Mora-Camino

From flight qualities to flight control law design : a simulation tool

come