# Rigorous uniform approximation of D-finite functions using Chebyshev expansions

Alexandre Benoit, Mioara Joldes, Marc Mezzarobba

► To cite this version:

Alexandre Benoit, Mioara Joldes, Marc Mezzarobba. Rigorous uniform approximation of D-finite functions using Chebyshev expansions. Mathematics of Computation, 2017, 86 (305), pp.1303-1341. hal-01022420

HAL Id: hal-01022420

https://hal.science/hal-01022420

Submitted on 10 Jul 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 - Public Domain Dedication 4.0 International License
# RIGOROUS UNIFORM APPROXIMATION OF D-FINITE FUNCTIONS USING CHEBYSHEV EXPANSIONS

ALEXANDRE BENOIT, MIOARA JOLDES, AND MARC MEZZAROBBA

ABSTRACT. A wide range of numerical methods exists for computing polynomial approximations of solutions of ordinary differential equations based on Chebyshev series expansions or Chebyshev interpolation polynomials. We consider the application of such methods in the context of rigorous computing (where we need guarantees on the accuracy of the result), and from the complexity point of view.

It is well-known that the order-$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.

## CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>1.1. Background</td><td>2</td></tr><tr><td>1.2. Setting</td><td>3</td></tr><tr><td>1.3. Summary of Results</td><td>4</td></tr><tr><td>1.4. Outline</td><td>5</td></tr><tr><td>2. Chebyshev Expansions of D-finite Functions</td><td>5</td></tr><tr><td>2.1. Chebyshev Series</td><td>5</td></tr><tr><td>2.2. The Chebyshev Recurrence Relation</td><td>7</td></tr><tr><td>2.3. Solutions of the Chebyshev Recurrence</td><td>9</td></tr><tr><td>3. Convergent and Divergent Solutions</td><td>11</td></tr><tr><td>3.1. Elements of Birkhoff-Trjitzinsky Theory</td><td>11</td></tr></table>

<footnote n="2000">2000 Mathematics Subject Classification. 68W30, 33C45, 65L05, 65G20,</footnote>

<footnote>Key words and phrases. Rigorous computing, computer algebra, complexity, D-finite functions, recurrence relation, Chebyshev series, Clenshaw method, Miller algorithm, asymptotics, functional enclosure.</footnote>

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

<footnote>This work is in the public domain. As such, it is not subject to copyright. Where it is not legally possible to consider this work as released into the public domain, any entity is granted the right to use this work for any purpose, without any conditions, unless such conditions are required by law.</footnote>
<table><tr><td>3.2. Newton Polygon of a Chebyshev Recurrence</td><td>13</td></tr><tr><td>4. Computing Approximation Polynomials</td><td>15</td></tr><tr><td>4.1. Clenshaw's Algorithm Revisited</td><td>15</td></tr><tr><td>4.2. Convergence</td><td>18</td></tr><tr><td>4.3. Examples</td><td>22</td></tr><tr><td>4.4. A Link with the Tau Method</td><td>24</td></tr><tr><td>5. Chebyshev Expansions of Rational Functions</td><td>25</td></tr><tr><td>5.1. Recurrence and Explicit Expression</td><td>26</td></tr><tr><td>5.2. Truncation Error</td><td>27</td></tr><tr><td>5.3. Computation</td><td>28</td></tr><tr><td>6. Validation</td><td>30</td></tr><tr><td>6.1. Principle</td><td>30</td></tr><tr><td>6.2. Algorithm</td><td>32</td></tr><tr><td>Acknowledgements</td><td>37</td></tr><tr><td>References</td><td>37</td></tr></table>

# 1. INTRODUCTION

## 1.1. Background.
Many of the special functions commonly used in areas such as mathematical physics are so-called D-finite functions, that is, solutions of linear ordinary differential equations (LODE) with polynomial coefficients [56]. This property allows for a uniform theoretic and algorithmic treatment of these functions, an idea that was recognized long ago in Numerical Analysis [34, p. 464], and more recently found many applications in the context of Symbolic Computation [67, 54, 32]. The present article is devoted to the following problem.

**Problem 1.1.** Let $y : [-1, 1] \to \mathbb{R}$ be a D-finite function specified by a linear differential equation with polynomial coefficients and initial conditions. Let $d \in \mathbb{N}$. Given $y$ and $d$, find the coefficients of a polynomial $p(x) = \sum_{n=0}^{d} c_n T_n(x)$ written on the Chebyshev basis ($T_n$), together with a “small” bound $B$ such that $|y(x) - p(x)| \le B$ for all $x \in [-1, 1]$.

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 [58, 17] is a popular numerical computation system based on this idea.

In a yet more general setting, Epstein, Miranker and Rivlin developed a so-called “ultra-arithmetic” for functions that parallels floating-point arithmetic for real numbers [20, 21, 30]. Various generalized Fourier series, including Chebyshev series, play the role of floating-point numbers. Ultra-arithmetic also comprises a function space counterpart of interval arithmetic, based on truncated series with interval coefficients and rigorous remainder bounds. This line of approach was revived with the introduction of “ChebModels” in recent work by Brisebarre and Joldes [10]. Part of the motivation for Problem 1.1 is to allow one to use arbitrary D-finite functions as “base functions” at the leaves of expression trees to be evaluated using ChebModels.

Finally, perhaps the main appeal of ultra-arithmetic and related techniques is the ability to solve functional equations rigorously using enclosure methods [44, 30, 37, 45, 60]. LODE with polynomial coefficients are among the simplest equations to which these tools apply. A third goal of this article is to begin the study of the complexity of validated enclosure methods, from a computer algebra point of view, using this simple family of problems as a prototype.

## 1.2. Setting.
To specify the D-finite function $y$, we fix a linear homogeneous differential equation of order $r$ with polynomial coefficients

$$ (1.1) \qquad L \cdot y = a_{r}y^{(r)} + a_{r-1}y^{(r-1)} + \cdots + a_{0}y = 0, \quad a_{i} \in \mathbb{Q}[x]. $$

We also write $L = a_r \partial^r + \cdots + a_1 \partial + 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 $\| \cdot \|_{\infty}$. We also assume that $a_r(x) \neq 0$ for $x \in [-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) \qquad y^{(i)}(0) = \ell_i, \quad 0 \le i \le 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 $\llbracket i, j \rrbracket = \{i, i+1, \dots, 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
<table><caption>TABLE 1. Notation.</caption><tr><td>a<sub>i</sub></td><td>coefficients of the operator L; a<sub>r</sub>(0) ≠ 0</td><td>p. 3</td></tr><tr><td>b<sub>i</sub></td><td>coefficients of the operator P</td><td>p. 7</td></tr><tr><td>C</td><td>sequences in C<sup>Z</sup> with exponential decrease</td><td>p. 6</td></tr><tr><td>∂</td><td>differentiation operator, ∂ = d/dx</td><td>p. 3</td></tr><tr><td>L</td><td>differential operator, L · y = 0</td><td>p. 3</td></tr><tr><td>l<sub>i</sub></td><td>initial values, y<sup>(i)</sup>(0) = l<sub>i</sub></td><td>p. 3</td></tr><tr><td>P</td><td>Chebyshev recurrence operator</td><td>p. 7, p. 9</td></tr><tr><td>p<sub>d</sub><sup>*</sup></td><td>degree-d minimax polynomial approximation of y</td><td>p. 6</td></tr><tr><td>π<sub>d</sub></td><td>truncated Chebyshev expansion operator</td><td>p. 6</td></tr><tr><td>s</td><td>order of L</td><td>p. 3</td></tr><tr><td>S</td><td>shift operator, S : (u<sub>n</sub>) ↦→ (u<sub>n+1</sub>)</td><td>p. 7</td></tr><tr><td>s</td><td>(usually) half-order of P</td><td>p. 7, p. 12</td></tr><tr><td>S</td><td>singularities of P, shifted by s</td><td>p. 16</td></tr><tr><td>T<sub>n</sub></td><td>Chebyshev polynomials of the first kind</td><td>p. 5</td></tr><tr><td>y</td><td>(usually) unknown function, L · y = 0</td><td>p. 2, p. 26</td></tr><tr><td>y<sup>(N)</sup></td><td>approximation of y computed by Algorithm 4.2</td><td>p. 18</td></tr><tr><td>f̂</td><td>inverse Joukowski transform of a function f</td><td>p. 6</td></tr><tr><td>[[i, j]]</td><td>integer interval, [[i, j]] = {i, i + 1, ..., j}</td><td>p. 3</td></tr></table>

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 [65]. Our experiments show a satisfactory numerical behaviour.

To account for this variability in the underlying arithmetic, we assess the complexity of the algorithms in the arithmetic model. In other words, we only count field operations in Q, 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 [9] 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 [3, 29, 42].

## 1.4. Outline.
This article is organized as follows. In Section 2, we review properties of Chebyshev series of D-finite functions and then study the recurrence relations satisfied by the coefficients of these series, whose use is key to the linear time complexity. Section 3 provides results on the asymptotics of solutions of these recurrences that will be critical for the computation of the coefficients. The actual algorithm for this task, described in Section 4, reminds of Fox and Parker's variant [23, Chap. 5] of Clenshaw's algorithm [16]. A short description of a prototype implementation and several examples follow.

The part dedicated to the validation step starts in Section 5 with a study of Chebyshev series expansions of rational functions. Most importantly, we state remainder bounds that are then used in Section 6, along with an enclosure method for differential equations, to validate the output of the first stage and obtain the bound $B$. We conclude with examples of error bounds obtained using our implementation of the validation algorithm, and some open questions.

## 2. CHEBYSHEV EXPANSIONS OF D-FINITE FUNCTIONS

### 2.1. Chebyshev Series.
Recall that the Chebyshev polynomials of the first kind are polynomials $T_n(x) \in \mathbb{Q}[x]$ defined for all $n \in \mathbb{Z}$ by the relation $T_n(\cos \theta) = \cos(n\theta)$. They satisfy $T_{-n} = T_n$ for all $n$. The family $(T_n)_{n \in \mathbb{N}}$ is a sequence of orthogonal polynomials over $[-1, 1]$ with respect to the weight function $w(x) = 1/\sqrt{1-x^2}$, and hence a Hilbert basis of the space $L^2(w)$. (We refer the reader to books such as Rivlin's [53] or Mason and Handscomb's [39] for proofs of the results collected in this section.)

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

$$ (2.1) \qquad \sum_{n} u_n T_n = \frac{u_0}{2} T_0 + u_1 T_1 + u_2 T_2 + \dots, $$
we write Chebyshev series as

$$ (2.2) \qquad \sum_{n=-\infty}^{\infty} c_n T_n(x), \qquad 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 = \frac{1}{2}u_n$ of the expansion of a function $f$ are given by

$$ (2.3) \qquad c_n = \frac{1}{\pi} \int_{-1}^{1} \frac{f(x)T_n(x)}{\sqrt{1-x^2}} dx. $$

for all $n \in \mathbb{Z}$. The series (2.2) converges to $f$ in the $L^2(w)$ sense for all $f \in L^2(w)$ (see for example [39, Chap. 5.3.1]). We denote by $\pi_d : f \mapsto \sum_{n=-d}^{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 \subset \mathbb{C}$ that does not contain any singular point of the equation. Let

$$ (2.4) \qquad E_r = \{x \in \mathbb{C} : |x + \sqrt{x^2 - 1}| < r\} $$

be the largest elliptic domain with foci in $\pm 1$ with this property. Since the singular points are in finite number, we have $1 < r \le \infty$. The coefficients $c_n$ then satisfy $c_n = O(\alpha^n)$ for all $\alpha > r^{-1}$; and the Chebyshev expansion (2.2) of $f$ converges uniformly to $f$ on $E_r$ [39, Theorem 5.16]. Letting $x = \cos\theta$ and $z = e^{i\theta}$, it is not hard to see that the $c_n$ are also the coefficients of the (doubly infinite) Laurent expansion of the function $\hat{f}(z) = f(\frac{z+z^{-1}}{2})$ around the unit circle. The transformation $x = \frac{z+z^{-1}}{2}$ sending $f(x)$ to $\hat{f}(z)$ is known as the inverse Joukowski transform. It maps the elliptic disk $E_r$ to the annulus

$$ A_r = \{z \in \mathbb{C} : r^{-1} < |z| < r\}. $$

The formula $T_n(\cos\theta) = \cos(n\theta)$ translates into $T_n(\frac{z+z^{-1}}{2}) = \frac{z^n+z^{-n}}{2}$. The coefficients $c_n$ are also related to those of the Fourier cosine expansion of $\theta \mapsto f(\cos\theta)$.

Let $\mathcal{C} \subset \mathbb{C}^Z$ be the vector space of doubly infinite sequences $(c_n)_{n \in \mathbb{Z}}$ such that

$$ (\forall n \in \mathbb{N})(c_n = c_{-n}) \quad \text{and} \quad (\exists \alpha < 1)(c_n = O_{n \to \infty}(\alpha^n)). $$

The sequence of Chebyshev coefficients of a function $f$ that is analytic on some complex neighborhood of $[-1, 1]$ belongs to $\mathcal{C}$. Conversely, for all $c \in \mathcal{C}$, the function series $\sum_{n=-\infty}^{\infty} 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 [59, Theorem 16.1]

$$ (2.5) \qquad \|f - \pi_d(f)\|_{\infty} \le \left(\frac{4}{\pi^2} \ln(d+1) + 4\right) \|f - p_d^*\|_{\infty} $$

where $p_d^*$ is the polynomial of degree at most $d$ that minimizes $\|f - p\|_{\infty}$.

Even though $p_d^*$ itself can be computed to arbitrary precision using the Remez algorithm [12, Chap. 3], Equation (2.5) shows that we do not lose much by replacing it by $\pi_d(f)$. Moreover, tighter approximations are typically hard to validate without resorting to intermediate approximations of higher degree [13]. The need for such intermediate approximations is actually part of the motivation that led to
the present work. There exist a variety of other near-minimax approximations with nice analytical properties, e.g., Chebyshev interpolation polynomials. Our choice of truncated Chebyshev expansions is based primarily on the existence of a recurrence relation on the coefficients ($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) \qquad 2xT_n(x) = T_{n-1}(x) + T_{n+1}(x), $$

as well as the mixed differential-difference relation

$$ (2.7) \qquad 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 $\sum c'_n T_n = (\sum 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 [22, 23] in special cases and later proved in general by Paszkowski [49]. Properties of this recurrence and generalizations to other orthogonal polynomial bases were explored in a series of papers by Lewanowicz starting 1976 (see in particular [35, 36]). The automatic determination of this recurrence in a symbolic computation system was first studied by Geddes [24].

The following theorem summarizes results regarding this recurrence, extracted from existing work [49, 35, 36, 52, 4] and extended to fit our purposes. Here and in the sequel, we denote by $Q(n)\langle S, S^{-1} \rangle$ the skew Laurent polynomial ring over $Q(n)$ in the indeterminate $S$, subject to the commutation rules

$$ (2.8) \qquad S\lambda = \lambda S \quad (\lambda \in \mathbb{Q}), \quad Sn = (n+1)S. $$

Likewise, $Q[n]\langle S, S^{-1} \rangle \subset Q(n)\langle S, S^{-1} \rangle$ is the subring of noncommutative Laurent polynomials in $S$ themselves with polynomial coefficients. The elements of $Q[n]\langle S, S^{-1} \rangle$ identify naturally with linear recurrence operators through the left action of $Q[n]\langle S, S^{-1} \rangle$ on $C^Z$ defined by $(n \cdot u)_n = nu_n$ and $(S \cdot u)_n = u_{n+1}$. Recall that $L$ denotes the differential operator appearing in Equation (1.1).

**Theorem 2.1.** [49, 35, 36, 52, 4] Let $u, v$ be analytic functions on some complex neighborhood of the segment $[-1, 1]$, with Chebyshev expansions

$$ u(x) = \sum_{n=-\infty}^{\infty} u_n T_n(x), \quad v(x) = \sum_{n=-\infty}^{\infty} v_n T_n(x). $$

There exist difference operators $P, Q \in Q[n]\langle S, S^{-1} \rangle$ with the following properties.

(i) The differential equation $L \cdot u(x) = v(x)$ is satisfied if and only if

$$ (2.9) \qquad P \cdot (u_n) = Q \cdot (v_n). $$

(ii) The left-hand side operator $P$ is of the form $P = \sum_{k=-s}^{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) \qquad \delta_r(n) = 2^r \prod_{i=-r+1}^{r-1} (n-i), \quad I = \frac{1}{2n}(S^{-1} - S), $$
we have $Q = Q_r = \delta_r(n)I^r$ (this expression is to be interpreted as a polynomial identity in $Q(n)\langle S, S^{-1} \rangle$). 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 \in \mathbb{Z}}$ to the sequences $(u_n)_{n \in \mathbb{Z}\setminus\{0\}}$ defined only for nonzero $n$. A sloppy but perhaps more intuitive statement of the main point of Theorem 2.1 would be: "$(\int)^r L \cdot u = w$ if and only if $\delta_r(n)P \cdot u = w$, up to some integration constants".

*Proof.* Assume $L \cdot u = v$. Benoit and Salvy [4, Theorem 1] give a simple proof that (2.9) holds for some $P, Q \in Q(n)\langle S, S^{-1} \rangle$. 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 [49, 35]. More precisely, multiplying both members of [4, Eq. (17)] by $\delta_r(n)$ yields a recurrence of the prescribed form. The recurrence has polynomial coefficients since $\delta_r(n)I^r \in Q\langle S, S^{-1} \rangle$. Rebillard's thesis [52, Section 4.1] contains detailed proofs of this last observation and of all assertions of Item ii. Note that, although Rebillard's and Benoit and Salvy's works are closest to the formalism we use, several of these results actually go back to [49, 35, 36].

There remains to prove the “if” direction. Consider sequences $u, v \in C$ such that $P \cdot u = Q \cdot v$, and let $y \in C$ be the Chebyshev coefficient sequence of the (analytic) function $L \cdot u$. We then have $P \cdot u = Q \cdot y$ by the previous argument. This implies $Q \cdot v = Q \cdot y$, whence finally $y = v$ by Lemma 2.2 below. $\square$

**Lemma 2.2.** *The restriction to $C$ of the operator $Q$ from Theorem 2.1 is injective.*

*Proof.* With the notation of Theorem 2.1, we show by induction on $r \ge 1$ that
$$ (2.11) \qquad (v \in C) \wedge (|n| \ge r \implies (Q_r \cdot v)_n = 0) \implies v = 0. $$

First, we have $(\ker Q_1) \cap C = \{0\}$ since any sequence belonging to $C$ converges to zero as $n \to \pm\infty$. Now assume that (2.11) holds, and let $v \in C$ be such that $(Q_{r+1} \cdot v)_n = 0$ for $|n| \ge r+1$. Let $w = Q_r \cdot v$. Observe that $C$ is stable under the action of $Q(n)\langle S, S^{-1} \rangle$, so $w \in C$. Since $r \ge 1$, we have

$$ \begin{aligned} 2n Q_{r+1} &= \delta_{r+1}(n)(S^{-1} - S)I^r \\ &= 2((n+r)(n+r-1)S^{-1}\delta_r(n) - (n-r)(n-r+1)S\delta_r(n))I^r \\ &= 2((n+r)(n+r-1)S^{-1} - (n-r)(n-r+1)S)Q_r. \end{aligned} $$

Hence, for $|n| \ge r+1$, it holds that

$$ (2.12) \qquad (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} \to 1$ as $n \to \infty$, which contradicts the fact that $w \in C$. It follows that $w_n = 0$ for $|n|$ large enough, and, using (2.12) again, that $w_n = 0$ as soon as $|n| \ge r$. Applying the hypothesis (2.11) concludes the induction. $\square$

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 = \frac{1}{2}(z + z^{-1})$ in the differential equation (1.1). Then, we compute a recurrence on the Laurent coefficients of $\hat{u}(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) + 2xy'(x) = 0$. The substitutions

$$x = \frac{z + z^{-1}}{2} \qquad \frac{d}{dx} = \frac{2z}{z - z^{-1}} \frac{d}{dz}$$

yield (after clearing common factors and denominators)

$$(z+1)(z-1)(z^4 + 18z^2 + 1) \hat{y}''(z) + 2(z^4 - 2z^2 - 19)z \hat{y}'(z) = 0.$$

We then set $\hat{y}(z) = \sum_{n=-\infty}^{\infty} c_n z^n$ and extract the coefficient of $z^n$ (which amounts to replacing $z$ by $S^{-1}$ and $z \frac{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 [4] give a unified presentation of several alternative algorithms, including Paszkowski's, by interpreting them as various ways to perform the substitution $x \mapsto \frac{1}{2}(S + S^{-1})$, $\frac{d}{dx} \mapsto (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 = \frac{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 \in Q(n)\langle S, S^{-1} \rangle$.

**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 \in \mathbb{Z}}$, we have the equalities

$$\forall n, \sum_k b_k(n)u_{n+k} = -\sum_k b_{-k}(-n)u_{n+k} = -\sum_k b_k(-n)u_{-n-k},$$

that is, $P \cdot (u_n)_{n \in \mathbb{Z}} = -P \cdot (u_{-n})_{n \in \mathbb{Z}}$. In particular, if $(u_n)_{n \in \mathbb{Z}}$ is a solution of a homogeneous Chebyshev recurrence, then so is $(u_{-n})_{n \in \mathbb{Z}}$, 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$.

**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, \dots$ 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” <fnref n="1" />.

Next, also in contrast to the case of power series, the leading and trailing coefficients $b_{\pm 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 = \sum y_{|n|}T_n$ of degree $d$ is a solution of $L \cdot 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 \in \mathbb{Z}$, the Chebyshev recurrence relation associated to the differential equation $y''(x) + (x^2 + 1)y'(x) - kxy(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) \qquad b_{j-i}(-j) = -b_{j+i}(-j), \quad |j| \le r-1, \quad i \in \mathbb{N}, $$

with $b_k = 0$ for $|k| > s$. In particular, $b_s(n)$ is zero for all $n \in [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 \ge 2$ and assume that the proposition holds when $L$ has order $r-1$. Write $L = L^\flat + \partial^r p_r(x)$ where $p_r \in \mathbb{Q}[x]$ and $L^\flat$ is a differential operator of order at most $r-1$. Letting $P^\flat = \sum_{k \in \mathbb{Z}} b_k^\flat(n)S^k$ be the Chebyshev recurrence operator associated to $L^\flat$, we then have [4]

$$ (2.14) \qquad \delta_r(n)^{-1}P = I\delta_{r-1}(n)^{-1}P^\flat + p_r(\frac{1}{2}(S + S^{-1})) $$

where the last term denotes the evaluation of $p_r$ at $x = \frac{1}{2}(S + S^{-1})$. Since

$$ I\delta_{r-1}(n)^{-1} = (n\delta_r(n))^{-1}((n-r+2)(n-r+1)S^{-1} - (n+r-2)(n+r-1)S) $$

<footnote n="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 [3] for some results.</footnote>
by the commutation rule (2.8), relation (2.14) rewrites as

$$P = \frac{1}{n} \sum_k ((n-r+2)(n-r+1)b_{k+1}^b(n-1) - (n+r-2)(n+r-1)b_{k-1}^b(n+1))S^k + \delta_r(n)p_r(\frac{1}{2}(S+S^{-1})).$$

The case $j = 0$ having already been dealt with, assume $0 < |j| < r$. Since $\delta_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

$$(2.15) \quad -jb_k(-j) = (j+r-2)(j+r-1)b_{k+1}^b(-j-1) - (j-r+2)(j-r+1)b_{k-1}^b(-j+1).$$

Now $b_{j-i}^b(-j) = -b_{j+i}^b(-j)$ for $|j| < r-1$ by the induction hypothesis, and the term involving $b_{k\pm 1}^b$ vanishes for $j = \mp(r-1)$ and $j = \mp(r-2)$. In each case, we obtain $b_{j-i}(-j) = -b_{j+i}(-j)$. $\square$

**Corollary 2.9.** Let $P$ be the Chebyshev recurrence operator associated to $L$. The image by $P$ of a symmetric sequence $(u|_n|)_{n \in \mathbb{Z}}$ satisfies $(P \cdot u)_n = 0$ for $|n| < r$.

*Proof.* Since

$$(P \cdot u)_n = \sum_{k \in \mathbb{Z}} b_k(n) u_{n+k} = \sum_{i \in \mathbb{Z}} b_{i-n}(n) u_i,$$

it follows from Proposition 2.8 with $j = -n$ and $|n| < r$ that

$$\sum_{i \in \mathbb{Z}} b_{i-n}(n) u_i = - \sum_{i \in \mathbb{Z}} b_{-i-n}(n) u_i = - \sum_{i \in \mathbb{Z}} b_{i-n}(n) u_i,$$

that is, $(P \cdot u)_n = -(P \cdot u)_n$. $\square$

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 \cdot u)_n = u(n+1) + 2n u(n) - u(n-1) = 0.$$

In terms of the modified Bessel functions $I_\nu$ and $K_\nu$, a basis of solutions of the recurrence is given by the sequences $(I_\nu(1))_{\nu \in \mathbb{Z}}$ and $(K_\nu(1))_{\nu \in \mathbb{Z}}$. The former is the coefficient sequence of the Chebyshev expansion of the exponential function and decreases as $\Theta(2^{-\nu} \nu!^{-1})$. The later satisfies $K_\nu(1) = \Theta(2^\nu (\nu - 1)!)$.

### 3. CONVERGENT AND DIVERGENT SOLUTIONS

#### 3.1. Elements of Birkhoff-Trjitzinsky Theory.
Before studying in more detail the convergent and divergent solutions of the Chebyshev recurrence relation, we recall some elements of the asymptotic theory of linear difference equations. Much of the presentation is based on Wimp's book [65, Appendix B], to which we refer the reader for more information.
**Definition 3.1.** For all $\rho \in \mathbb{N}\setminus\{0\}$, $J \in \mathbb{N}$, $\kappa \in \mathbb{Q}$, $\alpha, \pi_j, \theta, \beta_{j,i} \in \mathbb{C}$, we call the formal expansion

$$ (3.1) \qquad \bar{u}(n) = n!^{\kappa} \alpha^n e^{\pi(n)} \sum_{j=0}^{J} (\ln n)^j \sum_{i=0}^{\infty} \beta_{j,i} n^{\theta-i/\rho} $$

where

$$ \pi(n) = \pi_1 n^{1/\rho} + \dots + \pi_{\rho-1} n^{(\rho-1)/\rho} $$

a formal asymptotic series (FAS). The set of all FAS is denoted by $\mathcal{B}$.

Formal asymptotic series are to be interpreted as asymptotic expansions of sequences as $n \to \infty$. The product of two FAS is defined in the obvious way and is again an FAS. The same goes for the substitution $n \mapsto n + k$ for fixed $k \in \mathbb{Z}$, using identities such as $(n+k)^\theta = n^\theta(1+k\theta n^{-1}+\dots)$. The sum of two FAS is not always an FAS, but that of two FAS sharing the same parameters $\kappa, \alpha, \pi$ is. Thus, it makes sense to say that an FAS $\bar{u} \in \mathcal{B}$ satisfies a recurrence

$$ (3.2) \qquad \bar{b}_s(n)\bar{u}(n+s) + \dots + \bar{b}_0(n)\bar{u}(n) = 0 $$

with formal series coefficients of the form

$$ (3.3) \qquad \bar{b}_k(n) = n^{\tau_k/\omega} (\beta_{k,0} + \beta_{k,1} n^{-1/\omega} + \beta_{k,2} n^{-2/\omega} + \dots) \in C((n^{-1/\omega})). $$

Also, given s FAS $\bar{u}_0, \dots, \bar{u}_{s-1} \in \mathcal{B}$, the Casoratian

$$ C(n) = \det(\bar{u}_j(n+i))_{0 \le i,j < s} $$

belongs to $\mathcal{B}$ as well.

Following Wimp, we say that $\bar{u}_1, \dots, \bar{u}_s \in \mathcal{B}$ are formally linearly independent when their Casoratian is nonzero. Note that the elements of any subset of $\{\bar{u}_1, \dots, \bar{u}_s\}$ are then formally linearly independent as well. Indeed, it can be checked by induction on s that s FAS $\bar{u}_1, \dots, \bar{u}_s$ are formally linearly dependent if and only if there exists a relation of the form $\bar{\mu}_1(n)\bar{u}_1(n) + \dots + \bar{\mu}_n(s)\bar{u}_s(s) = 0$ where the $\bar{\mu}_k$ are FAS such that$^2 \bar{\mu}_k(n+1) = \bar{\mu}_k(n)$.

**Definition 3.2.** The FAS (3.1) is said to be an asymptotic expansion of a sequence $(u_n) \in \mathbb{C}^\mathbb{N}$, and we write $u_n \sim \bar{u}(n)$, when for any truncation order I, the relation

$$ u_n = n!^{\kappa} \alpha^n e^{\pi(n)} \sum_{j=0}^{J} (\ln n)^j \left( \sum_{i=0}^{I-1} \beta_{j,i} n^{\theta-i/\rho} + O(n^{\theta-I/\rho}) \right) $$

holds as $n \to \infty$.

The following fundamental result is known as the Birkhoff-Trjitzinsky theorem, or “main asymptotic existence theorem” for linear recurrences. It will be the starting point of our analysis of the computation of “convergent” solutions of the Chebyshev recurrence by backward recurrence.

**Theorem 3.3.** [6, 7, 61, 27] *Consider a linear recurrence*

$$ (3.4) \qquad b_s(n)u_{n+s} + \dots + b_0(n)u_n = 0 $$

whose coefficients $b_0, \dots, b_s$ admit asymptotic expansions (in the sense of Definition 3.2) $\bar{b}_0, \dots, \bar{b}_s$ of the form (3.3) for some integer $\omega \ge 1$. Then,

<footnote>$^2$Like Wimp, but unlike most authors, we consider recurrences rather than difference equations. Accordingly, we forbid factors of the form $e^{\pi \rho n}$ with $|\operatorname{Im} \pi \rho| > \pi$ in (3.1), so that the $\mu_k(n)$ are actually constants in our setting.</footnote>
(i) the (formal) recurrence (3.2) possesses a system of s formally linearly in-
dependent FAS solutions;

(ii) for any s formally linearly independent solutions $\bar{e}_1, \dots, \bar{e}_s \in \mathcal{B}$ of (3.2), there exists complex sequences $e_1 = (e_{1,n})_{n \ge N}, \dots, e_s = (e_{s,n})_{n \ge N}$ defined in some neighborhood of infinity, with the property that $e_k \sim \bar{e}_k$ for all $k$, and such that $(e_1, \dots, e_s)$ is a basis of the solution space of (3.4) for $n \ge 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 [27, 62].
Besides, under mild additional assumptions on the Chebyshev recurrence, all the
information needed in our analysis is already provided by the more elementary
Perron-Kreuser theorem<fnref n="3" /> (cf. [26, 41, 43]) or its extensions by Schäfke [55]. See
also Immink [28] and the references therein for an alternative approach in the
case of recurrences with polynomial coefficients, originating in unpublished work
by Ramis.

Also observe that for any subfamily $(f_1, \dots, f_{s'})$ of the sequences $e_i$ from The-
orem 3.3, the matrix $(f_{j,n+i})_{1 \le i,j \le 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, \dots, (e_{s-1,n})_n$ admit formally lin-
early independent asymptotic expansions of the form (3.1), with $\alpha_i \in C \setminus \{0\}$,
$\kappa_i \in Q$. Then the Casorati determinant

$$
C(n) = \begin{vmatrix}
e_{0,n} & e_{1,n} & \cdots & e_{s-1,n} \\
e_{0,n+1} & e_{1,n+1} & & e_{s-1,n+1} \\
\vdots & \vdots & & \vdots \\
e_{0,n+s-1} & e_{1,n+s-1} & \cdots & e_{s-1,n+s-1}
\end{vmatrix}
$$

satisfies

$$
C(n) = \beta e_{0,n} e_{1,n} \cdots e_{s-1,n} n^{\theta} ((\ln n)^{\lambda} + O((\ln n)^{\lambda-1})), \quad n \to \infty,
$$

for some $\beta \in C\setminus\{0\}$, $\theta \in C$, and $\lambda \in N$.

Proof. Write $C(n) = e_{0,n}e_{1,n} \cdots e_{s-1,n}C'(n)$. The formal linear independence hy-
pothesis means that $C(n)$, and hence $C'(n)$, admit nonzero FAS as asymptotic
expansions. Additionally,

$$
C'(n) = \det \left( \frac{e_{j,n+i}}{e_{j,n}} \right)_{0 \le i,j < s}
$$

has at most polynomial growth, so that the leading term of its asymptotic expansion
must be of the form $n^\theta (\ln n)^\lambda$. $\square$

**3.2. Newton Polygon of a Chebyshev Recurrence.** The formal solutions described in Theorem 3.3 may be constructed algorithmically using methods going back to Poincaré [50] and developed by many authors. See in particular Adams [1] and Birkhoff [6] for early history, Tournier [57] for a comparison of several methods

<footnote n="3"><sup>3</sup>The Perron-Kreuser theorem yields the existence of a basis of solutions such that $e_{i,n+1}/e_{i,n} \sim \alpha n^{\kappa_i}$, under the assumption that $\kappa_i = \kappa_j \Rightarrow |\alpha_i| \neq |\alpha_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.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. The Newton polygon of a Chebyshev recurrence.</figcaption></figure>

from a Computer Algebra perspective, and Balser and Bothner [2] for a modern algorithm as well as more references.

Here, we are mostly interested in the parameters $\kappa$ and $\alpha$ 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} + \dots + b_0(n) + \dots + b_s(n)S^s$$

from Section 2.2. The Newton polygon of $P$ is defined as the lower convex hull of the points $p_k = (k, -\deg b_k) \in \mathbb{R}^2$ (see Figure 1). To each edge $[p_i, p_j]$ ($i < j$) of the polygon is attached a characteristic equation

$$\chi_i(\alpha) = \sum_{k: p_k \in [p_i, p_j]} \ln(b_k)\alpha^{k-i},$$

where $\ln(p)$ denotes the leading coefficient of $p$. Note that the degrees of the $\chi_i$ sum to 2s. Let

$\alpha_s, \alpha_{s-1}, \dots, \alpha_1, \alpha_{-1}, \dots, \alpha_{-s+1}, \alpha_{-s}$

be the sequence of all roots of the polynomials $\chi_i$, with multiplicities, the roots corresponding to distinct edges being written in the order of increasing $i$ and the roots of each $\chi_i$ in that of increasing modulus. For all $k$, let $\kappa_k$ be the slope of the edge associated to $\alpha_k$. (Thus, each $\kappa_k$ is repeated a number of times equal to the horizontal length of the corresponding edge, and we have $\kappa_s \le \kappa_{s-1} \le \dots \le \kappa_{-s}$.)

How does this relate to the asymptotics of Chebyshev series? Assume that $n!\kappa^n$ is the leading factor of some FAS solution $\bar{u}$ of $P \cdot \bar{u} = 0$. It is not too hard to see that, in order for asymptotically dominant terms of $P \cdot \bar{u}$ to cancel out, $\kappa$ must be among the slopes of the Newton polygon of $P$, and $\alpha$ must be a root of the characteristic equation of the corresponding edge. This gives all possible values of $\kappa$ and $\alpha$. Conversely, the construction behind Theorem 3.3 (i) yields a number of linearly independent FAS with given $\kappa$ and $\alpha$ equal to the multiplicity of $\alpha$ as a root of the characteristic equation of the edge of slope $\kappa$. In the case of Chebyshev recurrences, the Newton polygon has the following symmetry property.

**Proposition 3.5.** The slopes $\kappa_i$ of the Newton polygon of $P$ and the roots $\alpha_i$ of its characteristic equations satisfy $\kappa_{-i} = -\kappa_i$ and $|\alpha_{-i}| = |\alpha_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
$\kappa_{-i} = -\kappa_i$ for all $i$. Now fix $i$, and let $\epsilon_i = [p_{\ell(i)}, p_{r(i)}]$ be the edge of slope $\kappa_i$. The characteristic equation of $\epsilon_i$ reads

$$ \chi_i(\alpha) = \sum_{k:p_k \in \epsilon_i} \mathrm{lc}(b_k) \alpha^{k-\ell(i)} = \sum_{k:p_k \in \epsilon_i} (-1)^{1+\deg b_k} \mathrm{lc}(b_{-k}) \alpha^{k-\ell(i)}, $$

where $\mathrm{lc}(b)$ denotes the leading coefficient of a polynomial $b$. Using the relation $\deg b_k - \deg b_{\ell(i)} = \kappa_i(k - \ell(i))$ for $p_k$ lying on $\epsilon_i$, we get

$$ \begin{aligned} \chi_i(\alpha) &= \pm \sum_{k:p_k \in \epsilon_{-i}} (-1)^{\kappa_i(-k+\ell(-i))} \mathrm{lc}(b_k) \alpha^{-k+\ell(-i)} \\ &= \pm \alpha^{\ell(i)-\ell(-i)} \chi_{-i}((-1)^{-\kappa_i} \alpha^{-1}), \end{aligned} $$

and hence $|\alpha_{-i}| = |\alpha_i|^{-1}$.

There remains to prove that $\kappa_i = 0$ implies $|\alpha_i| \neq 1$. Under the change of variable $x = \frac{1}{2}(z + z^{-1})$, the leading term with respect to $\theta = z \frac{d}{dz}$ of $(\frac{d}{dx})^k$ is $2^k(z - z^{-1})^{-k}$. (The leading term is well-defined because the commutation relation between $z$ and $\theta$ preserves degrees.) Therefore, the characteristic equation associated to the slope $\kappa = 0$ (when there is one) of the recurrence operator $P_1$ obtained by changing $z$ into $S^{-1}$ and $\theta$ into $n$ is

$$ \chi_{\mathrm{horiz}}(\alpha) := \sum_{k:\deg p_k = \max_i \deg p_i} \mathrm{lc}(b_k) \alpha^{k-i} = a_r \left( \frac{\alpha + \alpha^{-1}}{2} \right), $$

where $a_r$ is the leading coefficient of (1.1). Since $P$ is a right factor of $P_1$, the characteristic polynomial associated to $\kappa = 0$ in the Newton polygon of $P$ divides $\chi_{\mathrm{horiz}}$. But, due to the assumptions stated in Section 1.2, the polynomial $a_r(x)$ does not vanish for $x \in [-1, 1]$, hence $\chi_{\mathrm{horiz}}(\alpha) \neq 0$ for $|\alpha| = 1$. $\square$

Summing up, the asymptotic structure of the solutions of the Chebyshev recurrence may be described as follows. Similar observations were already made by Rebillard [52, Chap. 5].

**Corollary 3.6.** For large enough $N$, the space of sequences $(u_n)_{n \ge N}$ satisfying $(P \cdot u)_n = 0$ (“germs of solution at infinity of the Chebyshev recurrence”) has a basis comprising $s$ convergent sequences $e_1, \dots, e_s$ and $s$ divergent sequences $e_{-1}, \dots, e_{-s}$, all with formally linearly independent FAS expansions, such that

$$ e_{i,n} = n!^{\kappa_i} \alpha_i^n e^{o(n)}, \quad n \to \infty. $$

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. $\square$

# 4. 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 [16]. The link between Clenshaw's method<fnref n="4" /> and

<footnote n="4">⁴The Clenshaw method we are referring to in this text should not be confused with the Horner-like scheme for Chebyshev polynomials known as Clenshaw's algorithm [15].</footnote>
the Chebyshev recurrence was observed long ago by Fox and Parker [22, 23] and further discussed by Rebillard [52, 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 well-known backward recurrence technique [5, 65] to compute minimal (“slowest increasing”) solutions of three-term recurrences. Miller’s idea is to compute the coefficients $u_N, u_{N-1}, \dots, u_0$ of a linear recurrence sequence in the backward direction, starting from 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, \dots, 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 than 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 $\mathbf{S} = \{n \ge s : b_{-s}(n) = 0\}$. Let

$$ \mathcal{E} = \{(u_{|n|})_{n \in \mathbb{Z}} : n \in \mathbb{N} \setminus \mathbf{S} \Rightarrow (P \cdot u)_n = 0\} $$

be the space of symmetric sequences whose restriction to $n \in \mathbb{N}$ satisfies the Chebyshev recurrence, except possibly when $n \in \mathbf{S}$.

**Proposition 4.1.** (i) The space $\mathcal{S}$ of symmetric sequences $(u_{|n|})_{n \in \mathbb{Z}}$ such that $P \cdot u = 0$ has dimension $s + r$. Among the elements of $\mathcal{E}$, these sequences are characterized by the linear equations

$$ (4.1) \qquad (P \cdot u)_n = 0, \quad n \in [r, s-1] \cup \mathbf{S}. $$

These equations are linearly independent.

(ii) The space $\mathcal{S} \cap \mathcal{C}$ of symmetric, convergent sequences $(u_{|n|})_{n \in \mathbb{Z}}$ such that $P \cdot u = 0$ has dimension $r$. Its elements are the elements of $\mathcal{E}$ satisfying (4.1), and the equations (4.1) are independent as linear forms on $\mathcal{E} \cap \mathcal{C}$ as well.

*Proof.* First observe that a sequence $u \in \mathcal{E}$ automatically satisfies $(P \cdot u)_n = 0$ for $-n \in \mathbb{N} \setminus \mathbf{S}$ too, by Remark 2.6. Then $u$ belongs to $\mathcal{S}$ if and only if $(P \cdot u)_n = 0$ for $n \in \mathbf{S}$ and for $|n| < s$. But the equations $(P \cdot u)_n = 0$ with $|n| < r$ are trivial by Corollary 2.9, thus $u \in \mathcal{S}$ if and only if $u \in \mathcal{E}$ and (4.1) is satisfied.

Let $t = |\mathbf{S}|$, and let $p \le s - r + t$ denote the rank of (4.1), considered as a system of linear forms over $\mathcal{E}$. For large enough $n_0$, sequences $u \in \mathcal{E}$ are in bijection with their values $(u_n)_{n \in \mathbf{J}}$, where $\mathbf{J} = (\mathbf{S} - s) \cup [n_0, n_0 + 2s - 1]$, hence $\mathcal{E}$ has dimension $2s + t$. It follows that $\dim \mathcal{S} = \dim \mathcal{E} - p \ge s + r$.

Similarly, an element of $\mathcal{E} \cap \mathcal{C}$ is characterized, for large $n_0$, by a convergent sequence $(u_n)_{n \ge n_0}$ (“a convergent germ of solution”) and values $u_n$ for $n \in \mathbf{S}$. It belongs to $\mathcal{S}$ when additionally (4.1) is satisfied. Thus, by Corollary 3.6, we have
$dim(E \cap C) = s + t$ and $dim(S \cap C) = s + t - q$, where $q \le p$ is the rank of the system (4.1) restricted to $E \cap C$. But we already know from Theorem 2.1 that $S \cap C \simeq \ker L$, and hence $dim(S \cap C = r$. Therefore, we have $r = s+t-q \ge s+t-p \ge r$, and hence $p = q = s+t-r$. $\square$

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) = \ell_i$ as a special case of boundary conditions

$$ (4.2) \qquad \lambda_i(y) = \ell_i, \quad 1 \le i \le r, $$

each of the form $\lambda_i(y) = \sum_{j=1}^q \mu_j y^{(r_j)}(x_j)$ with $x_j \in [-1, 1]$ and $r_j \le 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 $\lambda_i : \ker L \to C$ are linearly independent.

Motivated by Proposition 4.1, we “unroll” $s+t$ linearly independent test sequences $f_i \in 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.

**Algorithm 4.2.** Input: a linear differential operator $L$ of order $r$, boundary conditions $\lambda_1(y) = \ell_1, \dots, \lambda_r(y) = \ell_r$ as in (4.2), a target degree $d > s$, an integer $N \ge \max(d, \max\{n : b_{-s}(n) = 0\})$. Output: an approximation $\tilde{y}(x) = \sum_{n=-d}^d \tilde{y}_n T_n(x)$ of the corresponding solution $y$ of $L \cdot y = 0$.

1 compute the Chebyshev recurrence operator $P = \sum_{k=-s}^s b_k(n)S^k$ associated to $L$

2 set $S = \{n \ge s : b_{-s}(n) = 0\}$ and $I = S \cup [N, N+s-1]$

3 for $n$ from $N+2s-1$ down to 1

4 for $i \in I$

5 if $n = i$ then set $f_{i,n-s} = 1$

6 else if $n \in I$ or $n \ge N+s$ then set $f_{i,n-s} = 0$

7 else compute $f_{i,n-s}$ using the relation $(P \cdot f)_n = 0$

8 using indeterminates $\eta_i$, $i \in I$, set

$$ \tilde{y}_n = \begin{cases} \sum_{i \in I} \eta_i f_{i,|n|}, & |n| \le N \\ \tilde{y}_n = 0, & |n| > N, \end{cases} \quad \text{and} \quad \tilde{y}(x) = \sum_{n=-N}^{N} \tilde{y}_n T_n(x) $$

9 solve for $(\eta_i)_{i \in I}$ the linear system

$$ (4.3) \qquad \left\{ \begin{array}{ll} \lambda_k(\tilde{y}) = \ell_k, & 1 \le k \le r, \\ b_{-s}(n)\tilde{y}_{n-s} + \dots + b_s(n)\tilde{y}_{n+s} = 0, & n \in [r, s-1] \cup S \end{array} \right. $$

10 return $\sum_{n=-d}^d \tilde{y}_n T_n(x)$

The complexity of Algorithm 4.2 is easy to estimate.
**Proposition 4.3.** For fixed $L$, $\lambda_i$, and $\ell_i$, Algorithm 4.2 runs in $O(N)$ arithmetic operations.

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 $\mathcal{E}$ by Proposition 4.1) always has a nontrivial intersection with $\text{Span}\{t_i : i \in \mathbf{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\|_{\infty}$) 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 [66, 65], 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 [64].

Retaining the notation from the previous subsection, assume that the operator $L$ and the boundary conditions $\lambda_i(y) = \ell_i$ are fixed. Write the Chebyshev expansion of $y$ as

$$y(x) = \sum_{n=-\infty}^{\infty} y_n T_n(x).$$

Let $y_n^{(N)} = \tilde{y}_n$, $|n| \le 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 \to \infty$, the polynomial output by Algorithm 4.2 converges at least exponentially fast to the truncated Chebyshev series $\pi_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 \to \infty$, its output satisfies

$$\max_{n=-N}^{N} |y_n^{(N)} - y_n| = O(N^{\tau} e_{1,N})$$

for some $\tau$ independent of $N$.

We write $f(N) = O_{\text{pol}}(g(N))$ when there exists $\tau \ge 0$ such that $f(N) = O(N^{\tau}g(N))$ as $n \to \infty$.

*Proof.* The finite sequence $(y_n^{(N)})_{n=-N}^N$ computed by Algorithm 4.2 extends to an element $(y_n^{(N)})_{n \in \mathbb{Z}}$ of $\mathcal{E}$ characterized by the conditions

$$y_N^{(N)} = \dots = y_{N+s-1}^{(N)} = 0$$

from Step 3, along with the linear system (4.3) solved in Step 9.
By writing the linear forms $\lambda_1, \dots, \lambda_r : \mathcal{C} \to \mathcal{C}$ that express the boundary conditions (4.2) as $\lambda_i(y) = \sum_{n=-\infty}^{\infty} \lambda_{i,n}y_n$, we define “truncations”

$$ \lambda_i^{(N)}(y) = \sum_{n=-N}^{N} \lambda_{i,n} y_n $$

that make sense even for divergent series. (Abusing notation slightly, we apply the $\lambda_i$ and $\lambda_i^{(N)}$ indifferently to functions, formal Chebyshev series or their coefficient sequences.) The system (4.3) consists of the equations $\lambda_i^{(N)}(y^{(N)}) = \ell_i$ and of the symmetry and extension-through-singularities constraints (4.1). We introduce additional linear forms $\lambda_{r+1} = \lambda_{r+1}^{(N)}, \dots, \lambda_{s+t} = \lambda_{s+t}^{(N)}$ to write these last $s-r+t$ equations in the same form as the first $r$, so that (4.3) becomes

$$ (4.4) \qquad \lambda_i^{(N)}(y^{(N)}) = \sum_{n=-N}^{N} \lambda_{i,n} y_n^{(N)} = \ell_i, \quad 1 \le i \le s+t. $$

Now let $(e_1, \dots, e_s, e_{-1}, \dots, e_{-s})$ be a basis of the solutions of $P \cdot u = 0$ in the neighborhood of $+\infty$ of the form provided by Corollary 3.6. Extend each $e_i$ to an element of $\mathcal{E}$, and then the tuple to a basis of $\mathcal{E}$ by setting $e_{s+1} = f_{n_1}, \dots, e_{s+t} = f_{n_t}$ where $n_1 < n_2 < \dots < n_t$ are the elements of $\mathbf{S}$. Let

$$ (4.5) \qquad \Delta^{(N)} = \begin{vmatrix} e_{1,N} & \cdots & e_{s+t,N} & e_{-1,N} & \cdots & e_{-s,N} \\ \vdots & & \vdots & \vdots & & \vdots \\ e_{1,N+s-1} & \cdots & e_{s+t,N+s-1} & e_{-1,N+s-1} & \cdots & e_{-s,N+s-1} \\ \lambda_1^{(N)}(e_1) & \cdots & \lambda_1^{(N)}(e_{s+t}) & \lambda_1^{(N)}(e_{-1}) & \cdots & \lambda_1^{(N)}(e_{-s}) \\ \vdots & & \vdots & \vdots & & \vdots \\ \lambda_{s+t}^{(N)}(e_1) & \cdots & \lambda_{s+t}^{(N)}(e_{s+t}) & \lambda_{s+t}^{(N)}(e_{-1}) & \cdots & \lambda_{s+t}^{(N)}(e_{-s}) \end{vmatrix}, $$

Let $\Delta_j^{(N)}$ be the same determinant with the column involving $e_j$ replaced by

$$ \underbrace{(0, \dots, 0, \ell_1, \dots, \ell_r)}_{s \text{ times}} \underbrace{, 0, \dots, 0}_{s-r+t \text{ times}}^{\mathrm{T}}. $$

By Cramer's rule, provided $\Delta^{(N)} \neq 0$, the sequence $(y_n^{(N)})_n$ decomposes on the basis $(e_j)_{j=-s}^{s+t}$ of $\ker P \subset \mathbb{C}^L$ as

$$ (4.6) \qquad y^{(N)} = \sum_{k=-s}^{s+t} \gamma_k^{(N)} e_k, \qquad \gamma_k^{(N)} = \frac{\Delta_k^{(N)}}{\Delta^{(N)}}. $$

Algorithm 4.2 fails if and only if $\Delta^{(N)} = 0$.

The sequence of “exact” Chebyshev coefficients of the function $y$ defined by the input is likewise given by

$$ (4.7) \qquad y = \sum_{k=1}^{s+t} \gamma_k e_k, \qquad \gamma_k = \frac{\Delta_k}{\Delta}, $$

where

$$ \Delta = \det(\lambda_i(e_j))_{1 \le i,j \le s+t} $$

and $\Delta_j$ denotes the determinant $\Delta$ with the $j$-th column replaced by $(\ell_1, \dots, \ell_{s+t})^\mathrm{T}$.
Our goal is now to prove that $\gamma_k^{(N)} \to \gamma_k$ fast as $N \to \infty$. To do that, we study the asymptotic behaviours of the determinants $\Delta_k^{(N)}$ and $\Delta_k^{(N)}$.

We decompose $\Delta_k^{(N)}$ into the four blocks indicated by Eq. (4.5) as follows:

$$ \Delta_k^{(N)} = \begin{vmatrix} A & B \\ C & D \end{vmatrix}. $$

The corresponding modified blocks in $\Delta_k^{(N)}$ 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 \to \Delta \neq 0$ as $N \to \infty$. The Schur complement formula implies

$$ \Delta_k^{(N)} = -\det(B)\det(C)\det(I - C^{-1}DB^{-1}A). $$

Setting $\mathbf{e}_j = (e_{j,N}, \dots, 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

$$ \begin{aligned} (B^{-1}A)_{i,j} &= \frac{\det(\mathbf{e}_{-1}, \dots, \mathbf{e}_{-i+1}, \mathbf{e}_j, \mathbf{e}_{-i-1}, \dots, \mathbf{e}_{-s})}{\det B} \\ &= \frac{(-1)^{i-1} \det(\mathbf{e}_j, \mathbf{e}_{-1}, \dots, \widehat{\mathbf{e}_{-i}}, \dots, \mathbf{e}_{-s})}{\det B} \\ &= O_{\text{pol}}\left(\frac{e_{j,N}}{e_{-i,N}}\right) \end{aligned} $$

(where the notation $\widehat{\cdot}$ indicates the omission of the corresponding term) as $N \to \infty$ by Lemma 3.4.

In view of our assumptions on the boundary conditions (4.2), we have $\lambda_{i,n} = O_{n \to \pm\infty}(n^r)$ for all $i \le \sigma$, where $\sigma$ is the maximum derivation order appearing in (4.2). Additionally, the sequences $\lambda_{i,n}$ with $r+1 \le i \le s+t$ are ultimately zero. Therefore the entries of $D$ satisfy

$$ D_{i,j} = \lambda_i^{(N)}(e_{-j}) = O_{\text{pol}}(e_{-j,N}). $$

This yields the estimate

$$ (DB^{-1}A)_{i,j} = O_{\text{pol}}(e_{j,N}) $$

for the $j$-th column of $DB^{-1}A$. Since

$$ C_{i,j} = \lambda_i^{(N)}(e_j) = \lambda_i(e_j) + O_{\text{pol}}(e_{j,N}), $$

we get $(C^{-1}DB^{-1}A)_{i,j} = O_{\text{pol}}(e_{j,N})$ as well, and

$$ \begin{aligned} \Delta_k^{(N)} &= -\det(B)\det(C)(1 - \text{tr}(C^{-1}DB^{-1}A) + O(\|C^{-1}DB^{-1}A\|^2)) \\ &= -\det(B)(\Delta + O_{\text{pol}}(e_{1,N})). \end{aligned} $$

In particular, $\Delta_k^{(N)} \neq 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

$$ \Delta_k^{(N)} = \begin{vmatrix} A_k & B_k \\ C_k & D_k \end{vmatrix}. $$
For $k > 0$, the same reasoning as above (except that $C_k$ may now be singular) leads to

$$ 
\begin{aligned}
\Delta_k^{(N)} &= -\det(B) \det(C_k - DB^{-1}A_k) \\
&= -\det(B)(\det(C_k) + O_{\text{pol}}(e_{1,N})) \\
&= -\det(B)(\Delta_k + O_{\text{pol}}(e_{1,N})),
\end{aligned}
$$

hence

$$ (4.8) \qquad \gamma_k^{(N)} = \frac{\Delta_k^{(N)}}{\Delta^{(N)}} = \gamma_k + O_{\text{pol}}(e_{1,N}), \quad k > 0. $$

In the case $k < 0$, write

$$ \Delta_k^{(N)} = -\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

$$ (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} \sim e_{-j,N+i-1}, \quad j \neq -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}), \quad j = -k. $$

It follows that

$$ \det(B_k + AC^{-1}D_k) = O_{\text{pol}}(e_{-1,N} \cdots \widehat{e_{k,N}} \cdots e_{-s,N}e_{1,N}) = O_{\text{pol}}\left(\frac{\det(B)}{e_{k,N}}e_{1,N}\right), $$

whence

$$ (4.9) \qquad \gamma_k^{(N)} = \frac{\Delta_k^{(N)}}{\Delta^{(N)}} = \frac{-\det(B)\det(C)O_{\text{pol}}(e_{1,N}/e_{k,N})}{-\det(B)\det(C)(1+O(e_{1,N}))} = O_{\text{pol}}\left(\frac{e_{1,N}}{e_{k,N}}\right), \quad k < 0. $$

Combining (4.6), (4.7) with (4.8), (4.9) finally yields

$$ y_n^{(N)} = y_n + O_{\text{pol}}\left(e_{1,N} \sum_{k=1}^{s} \left(e_{k,n} + \frac{e_{-k,n}}{e_{-k,N}}\right)\right) $$

as $N \to \infty$, uniformly in $n$. $\square$

*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)} = \tilde{y}$ returned by Algorithm 4.2 satisfies

$$ (4.10) \qquad \|y^{(N)} - \pi_d(y)\|_{\infty} \le \phi(N)N^{\kappa_1} \alpha_1^N, $$

where $\kappa_1$ and $\alpha_1$ are the asymptotic growth parameters defined in Section 3.2, for some $\phi$ with $\ln \phi(N) = o(N)$. Thus, given $\epsilon > 0$, it suffices to take $N = O(\ln(\epsilon^{-1}))$ in order to obtain $\|y^{(N)} - \pi_d(y)\|_{\infty} \le \epsilon$. The constant hidden in the $O(\cdot)$ depends on $y$. The estimate goes down to $O(\ln(\epsilon^{-1})/\ln(\epsilon^{-1}))$ when $\kappa_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.
**Corollary 4.6.** Let $L$ and $(\ell_k)_{k=1}^r$ be fixed. Given $d \in \mathbb{N}$, there exists $N$ such that Algorithm 4.2, called with parameters $L$, $(\ell_k)$, $d$, and $N$, computes a polynomial $p_d$ of degree at most $d$ satisfying $\|p_d - y\|_{\infty} \le (4\pi^{-2}\ln(d+1)+5) \|p_d^* - y\|_{\infty}$ in $O(\ln \|p_d - y\|_{\infty}^{-1})$ arithmetic operations.

There is a different way of looking at this, starting with the lower bound [12, Sec. 4.4, Theorem 5(i)]

$$ (4.11) \qquad \|p_d^* - y\|_{\infty} \ge \frac{\pi}{2} \max_{n>d} |y_n| $$

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) \qquad \max_{k=0}^{r-1} |y_{n+k}| \ge n!^{\kappa_1} |\alpha_1|^n \psi(n), \quad \psi(n) = e^{o(n)}, $$

we see by comparing with (4.10) that choosing $N = d + o(d)$ is enough to get

$$ \|y^{(N)} - \pi_d(y)\|_{\infty} \le \|p_d^* - y\|_{\infty}. $$

This last inequality in turn implies

$$ \|y^{(N)} - y\|_{\infty} \le (4\pi^{-2}\ln(d+1)+5) \|p_d^* - y\|_{\infty}. $$

When, in (4.12), $\kappa_1$ and $\alpha_1$ are replaced by $\kappa_i$ and $\alpha_i$ for some $i > 1$, the estimate $N = d + o(d)$ still holds in the case where $\kappa_i = \kappa_1 < 0$. It becomes $N = O(d)$ when either $\kappa_i = \kappa_1 = 0$ or $\kappa_i < \kappa_1 < 0$, and $N = O(d \ln d)$ in general (that is, when $\kappa_i < \kappa_1 = 0$).

**Remark 4.7.** Assuming only $e_{i,n+1}/e_{i,n} \sim \alpha_i n^{\kappa_i} e_{i,n}$ instead of full asymptotic expansions in Lemma 3.4, but with the additional hypothesis $\kappa_i = \kappa_j \Rightarrow \alpha_i \ne \alpha_j$, one can prove that

$$ C(n) \sim_{n \to \infty} e_{0,n} e_{1,n+1} \cdots e_{s-1,n+s-1} \prod_{\substack{i<j \\ \kappa_i = \kappa_j}} \left( \frac{\alpha_i}{\alpha_j} - 1 \right), $$

if the $e_i$ are sorted so that $\kappa_0 \le \kappa_1 \le \cdots \le \kappa_{s-1}$. This leads to a weaker variant of Theorem 4.4 that does not rely on the Birkhoff-Trjitzinsky theorem.

**4.3. 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 [3] 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.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Plot of the error $p_d - y$ between a degree-$d$ approximation $p_d$ computed by Algorithm 4.2 and the exact solution, for each of the problems listed in Section 4.3 and for $d \in \{30, 60, 90\}$.</figcaption></figure>

(i) The first example is adapted from Kaucher and Miranker [30, p. 222]. It concerns the hyperexponential function

$$y(x) = \frac{e^{x/2}}{\sqrt{x+16}},$$

which can be defined by the differential equation

$$2(x+16)y'(x) - (x+15)y(x) = 0, \quad y(0) = \frac{1}{4}.$$

(ii) Next, we consider the fourth order initial value problem (taken from Geddes [24, p. 31])

$$y^{(4)}(x) - y(x) = 0, \quad y(0) = -y''(0) = \frac{3}{2}, \quad -y'(0) = y'''(0) = \frac{1}{2},$$

with the exact solution

$$y(x) = \frac{3}{2} \cos(x) - \frac{1}{2} \sin(x).$$
(iii) Finally, the second-order differential equation

$$ (2x^2 + 1)y''(x) + 8xy'(x) + (2x^2 + 5)y(x), \quad y(0) = 1, y'(0) = 0, $$

has complex singular points at $z = \pm i/\sqrt{2}$, relatively close to the interval $[-1, 1]$, and admits the exact solution

$$ y(x) = \frac{\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 near-uniform 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. 35) gives numerical values of $\|p - y\|_{\infty}$ and $\|p^* - y\|_{\infty}$ in each case. We will later extend these examples to include in the comparison the bounds on $\|p - y\|_{\infty}$ 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 [33, 34]. It has been observed by Fox [22] and later in greater generality (and different language) by El Daou, Ortiz and Samara [19] 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 \cdot 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) \cap C[x] = \{0\}$.

In a nutshell, the tau method works as follows. The first step is to compute $L \cdot p$ where $p$ is a polynomial of degree $d$ with indeterminate coefficients. Since $(\ker L) \cap C[x] = \{0\}$, the result has degree greater than $d$. One then introduces additional unknowns $\tau_{d+1}, \dots, \tau_{d+m}$ in such number that the system

$$ (4.13) \qquad \begin{cases} L \cdot p = \tau_{d+1} T_{d+1} + \dots + \tau_{d+m} T_{d+m} \\ \lambda_i(p) = \ell_i \end{cases} \qquad (1 \le i \le r) $$

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 $\pi_d(L \cdot y) = 0$ of the original differential equation.

Now let $p = \sum_{n=-d}^{d} p_n T_n$ and extend the sequence $(\tau_n)$ by putting $\tau_n = 0$ for $|n| \notin \lfloor d+1, d+m \rfloor$ and $\tau_{-n} = \tau_n$. It follows from (4.13) that $P \cdot (p_n) = \frac{1}{2} Q \cdot (\tau_n)$ where $P$ and $Q$ are the recurrence operators given by Theorem 2.1. Denoting $\text{Supp } u = \{|n| : u_n \neq 0\}$, we also see from the explicit expression of $Q$ that $\text{Supp}(Q \cdot \tau) \subset \lfloor d, d+m+1 \rfloor$. 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 $\tilde{y}$ computed in Algorithm 4.2 with $N = d$, and let $v = \sum_n v_n T_n = L \cdot \tilde{y}$. We have $P \cdot \tilde{y} = Q \cdot v$ by Theorem 2.1. But the definition of $\tilde{y}$ in the algorithm also implies that $(P \cdot \tilde{y})_n = 0$ when $|n| \le N - s$ (since the $\tilde{y}_n$, $|n| \le N$ are linear combinations of sequences $(f_{i,n})_{|n| \le N}$ recursively
computed using the recurrence $P \cdot f_i = 0$) or $|n| > N + s$ (since $\tilde{y}_n = 0$ for $|n| > N$), so that $\text{Supp}(Q \cdot v) \subset [N - s, N + s - 1]$. It can be checked that the Chebyshev recurrence associated to $L = (\frac{d}{dx})^r$ is $P = \delta_r(n)$: indeed, in the language of [4], it must be the first element of a pair $(P_1, Q_1)$ satisfying $Q_1^{-1}P_1 = I^{-r}$. Thus $\delta_r(n) \cdot u = Q \cdot v$ is equivalent to $u^{(r)} = v$, whence

$$ (4.14) \qquad v(x) = \frac{d^r}{dx^r} \sum_{|n|>r} \frac{(P \cdot \tilde{y})_n}{\delta_r(n)} T_n(x) = \sum_{N-s \le |n| < N+s} \frac{(P \cdot \tilde{y})_n}{\delta_r(n)} T_n^{(r)}(x). $$

We see that the output $\tilde{y}(x)$ of Algorithm 4.2 satisfies an inhomogeneous differential equation of the form $L \cdot \tilde{y} = \tau_{N-s} T_{N-s}^{(r)}(x) + \dots + \tau_{N+s-1} T_{N+s-1}^{(r)}(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.

**Proposition 4.8.** Assume that Equation 1.1 has no polynomial solution. The expression on the monomial basis of the polynomial $\tilde{y}(x)$ returned by Algorithm 4.2 with $N = d$ can be computed in $O(d)$ arithmetic operations, all other parameters being fixed.

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 [48, 8].

*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 \cdot 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 $\tilde{y}$ satisfies $L \cdot \tilde{y} = v$ where $v$ is given by (4.14). The coefficients $(P \cdot \tilde{y})_n/\delta_r(n)$ of (4.14) are easy to compute from the last few Chebyshev coefficients of $\tilde{y}$. One deduces the coefficients $v_n$ in linear time by applying repeatedly the non-homogeneous recurrence relation

$$ (4.15) \qquad 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 $\tilde{y}$ on the monomial basis using the recurrence relation they satisfy. $\square$

## 5. 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 = \frac{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 [18, 40]. Indeed, the sequence $(y_n)_{n \in \mathbb{N}}$ 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.

**5.1. Recurrence and Explicit Expression.** Let $y(x) = a(x)/b(x) \in \mathbb{Q}[x]$ be a rational function with no pole in $[-1, 1]$. As usual, we denote by $(y_n)_{n \in \mathbb{Z}}$, $(a_n)_{n \in \mathbb{Z}}$ and $(b_n)_{n \in \mathbb{Z}}$ the symmetric Chebyshev coefficient sequences of $y$, $a$ and $b$.

**Proposition 5.1.** The Chebyshev coefficient sequence $(y_n)_{n \in \mathbb{Z}}$ obeys the recurrence relation with constant coefficients $b(\frac{1}{2}(S + S^{-1})) \cdot (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

$$ \sum_{i=-\deg b}^{\deg b} b_i z^i \sum_{n=-\infty}^{\infty} y_n z^n = \sum_{n=-\infty}^{\infty} \left( \sum_{i=-\infty}^{\infty} b_i y_{n-i} \right) z^n = \sum_{n=-\infty}^{\infty} a_n z^n, \quad x = \frac{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) \qquad \hat{y}(z) = y\left(\frac{z+z^{-1}}{2}\right) = \sum_{n=-\infty}^{\infty} y_n z^n, \quad \rho^{-1} < |z| < \rho, $$

using partial fraction decomposition. From the computational point of view, it is better to start with the squarefree factorization of the denominator of $\hat{y}$:

$$ (5.2) \qquad \beta(z) = z^{\deg b} b\left(\frac{z+z^{-1}}{2}\right) = \beta_1(z)\beta_2(z)^2 \cdots \beta_k(z)^k $$

and write the full partial fraction decomposition of $\hat{y}(z)$ in the form

$$ (5.3) \qquad \hat{y}(z) = q(z) + \sum_{i=1}^{k} \sum_{\beta_i(\zeta)=0} \sum_{j=1}^{i} \frac{h_{i,j}(\zeta)}{(\zeta - z)^j}, \quad q(z) = \sum_{n} q_n z^n \in \mathbb{Q}[z], $$

where $h_{i,j} \in \mathbb{Q}(z)$. The $h_{i,j}$ may be computed efficiently using the Bronstein-Salvy algorithm [11] (see also [25]).

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

$$ \frac{h_{i,j}(\zeta)}{(\zeta - z)^j} = \frac{(-1)^j z^{-j} h_{i,j}(\zeta)}{(1 - \zeta z^{-1})^j} $$

converges for $|z| > |\zeta|$ and does not contribute to the coefficients of $z^n$, $n \ge 0$ in the complete Laurent series. It follows from the uniqueness of the Laurent expansion
of $\hat{y}$ on the annulus $\rho^{-1} < |z| < \rho$ that<fnref n="5" />

$$ (5.4) \qquad \sum_{n=0}^{\infty} y_n z^n = q(z) + \sum_{i=1}^{k} \sum_{\substack{\beta_i(\zeta) = 0 \\ |\zeta| > 1}} \sum_{j=1}^{i} \frac{h_{i,j}(\zeta)}{(\zeta - z)^j}. $$

We now extract the coefficient of $z^n$ in (5.4) and use the symmetry of $(y_n)_{n \in \mathbb{Z}}$ to get an explicit expression of $y_n$ in terms of the roots of $b(\frac{1}{2}(z + z^{-1}))$.

**Proposition 5.2.** The coefficients of the Chebyshev expansion $y(x) = \sum_n y_{|n|} T_n(x)$ are given by

$$ (5.5) \qquad y_n = q_n + \sum_{i=1}^{k} \sum_{j=1}^{i} \sum_{\substack{\beta_i(\zeta) = 0 \\ |\zeta| > 1}} \binom{n+j-1}{j-1} h_{i,j}(\zeta) \zeta^{-n-j} \quad (n \ge 0) $$

where the $q_n \in \mathbb{Q}$, $\beta_i \in \mathbb{Q}[z]$ and $h_{i,j} \in \mathbb{Q}(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 \in \mathbb{N}}$, instead of $2\deg b$ for that from Proposition 5.1, but now with algebraic instead of rational coefficients in general.

**5.2. Truncation Error.** We can now explicitly bound the error in truncating the Chebyshev expansion of $y$.

**Proposition 5.3.** Let $y \in \mathbb{Q}(x)$ have no pole within the elliptic disk $E_\rho$ (see (2.4)). Assume again the notations from (5.2) and (5.3). For all $d \ge \deg q$, it holds that

$$ \left\| \sum_{n>d} y_n T_n \right\|_{\infty} \le \sum_{i=1}^{k} \sum_{j=1}^{i} \sum_{\substack{\beta_i(\zeta) = 0 \\ |\zeta| > 1}} \frac{|h_{i,j}(\zeta)|(d+2)^{j-1}}{(|\zeta|-1)^j} |\zeta|^{-d-1} = O(d^{\deg b} \rho^{-d}) $$

as $d \to \infty$.

*Proof.* We have $\| \sum_{n>d} y_n T_n \|_{\infty} \le \sum_{n>d} |y_n|$ because $\|T_n\|_{\infty} \le 1$ for all $n$. Using the inequality

$$ \sum_{n>d} \binom{n+j-1}{j-1} t^{n+j} \le (d+2)^{j-1} t^{d+1} \sum_{n=0}^{\infty} \binom{n+j-1}{j-1} t^{n+j} = \frac{(d+2)^{j-1} t^{d+j+1}}{(1-t)^j} $$

$$ \hat{y}(z) = q(z) + \sum_{i=1}^{k} \sum_{\substack{\beta_i(\zeta) = 0 \\ |\zeta| > 1}} \sum_{j=1}^{i} \left( \frac{h_{i,j}(\zeta)}{(\zeta - z)^j} + \frac{h_{i,j}(\zeta^{-1})}{(\zeta^{-1} - z)^j} \right) $$

the Laurent expansion of a single term of the form $\frac{h_{i,j}(\zeta)}{(\zeta - z)^j} + \frac{h_{i,j}(\zeta^{-1})}{(\zeta^{-1} - z)^j}$ is not symmetric for $j > 1$, even if $q(z) = 0$.

<footnote n="5">⁵To prevent confusion, it may be worth pointing out that in the expression</footnote>
for $t < 1$, the explicit expression from Proposition 5.2 yields

$$ 
\begin{aligned}
\sum_{n>d} |y_n| &\le \sum_{n>d} \sum_{i=1}^k \sum_{j=1}^i \sum_{\substack{|\zeta|>1 \\ \beta_i(\zeta)=0}} \binom{n+j-1}{j-1} |h_{i,j}(\zeta)|^{|\zeta|^{-n-j}} \\
&\le \sum_{i=1}^k \sum_{j=1}^i \sum_{\substack{|\zeta|>1 \\ \beta_i(\zeta)=0}} \frac{|h_{i,j}(\zeta)|(d+2)^{j-1}}{(|\zeta|-1)^j} |\zeta|^{-d-1}.
\end{aligned}
$$

Since $|\zeta| > 1$ actually implies $|\zeta| > \rho$ when $b(\frac{1}{2}(\zeta + \zeta^{-1})) = 0$, the asymptotic estimate follows. □

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

**Lemma 5.4.** The product $ab$ where the operands $a, b \in \mathbb{Q}[x]$ and the result are written in the Chebyshev basis may be computed in $O((\deg a)(\deg b))$ operations.

*Proof.* It suffices to loop over the indices $(i, j)$ and, at each step, add to the coefficients of $T_{|i\pm 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_iT_j = T_{i+j} + T_{i-j}$. □

The naïve Euclidean division algorithm [63, 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.

**Lemma 5.5.** The division with remainder $a = bq + r$ ($\deg r < \deg b$) where $a, b, q, r \in \mathbb{Q}[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_m^{-1}a_n x^{n-m}b) < n$ where $a = \sum_i a_i x^i$ and $b = \sum_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_m^{-1}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 \leftarrow a - b_m^{-1}a_n x^{n-m}b$ repeatedly done by the classical algorithm by $a \leftarrow a - 2b_m^{-1}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)$. □

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.

**Proposition 5.7.** Algorithm 5.6 is correct. As $d \to \infty$ and $\epsilon \to 0$ with all other parameters fixed, it runs in $O(d + \ln(\epsilon^{-1}))$ arithmetic operations and returns a polynomial of degree $d' \le \max(d, K \ln(\epsilon^{-1}))$, where $K$ depends on $y$, but not on $f$ or $d$.
**Algorithm 5.6.** Input: a rational fraction $y(x) = a(x)/b(x)$, the Chebyshev coefficients of a polynomial $f = \sum_{n=-d}^{d} f_n T_n(x)$, an error bound $\epsilon$. Output: the Chebyshev coefficients of an approximation $\tilde{y}(x)$ of $fy$ such that $\|\tilde{y} - fy\|_{\infty} \le \epsilon$.

1 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 $\hat{w}(z) = w(x) = r(x)/b(x)$, where $x = \frac{1}{2}(z + z^{-1})$, using the Bronstein-Salvy algorithm

5 find $d' \ge \deg q$ such that $\left\| \sum_{n>d'} y_n T_n \right\| \le \epsilon/4$ using Proposition 5.3

6 compute $\rho_-$ and $\rho_+$ such that $\beta_i(\zeta) = 0 \wedge |\zeta| > 1 \Rightarrow 1 < \rho_- \le |\zeta| \le \rho_+$

7 compute $M \ge \sum_{i=1}^{k} \sum_{j=1}^{i} j(\deg \beta_i) \sup_{\rho_- \le |\zeta| \le \rho_+} (|h'_{i,j}(\zeta)| + |\zeta^{-1}h_{i,j}(\zeta)|) \rho_-^{-j}$

8 set $\epsilon' := \min(\rho_- - 1, M^{-1}(1 - \rho_-^{-1})^{D+1} \frac{\epsilon}{4})$, with $D = \deg b$

9 compute approximations $\tilde{\zeta} \in Q[i]$ of the roots $\zeta$ of $\beta_i$ such that $|\tilde{\zeta} - \zeta| < \epsilon'$

10 for $0 \le n \le d'$

11 set $\tilde{y}_n = q_n + \text{Re} \left( \sum_{i=1}^{k} \sum_{j=1}^{i} \sum_{\substack{\beta_i(\zeta)=0 \\ |\zeta|>1}} \binom{n+j-1}{j-1} h_{i,j}(\zeta) \tilde{\zeta}^{-n-j} \right)$

12 return $\tilde{y}(x) = \sum_{n=-d'}^{d'} \tilde{y}_n T_n(x)$

*Proof.* Firstly, we prove the error bound $\|\tilde{y} - y\|_{\infty} \le \epsilon$. Let $A = \{\zeta : \rho_- \le |\zeta| \le \rho_+\}$ and

$$M_0 = \sup_{\zeta \in A} |h'_{i,j}(\zeta)|, \quad M_1 = \sup_{\zeta \in A} |\zeta^{-1} h_{i,j}(\zeta)|.$$

On $A$, we have

$$\left| \frac{d}{d\zeta} (h_{i,j}(\zeta) \zeta^{-n-j}) \right| \le (M_0 + (n+j)M_1) |\zeta|^{-n-j} \le (n+j)(M_0 + M_1) \rho_-^{-n-j}.$$

By Proposition 5.2, observing that the condition $|\zeta - \tilde{\zeta}| < \rho_- - 1$ from Step 8 implies $[\zeta, \tilde{\zeta}] \subset A$, we have the inequalities

$$\begin{aligned} |y_n - \tilde{y}_n| &\le \sum_{i=1}^{k} \sum_{j=1}^{i} \sum_{\substack{\beta_i(\zeta)=0 \\ |\zeta|>1}} \binom{n+j-1}{j-1} |h_{i,j}(\zeta)\zeta^{-n-j} - h_{i,j}(\tilde{\zeta})\tilde{\zeta}^{-n-j}| \\ &\le \sum_{i=1}^{k} \sum_{j=1}^{i} j(\deg \beta_i) \binom{n+j}{j} (M_0 + M_1) \rho_-^{-n-j} \epsilon' \le M \binom{n+D}{D} \rho_-^{-n} \epsilon'. \end{aligned}$$

Therefore, the output of Algorithm 5.6 satisfies

$$\|y_n - \tilde{y}_n\|_{\infty} \le \sum_{n=-d'}^{d'} |y_n - \tilde{y}_n| + 2 \left\| \sum_{n>d'} y_n T_n \right\| \le \frac{2M\epsilon'}{(1 - \rho_-^{-1})^{D+1}} + 2\frac{\epsilon}{4} \le \epsilon.$$
By Proposition 5.3, for all $\epsilon > 0$, there exists $d' \le K \ln(\epsilon^{-1})$ with $\|\sum_{n>d} y_n T_n\| \le \epsilon$. In addition, we have $\deg q \le d$, hence the degree of the output, as computed in Step 5, satisfies $d' \le \max(d, K \ln(\epsilon^{-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 $\tilde{\zeta}^{-1}$ are computed incrementally. Steps 2 and 3 take $O(d)$ operations by Lemmas 5.4 and 5.5, and do not depend on $\epsilon$. Regarding Step 9, it is known [47, Theorem 1.1(d)] that the roots of a polynomial with integer coefficients can be approximated with of absolute accuracy $\eta$ in $O(\eta^{-1})$ arithmetic operations. (In fact, the bit complexity of the algorithm beyond this statement is also softly linear in $\ln(\epsilon^{-1})$.) Since $M$ depends neither on $\epsilon$ nor on $d$, we have $\epsilon' = \Omega(\epsilon)$, and hence the cost of Step 9 is in $O(\ln(\epsilon^{-1}))$. $\square$

## 6. VALIDATION

Assume that we have computed a polynomial of degree $d$

$$p(x) = \sum_{n=-d}^{d} \tilde{y}_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\|_{\infty} \le 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 [44]. In most cases (with the notable exception of [21, 30, 31]), these methods are based on interval Taylor series expansions such as so-called Taylor models [37, 45]. 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 $\alpha_0, \alpha_1, \dots, \alpha_r \in Q[x]$ be such that

$$L = a_r \partial^r + \dots + a_1 \partial + a_0 = \partial^r \alpha_r + \dots + \partial \alpha_1 + \alpha_0,$$

and, for $k \in [0, r]$, define

$$L_k = \partial^{r-k} \alpha_r + \dots + \partial \alpha_{k+1} + \alpha_k \in Q[x] \langle \partial \rangle.$$

Observe that for all $a \in Q[x]$ and $k \in \mathbb{N}$, the operator $\partial^k a$ is of the form

$$\partial^k a = a \partial^k + * \partial^{k-1} + \dots,$$

so that $\alpha_r = a_r$. In particular, $\alpha_r$ does not vanish on $[-1, 1]$.

**Lemma 6.1.** The solution $y \in X$ of $L \cdot y = 0$ such that $y^{(k)}(0) = \ell_k$ for $k \in [0, r-1]$ satisfies

$$(6.1) \qquad \alpha_r(x)y(x) = g(x) + \int_0^x K(x, t)y(t)dt,$$
where

$$
(6.2) \qquad K(x,t) = - \sum_{k=0}^{r-1} \frac{(x-t)^k}{k!} \alpha_{r-1-k}(t) = \sum_{j=0}^{r-1} \beta_j(t)x^j,
$$
$$
\beta_j(t) = \sum_{i=0}^{r-1-j} \frac{(-1)^{i+1}}{i!j!} t^i \alpha_{r-1-j-i}(t),
$$

and

$$
(6.3) \quad g(x) = \sum_{k=0}^{r-1} (L_{r-k} \cdot y)(0) \frac{x^k}{k!} = \sum_{k=0}^{r-1} \sum_{j=0}^{k} \sum_{i=j}^{k} \binom{i}{j} \alpha_{r-k+i}^{(i-j)} (0) \ell_j \frac{x^k}{k!}
$$

Proof. We have

$$
\int_0^x L_k \cdot y = (L_{k+1} \cdot y)(x) - (L_{k+1} \cdot y)(0) + \int_0^x \alpha_k y.
$$

Applying this formula to integrate $r$ times the equation $L \cdot y = 0$ yields

$$
\begin{align*}
& (\alpha_r y)(x) + \int_0^x (\alpha_{r-1} y)(x_1) dx_1 + \dots + \int_0^x \int_0^{x_1} \dots \int_0^{x_{r-1}} (\alpha_0 y)(x_r) dx_r \dots dx_2 dx_1 \\
& = (L_1 \cdot y)(0) \frac{x^{r-1}}{(r-1)!} + \dots + (L_{r-1} \cdot y)(0) x + (L_r \cdot y)(0),
\end{align*}
$$

and (6.1) follows using the relation

$$
\int_0^x \int_0^{x_1} \dots \int_0^{x_{k-1}} f(x_k) dx_k \dots dx_2 dx_1 = \int_0^x \frac{(x-t)^{k-1}}{(k-1)!} f(t) dt
$$

(which can be obtained by repeated integrations by parts).

<figure><img src="image_4.png" /></figure>

Let $X$ denote the Banach space of continuous functions from $[0,1]$ to $\mathbb{R}$ (or $\mathbb{C}$), equipped with the uniform norm. With $K$ and $g$ as in Lemma 6.1, define an operator $T: X \to X$ by

$$
(6.4) \quad T(f)(x) = \alpha_r(x)^{-1} \left( g(x) + \int_0^x K(x,t)f(t)dt \right),
$$

so that (6.1) becomes $T(y) = y$. 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\|_{\infty}$.

**Proposition 6.2.** *Assume that A is an upper bound on $|\alpha_r(x)^{-1}K(x,t)|$ for $x \in [-1,1]$ and $t$ between 0 and $x$. Then the bound*

$$
\|p - y\|_{\infty} \leq \gamma_i \|T^i(p) - p\|_{\infty}, \quad \gamma_i = \sum_{j=0}^{\infty} \frac{A^{ij}}{(ij)!},
$$

holds for all $i \ge 1$.

Proof. Let $V : f \mapsto \alpha_r^{-1} \int K f$ denote the linear part of $T$. Equation (6.4) rewrites again as $(\mathrm{Id} - V) \cdot y = g$. More generally, we have

$$
(\mathrm{Id} - V^i) \cdot y = T^i(y) - V^i \cdot y = T^i(f) - V^i \cdot f = (\mathrm{Id} + V + \dots + V^{i-1}) \cdot g
$$
for all $f \in X$. The operator $V$ is continuous, and since

$$|V^i \cdot y(x)| \le \int_0^x A \int_0^{x_1} A \cdots \int_0^{x_{i-1}} A \|y\|_{\infty} dx_i \cdots dx_2 dx_1 = \frac{A^i}{i!} \|y\|_{\infty},$$

the (subordinate) norm of $V^i$ satisfies

$$\|V^i\| \le \frac{A^i}{i!}.$$

Therefore, the series $\sum_j V^{ij}$ converges, $(\mathrm{Id} - V^i)$ is invertible, and $\|(\mathrm{Id} - V^i)^{-1}\| \le \gamma_i$. Writing

$$p - y = (\mathrm{Id} - V^i)^{-1} \cdot [(\mathrm{Id} - V^i) \cdot p - (T^i(p) - V^i \cdot p)] = (\mathrm{Id} - V^i)^{-1} \cdot (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 $\gamma_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 $\gamma_i \to 1$ as $i \to \infty$. In particular, we have

$$(6.5) \qquad \gamma_i \le \frac{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 $\gamma_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$.

**6.2. Algorithm.** The actual computation of $\|T^i(p) - p\|_{\infty}$ relies on the following lemmas.

**Lemma 6.3.** One can compute an antiderivative of a polynomial<fnref n="6" /> of degree at most $d$ written on the Chebyshev basis in $O(d)$ arithmetic operations.

*Proof.* If $f = \sum_n c_n T_n$ and $f' = \sum_n c'_n T_n$, then we have $2nc_n = c'_{n-1} + c'_{n+1}$ according to Equation (2.7). □

**Lemma 6.4.** Let $f(x) = \sum_{n=-d}^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 \ge 0$ such that

$$\|f\|_{\infty} \le M \le \sqrt{d+1} \|f\|_{\infty}$$

in $O(d)$ arithmetic operations.

*Proof.* It suffices to take $M = \sum_{n=-d}^d |c_n|$. Indeed, we have

$$\|f\|_2 = \left( \frac{1}{\pi} \int_{-1}^{1} \frac{f(t)^2}{\sqrt{1-t^2}} dt \right)^{1/2} = \left( |c_0|^2 + 4 \sum_{n=1}^{d} |c_n|^2 \right).$$

It follows that

$$\|f\|_2 \le \|f\|_{\infty} \le M \le \sqrt{d+1} \|f\|_2,$$

<footnote n="6">⁶Note however that derivation does not commute with the truncation of Chebyshev series.</footnote>
where the first inequality results from the integral expression of $\|f\|_2$, the second one, from the fact that $\|T_n\|_\infty \le 1$ for all $n$, and the last one from the Cauchy-Schwarz inequality. $\square$

This results in Algorithm 6.5. Again here, the suggested bounds for $\|\alpha_r^{-1}\|_\infty$, $\|p - p_i\|_\infty$ and $\gamma_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\partial^r + \dots + a_1\partial + a_0 \in Q[x]\langle\partial\rangle$ of order $r$ such that $a_r(x) \ne 0$ for $x \in [-1, 1]$, initial values $\ell_0, \ell_1, \dots, \ell_{r-1}$. A degree-$d$ polynomial $p$ written on the Chebyshev basis. An accuracy parameter $\epsilon > 0$. Output: A real number $B > 0$ such that $\|y - p\|_\infty \le B$, where $y$ is the unique solution of $L \cdot y = 0$ satisfying $y(0) = \ell_0, y'(0) = \ell_1, \dots, y^{(r-1)}(0) = \ell_{r-1}$.

1. using the commutation rule $x\partial = \partial x - 1$, compute polynomials $\alpha_0, \alpha_1, \dots, \alpha_r \in Q[x]$ such that $L = \partial^r \alpha_r + \dots + \partial \alpha_1 + \alpha_0$

2. define $K \in Q[t, x], (\beta_j)_{j=0}^{r-1} \in Q[x]^r$ and $g \in Q[x]$ as in Lemma 6.1

3. compute $A \ge \max\{|\alpha_r(x)^{-1}K(x, t)|: 0 \le t \le x \le 1 \text{ or } -1 \le x \le t \le 0\}$ (e.g., using Algorithm 5.6 to expand $\alpha_r(x)^{-1}$ in Chebyshev series), and define $(\gamma_i)_{i=1}^\infty$ as in Proposition 6.2

4. compute the minimum $i$ such that $A^i/i! \le 1/2$

5. set $p_0 = p$

6. for $k = 0, 1, \dots, i-1$

7. compute $q_{k+1}(x) = g(x) + \sum_{j=0}^{r-1} x^j \int_0^x \beta_j(t)p_k(t)dt \in Q[x]$

8. compute $p_{k+1} \in Q[x]$ such that $\|p_{k+1} - \alpha_r(x)^{-1}q_{k+1}(x)\|_\infty \le \epsilon$ using Algorithm 5.6

9. compute $\delta \ge \|p - p_i\|_\infty$ using Lemma 6.4

10. return $B = \gamma_i(\delta + e^A\epsilon)$

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 $\|y - p\|_\infty$. For fixed $L$, as $d \to \infty$ and $\epsilon \to 0$, the bound $B$ satisfies

$$B = O((\|y - p\|_{\infty} + \epsilon)(d + \ln(\epsilon^{-1}))^{1/2})$$

and the algorithm performs $O(d + \ln(\epsilon^{-1}))$ arithmetic operations.

*Proof.* Denote by $V$ the linear part of $T$ and recall from the proof of Proposition 6.2 that $\|V^i\| \le A^i/i!$. For all $k$, the polynomial $p_{k+1}$ computed on line 8 satisfies

$$\|p_{k+1} - T(p_k)\|_{\infty} \le \epsilon,$$

hence we have

$$\begin{align*}
& \|p_i - T^i(p)\|_{\infty} \\
& \le \|p_i - T(p_{i-1})\|_{\infty} + \|T(p_{i-1}) - T^2(p_{i-2})\|_{\infty} + \dots + \|T^{i-1}(p_1) - T^i(p_0)\|_{\infty} \\
& \le \|p_i - T(p_{i-1})\|_{\infty} + \|V\| \|p_{i-1} - T(p_{i-2})\|_{\infty} + \dots + \|V^{i-1}\| \|p_1 - T(p_0)\|_{\infty} \\
& \le e^A \epsilon.
\end{align*}$$
By Proposition 6.2, it follows that

$$ (6.6) \quad \|p-y\|_{\infty} \le \gamma_i \|p-T^i(p)\|_{\infty} \le \gamma_i (\|p-p_i\|_{\infty} + \|p_i-T^i(p)\|_{\infty}) \le \gamma_i (\delta + e^A \epsilon). $$

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) \qquad \delta \le \sqrt{D} \|p - p_i\|_{\infty} $$

where

$$ \begin{align} \|p - p_i\|_{\infty} &\le \|p - y\|_{\infty} + \|y - T^i(p)\|_{\infty} + \|T^i(p) - p_i\|_{\infty} \nonumber \\ (6.8) \qquad &\le (1 + \|V^i\|) \|p - y\|_{\infty} + e^A \epsilon \nonumber \\ &\le \frac{3}{2} \|p - y\|_{\infty} + e^A \epsilon. \end{align} $$

Looking at the definition of $p_{k+1}$ in step 7, we see that $\deg q_{k+1} \le \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} \le \max(\deg q_k, C_2 \ln(\epsilon^{-1}))$. It follows by induction that $\deg p_k \le \max(d, C_2 \ln(\epsilon^{-1})) + C_1k$ for all $k$, whence

$$ (6.9) \qquad D \le \max(d, C_2 \ln(\epsilon^{-1})) + C_1 i + 1. $$

Plugging (6.8) and (6.9) into (6.7) yields the estimate

$$ (6.10) \qquad \delta \le \sqrt{D} \left( \frac{3}{2} \|p - y\|_{\infty} + e^A \epsilon \right) = O\left( (\|y - p\| + \epsilon)(d + \ln(\epsilon^{-1}))^{1/2} \right) $$

and the result then follows from the definition of $\delta$ since $\gamma_i \le 2$.

Finally, the only steps whose cost depends on $d$ or $\epsilon$ 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(\epsilon^{-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. $\square$

Another way to put this is to say that Algorithm 6.5 can be modified to provide an enclosure of $\|p - y\|_{\infty}$. Indeed, Equations (6.7) and (6.8) imply

$$ (6.11) \qquad \|p - y\|_{\infty} \ge b = \frac{2}{3} \left( \frac{\delta}{\sqrt{D}} - e^A \epsilon \right), $$

and this $b$ is a computable lower bound for $\|p - y\|_{\infty}$. Furthermore, using (6.6), we have $\delta \ge \gamma_i^{-1} \|p - y\| - e^A \epsilon$, and hence

$$ b \ge \frac{1}{3\sqrt{D}} (\|p - y\|_{\infty} - (\sqrt{D} + 2)e^A \epsilon). $$

Comparing with the upper bound on $B$ resulting from (6.10), we deduce

$$ \frac{B}{b \ln(b^{-1})} \le \frac{9D \|p - y\|_{\infty} + 6(D + \sqrt{D})e^A \epsilon}{(\|p - y\|_{\infty} - (\sqrt{D} + 2)e^A \epsilon) \ln(\|p - y\|_{\infty}^{-1})}. $$

In particular, if we restrict ourselves to polynomials $p$ satisfying $\|p - y\|_{\infty} \le e^{-\Gamma d}$ for some fixed $\Gamma$, and if $\epsilon$ is chosen such that

$$ (6.12) \qquad \|p - y\|_{\infty}^{2E} \le \epsilon \le \|p - y\|_{\infty}^{E}, \quad E > 1, $$
<table><caption>TABLE 2. 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}/\sqrt{x+16}$, (ii) $y(x) = \frac{3}{2}\cos x + \frac{1}{2}\sin x$, and (iii) $y(x) = (\cos x)/(2x^2 + 1)$.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">d</th><th rowspan="2">ε</th><th rowspan="2">enclosure of ||y - p||<sub>∞</sub><br>computed by Algo. 6.5</th><th rowspan="2">||y - p||<sub>∞</sub></th><th rowspan="2">||y - p*||<sub>∞</sub></th><th colspan="4">time (s)</th></tr><tr><th>D</th><th>i</th><th>4.2</th><th>6.5</th></tr></thead><tbody><tr><td rowspan="3">(i)</td><td>30</td><td>10<sup>-104</sup></td><td>[2.3 · 10<sup>-53</sup>, 4.3 · 10<sup>-52</sup>]</td><td>3.4 · 10<sup>-52</sup></td><td>3.4 · 10<sup>-52</sup></td><td>102</td><td>2</td><td>0.05</td><td>0.54</td></tr><tr><td>60</td><td>10<sup>-194</sup></td><td>[9.0 · 10<sup>-99</sup>, 2.4 · 10<sup>-97</sup>]</td><td>2.0 · 10<sup>-97</sup></td><td>1.9 · 10<sup>-97</sup></td><td>192</td><td>2</td><td>0.05</td><td>1.07</td></tr><tr><td>90</td><td>10<sup>-284</sup></td><td>[4.6 · 10<sup>-144</sup>, 1.5 · 10<sup>-142</sup>]</td><td>1.2 · 10<sup>-142</sup></td><td>1.1 · 10<sup>-142</sup></td><td>282</td><td>2</td><td>0.06</td><td>1.87</td></tr><tr><td rowspan="3">(ii)</td><td>30</td><td>10<sup>-88</sup></td><td>[6.0 · 10<sup>-45</sup>, 9.8 · 10<sup>-44</sup>]</td><td>5.9 · 10<sup>-44</sup></td><td>5.6 · 10<sup>-44</sup></td><td>42</td><td>3</td><td>0.06</td><td>0.06</td></tr><tr><td>60</td><td>10<sup>-206</sup></td><td>[6.7 · 10<sup>-104</sup>, 1.5 · 10<sup>-102</sup>]</td><td>8.8 · 10<sup>-103</sup></td><td>8.5 · 10<sup>-103</sup></td><td>72</td><td>3</td><td>0.07</td><td>0.10</td></tr><tr><td>90</td><td>10<sup>-334</sup></td><td>[2.0 · 10<sup>-169</sup>, 5.1 · 10<sup>-168</sup>]</td><td>3.1 · 10<sup>-168</sup></td><td>3.0 · 10<sup>-168</sup></td><td>102</td><td>3</td><td>0.08</td><td>0.23</td></tr><tr><td rowspan="3">(iii)</td><td>30</td><td>10<sup>-18</sup></td><td>[1.2 · 10<sup>-10</sup>, 2.4 · 10<sup>-9</sup>]</td><td>1.6 · 10<sup>-9</sup></td><td>1.1 · 10<sup>-9</sup></td><td>79</td><td>3</td><td>0.05</td><td>0.74</td></tr><tr><td>60</td><td>10<sup>-36</sup></td><td>[2.2 · 10<sup>-19</sup>, 6.1 · 10<sup>-18</sup>]</td><td>4.1 · 10<sup>-18</sup></td><td>3.0 · 10<sup>-18</sup></td><td>151</td><td>3</td><td>0.06</td><td>1.6</td></tr><tr><td>90</td><td>10<sup>-54</sup></td><td>[4.8 · 10<sup>-28</sup>, 1.7 · 10<sup>-26</sup>]</td><td>1.1 · 10<sup>-26</sup></td><td>7.7 · 10<sup>-27</sup></td><td>223</td><td>3</td><td>0.10</td><td>2.7</td></tr></tbody></table>

then $D \le \max(\Gamma^{-1}, 2EC_2) \ln \|p - y\|_{\infty}^{-1} + O(1)$ as $\|p - y\|_{\infty}$ tends to zero, so there exists $K$ (computable as a function of $y$, $\Gamma$, and $E$) such that

$$ (6.13) \qquad B \le K b \ln(b^{-1}) \le K \|p - y\|_{\infty} \ln(\|p - y\|_{\infty}^{-1}) $$

for small $\|p - y\|_{\infty}$. Of course, since $\|p - y\|_{\infty}$ is what we want to estimate, we do not know the “correct” choice of $\epsilon$ beforehand. But, assuming $\|p - y\|_{\infty}$ is indeed small enough, we can search for a suitable $\epsilon$ iteratively, starting, say, with $\epsilon = 2^{-d}$ and checking whether (6.13) holds at each step. As our hypotheses imply $d = O(\ln \|p - y\|_{\infty}^{-1})$, the whole process requires at most $O(\ln \|p - y\|_{\infty}^{-1})$ operations.

By combining these tightness guarantees with Corollary 4.6 and lower bounds on $\|p_d^* - y\|_{\infty}$ 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\|_{\infty}$.

**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 $\epsilon$ was manually set to approximately $\|p - y\|_{\infty}^2$ based on a heuristic estimate of $\|p - y\|_{\infty}$. In practice, this makes the term $e^A\epsilon$ in the error bound (6.6) small, so that the main contribution to the error bound $B$ in practice is $\|p_i - p\|_{\infty}$.

Besides the upper bound $B$, the table gives a lower bound $b \le \|p - y\|_{\infty}$ obtained as discussed above. For comparison, we include the “true” value of $\|p - y\|_{\infty}$, as well as the error $\|y - p^*\|_{\infty}$ corresponding to the minimax polynomial of degree $d$, computed using Sollya [14].

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 $\epsilon$ 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. Nevertheless, 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\|_{\infty} \le \Lambda(d) \|p_d^* - y\|_{\infty} $$

where $\Lambda(d)$ does not depend on $y$. It is then natural to ask for polynomial approximations where $\|p_d^* - y\|_{\infty}$ not only satisfies the above inequality, but also comes with an explicit upper bound satisfying a similar inequality, that is

$$ (6.14) \qquad \|p_d - y\|_{\infty} \le B \le \Lambda(d) \|p_d^* - y\|_{\infty}, \quad \Lambda \text{ independent of } y. $$

We thus leave open the following question.

**Question 6.8.** 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 \le d$ satisfying (6.14) for some $\Lambda(d)$? For instance, can it be done in $O(d) + \ln(\|p_d^* - y\|_{\infty}^{-1})$ arithmetic operations when $y$ is fixed?

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 [46] 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.
## 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.

## REFERENCES

[1] C. R. Adams. On the irregular cases of the linear ordinary difference equation. *Transactions of the American Mathematical Society*, 30(3):507–541, 1928.

[2] W. Balser and T. Bothner. Computation of formal solutions of systems of linear difference equations. *Advances in Dynamical Systems and Applications*, 5(1):29–47, 2010.

[3] A. Benoit. *Algorithmique semi-numérique rapide des séries de Tchebychev*. Thèse de doctorat, École polytechnique, 2012.

[4] A. Benoit and B. Salvy. Chebyshev expansions for solutions of linear differential equations. In J. P. May, editor, *ISSAC '09*, page 23–30. ACM, 2009.

[5] W. G. Bickley, L. J. Comrie, J. C. Miller, D. H. Sadler, and A. J. Thompson. *Bessel functions. Part II. Functions of positive integer order*, volume X of *Mathematical Tables*. British Association for the Advancement of Science, 1952.

[6] G. D. Birkhoff. Formal theory of irregular linear difference equations. *Acta Mathematica*, 1930.

[7] G. D. Birkhoff and W. J. Trjitzinsky. Analytic theory of singular difference equations. *Acta Mathematica*, 60:1–89, 1933.

[8] A. Bostan, B. Salvy, and Éric Schost. Power series composition and change of basis. In J. R. Sendra and L. González-Vega, editors, *ISSAC '08*, page 269–276. ACM, 2008.

[9] R. P. Brent and P. Zimmermann. *Modern Computer Arithmetic*. Cambridge University Press, 2010.

[10] N. Brisebarre and M. Joldeş. Chebyshev interpolation polynomial-based tools for rigorous computing. In S. M. Watt, editor, *ISSAC '10*, page 147–154. ACM, 2010.

[11] M. Bronstein and B. Salvy. Full partial fraction decomposition of rational functions. In M. Bronstein, editor, *ISSAC '93*, page 157–160. ACM, 1993.

[12] E. W. Cheney. *Introduction to approximation theory*. American Mathematical Society, 1998.

[13] S. Chevillard, J. Harrison, M. Joldeş, and C. Lauter. Efficient and accurate computation of upper bounds of approximation errors. *Theoretical Computer Science*, 16(412):1523–1543, 2011.

[14] S. Chevillard, M. Joldeş, and C. Lauter. Sollya: An environment for the development of numerical codes. In K. Fukuda, J. van der Hoeven, M. Joswig, and N. Takayama, editors, *Mathematical Software - ICMS 2010*, volume 6327 of *Lecture Notes in Computer Science*, pages 28–31, Heidelberg, Germany, September 2010. Springer.

[15] C. W. Clenshaw. A note on the summation of Chebyshev series. *Mathematics of Computation*, 9:118–120, 1955.

[16] C. W. Clenshaw. The numerical solution of linear differential equations in Chebyshev series. *Proceedings of the Cambridge Philosophical Society*, 53(1):134–149, 1957.

[17] T. A. Driscoll, F. Bornemann, and L. N. Trefethen. The chebop system for automatic solution of differential equations. *BIT Numerical Mathematics*, 48(4):701–723, 2008.

[18] T. H. Einwohner and R. J. Fateman. A MACSYMA package for the generation and manipulation of Chebyshev series. In G. H. Gonnet, editor, *ISSAC '89*, page 180–185. ACM, 1989.

[19] M. K. El-Daou, E. L. Ortiz, and H. Samara. A unified approach to the tau method and Chebyshev series expansion techniques. *Computers & Mathematics with Applications*, 25(3):73–82, 1993.

[20] C. Epstein, W. L. Miranker, and T. J. Rivlin. Ultra-arithmetic I: function data types. *Mathematics and Computers in Simulation*, 24(1):1–18, 1982.

[21] C. Epstein, W. L. Miranker, and T. J. Rivlin. Ultra-arithmetic II: intervals of polynomials. *Mathematics and Computers in Simulation*, 24(1):19–29, 1982.

[22] L. Fox. Chebyshev methods for ordinary differential equations. *The Computer Journal*, 4(4):318, 1962.
[23] L. Fox and I. B. Parker. *Chebyshev polynomials in numerical analysis*. Oxford University Press, 1968.

[24] K. O. Geddes. Symbolic computation of recurrence equations for the Chebyshev series solution of linear ODE's. In *Proceedings of the 1977 MACSYMA User's Conference*, page 405-423, July 1977.

[25] X. Gourdon and B. Salvy. Effective asymptotics of linear recurrences with rational coefficients. *Discrete Mathematics*, 153(1-3):145-163, 1996.

[26] A. O. Guelfond. *Calcul des différences finies*. Collection Universitaire de Mathématiques, XII. Dunod, Paris, 1963. Traduit par G. Rideau.

[27] G. K. Immink. Reduction to canonical forms and the Stokes phenomenon in theory of linear difference equations. *SIAM Journal on Mathematical Analysis*, 22(1):238-259, 1991.

[28] G. K. Immink. On the relation between global properties of linear difference and differential equations with polynomial coefficients. II. *Journal of Differential Equations*, 128(1):168-205, 1996.

[29] M. Joldes. *Rigorous polynomial approximations and applications*. Thèse de doctorat, École normale supérieure de Lyon, 2011.

[30] E. W. Kaucher and W. L. Miranker. *Self-validating numerics for function space problems*. Academic Press, 1984.

[31] E. W. Kaucher and W. L. Miranker. Validating computation in a function space. In *Reliability in computing: the role of interval methods in scientific computing*, page 403-425. Academic Press Professional, Inc., 1988.

[32] M. Kauers. The holonomic toolkit. *Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions, Texts and Monographs in Symbolic Computation*, 2013.

[33] C. Lanczos. Trigonometric interpolation of empirical and analytical functions. *Journal of Mathematical Physics*, 17:123-199, 1938.

[34] C. Lanczos. *Applied analysis*. Prentice-Hall, 1956.

[35] S. Lewanowicz. Construction of a recurrence relation of the lowest order for coefficients of the Gegenbauer series. *Zastosowania Matematyki*, XV(3):345-395, 1976.

[36] S. Lewanowicz. A new approach to the problem of constructing recurrence relations for the Jacobi coefficients. *Zastosowania Matematyki*, 21:303-326, 1991.

[37] K. Makino and M. Berz. Taylor models and other validated functional inclusion methods. *International Journal of Pure and Applied Mathematics*, 4(4):379-456, 2003.

[38] Maplesoft (Waterloo Maple, Inc.). Maple, 1980-2014.

[39] J. C. Mason and D. C. Handscomb. *Chebyshev polynomials*. CRC Press, 2003.

[40] R. J. Mathar. Chebyshev series expansion of inverse polynomials. *Journal of Computational and Applied Mathematics*, 196(2):596-607, 2006.

[41] H. Meschkowski. *Differenzengleichungen*. Vandenhoeck & Ruprecht, 1959.

[42] M. Mezzarobba. *Autour de l'évaluation numérique des fonctions D-finies*. Thèse de doctorat, École polytechnique, Nov. 2011.

[43] L. M. Milne-Thomson. *The calculus of finite differences*. Macmillan, London, 1933.

[44] R. E. Moore. *Methods and applications of interval analysis*. Society for Industrial and Applied Mathematics, 1979.

[45] A. Neumaier. Taylor forms – Use and limits. *Reliable Computing*, 9(1):43–79, 2003.

[46] S. Olver and A. Townsend. A fast and well-conditioned spectral method. *SIAM Review*, 55(3):462-489, 2013.

[47] V. Y. Pan. Optimal and nearly optimal algorithms for approximating polynomial zeros. *Computers & Mathematics with Applications*, 31(12):97-138, 1996.

[48] V. Y. Pan. New fast algorithms for polynomial interpolation and evaluation on the Chebyshev node set. *Computers & Mathematics with Applications*, 35(3):125-129, 1998.

[49] S. Paszkowski. Zastosowania numeryczne wielomianow i szeregow Czebyszewa. *Podstawowe Algorytmy Numeryczne*, 1975.

[50] H. Poincaré. Sur les équations linéaires aux différentielles ordinaires et aux différences finies. *American Journal of Mathematics*, 7(3):203-258, 1885.

[51] L. B. Rall. *Computational solution of nonlinear operator equations*. With an appendix by Ramon E. Moore. John Wiley & Sons Inc., New York, 1969.

[52] L. Rebillard. *Etude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes*. Thèse de doctorat, Institut national polytechnique de Grenoble, 1998.
[53] T. J. Rivlin. *The Chebyshev polynomials*. Wiley, 1974.

[54] B. Salvy. D-finiteness: Algorithms and applications. In M. Kauers, editor, *ISSAC '05*, page 2–3. ACM, 2005. Abstract for an invited talk.

[55] F. W. Schäfke. Lösungstypen von Differenzengleichungen und Summengleichungen in normierten abelschen Gruppen. *Mathematische Zeitschrift*, 88(1):61–104, Feb. 1965.

[56] R. P. Stanley. Differentiably finite power series. *European Journal of Combinatorics*, 1(2):175–188, 1980.

[57] É. Tournier. *Solutions formelles d'équations différentielles*. Doctorat d'État, Université scientificque, technologique et médicale de Grenoble, 1987.

[58] L. N. Trefethen. Computing numerically with functions instead of numbers. *Mathematics in Computer Science*, 1(1):9–19, 2007.

[59] L. N. Trefethen. *Approximation theory and approximation practice*. SIAM, 2013.

[60] W. Tucker. *Validated Numerics: A Short Introduction to Rigorous Computations*. Princeton University Press, 2011.

[61] H. L. Turrittin. The formal theory of systems of irregular homogeneous linear difference and differential equations. *Boletín de la Sociedad Matemática Mexicana*, 5:255–264, 1960.

[62] M. van der Put and M. F. Singer. *Galois theory of difference equations*, volume 1666 of *Lecture Notes in Mathematics*. Springer, Berlin, 1997.

[63] J. von zur Gathen and J. Gerhard. *Modern Computer Algebra*. Cambridge University Press, 2nd edition, 2003.

[64] J. Wimp. On recursive computation. Technical Report ARL 69-0186, Aerospace Research Laboratories, 1969.

[65] J. Wimp. *Computation with Recurrence Relations*. Pitman, Boston, 1984.

[66] R. V. M. Zahar. A mathematical analysis of Miller's algorithm. *Numerische Mathematik*, 27(4):427–447, 1976.

[67] D. Zeilberger. A holonomic systems approach to special functions identities. *Journal of Computational and Applied Mathematics*, 32(3):321–368, 1990.

ALEXANDRE BENOIT; ÉDUCATION NATIONALE; FRANCE

E-mail address: alexandrebenoit@yahoo.fr

URL: http://alexandre.benoit.83.free.fr/

MIOARA JOLDES; CNRS, LAAS, 7 AVENUE DU COLONEL ROCHE, 31077 TOULOUSE, CEDEX 4,
FRANCE

E-mail address: joldes@laas.fr

URL: http://homepages.laas.fr/mmjoldes/

MARC MEZZAROBBA; CNRS, UMR 7606, LIP6, F-75005, PARIS, FRANCE; SORBONNE UNIVERSITÉS, UPMC UNIV PARIS 06, UMR 7606, LIP6, F-75005, PARIS, FRANCE

E-mail address: marc@mezzarobba.net

URL: http://marc.mezzarobba.net/