Florent Bréhard
email: florent.brehard@math.uu.se

A Symbolic-Numeric Validation Algorithm for Linear ODEs with Newton-Picard Method

Keywords: Mathematics Subject Classification (2010). 65G20, 65L60, 65L70, 65Y20. Keywords. validated numerics, a posteriori validation, Chebyshev spectral methods, linear differential equations, D-finite functions

A symbolic-numeric validation algorithm is developed to compute rigorous and tight uniform error bounds for polynomial approximate solutions to linear ordinary differential equations, and in particular D-finite functions. It relies on an a posteriori validation scheme, where such an error bound is computed afterwards, independently from how the approximation was built. Contrary to Newton-Galerkin validation methods, widely used in the mathematical community of computer-assisted proofs, our algorithm does not rely on finite-dimensional truncations of differential or integral operators, but on an efficient approximation of the resolvent kernel using a Chebyshev spectral method. The result is a much better complexity of the validation process, carefully investigated throughout this article. Indeed, the approximation degree for the resolvent kernel depends linearly on the magnitude of the input equation, while the truncation order used in Newton-Galerkin may be exponential in the same quantity. Numerical experiments based on an implementation in C corroborate this complexity advantage over other a posteriori validation methods, including Newton-Galerkin.

Introduction

Many functions arising in scientific computing are implicitly defined by functional equations, of which differential equations are probably the most commonly encountered ones. Yet, one frequently requires to explicitly manipulate and evaluate these functions. Polynomials (in a broad sense) often provide accurate, global and smooth approximants, and are easy to manipulate on a computer. Moreover, in some applicative contexts, like computer-assisted proofs in mathematics, notably in dynamical systems, or safety-critical industrial applications (e.g., aerospace challenges, automatic vehicles or remote medicine), those computations must be accompanied with strong reliability guarantees.

Background and Related Works

A wide range of techniques in computer algebra allow for a symbolic, exact manipulation of certain classes of functions defined implicitly. A major example are the D-finite functions, that is the solutions of linear ordinary differential equations (ODEs) with polynomial coefficients [START_REF] Stanley | Differentiably finite power series[END_REF]. Despite this seemingly restricted definition, D-finite functions account for about 60% [START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF] of the functions listed in Abramowitz and Stegun's Handbook of Mathematical Functions [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. Yet, many applications require to move from symbolic representations, based on linear differential operators [START_REF] Salvy | Linear differential equations as a data-structure[END_REF], to numeric or symbolicnumeric ones providing quantitative information. While numerical analysis focuses on algorithms to compute approximations, error estimates are most of the time asymptotic. For example, the Chebfun library 1 in Matlab allows for powerful numerical computations with functions using polynomials obtained by truncating Chebyshev expansions, but it does not provide rigorous error bounds. On the other side, validated numerics [START_REF] Tucker | Validated Numerics: A Short Introduction to Rigorous Computations[END_REF] (also called rigorous or set-valued numerics) is concerned with the computation of such error bounds, encompassing both method and rounding errors. Such a challenge, and more specifically the following problem, is at the core of the present work.

Problem 1.1. Given a polynomial y • approximating the solution y * of the Initial Value Problem (IVP):

y (r) (x) + a r-1 (x)y (r-1) (x) + • • • + a 0 (x)y(x) = h(x),
x ∈ [x l , x r ],

y (i) (x 0 ) = v i ∈ R, 0 i < r, x 0 ∈ [x l , x r ], (1)
compute a bound ε y • -y * ∞ = sup x∈[x l ,xr] |y • (x) -y * (x)|. The branch of validated numerics dedicated to the validation of functional problems has its origins in the early 80s with ultra-arithmetics [START_REF] Epstein | Ultra-arithmetic. I. Function data types[END_REF][START_REF] Epstein | Ultra-arithmetic. II. Intervals of polynomials[END_REF][START_REF] Kaucher | Validating computation in a function space[END_REF]. By analogy with floating-point numbers and interval arithmetics [START_REF] Moore | Interval Analysis[END_REF] to rigorously approximate real numbers on computers, ultra-arithmetics proposes a quite general framework to provide similar operations on functions. Solutions are approximated using polynomials (in a broad sense), and rigorously enclosed in balls or other domains in function spaces using rigorous polynomial approximations (RPAs), typically pairs (p, ε) with an approximation p and a rigorous error bound ε. Notable implementations of RPAs include Taylor models [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF][START_REF] Neumaier | Taylor forms-use and limits[END_REF] and Chebyshev models [START_REF] Brisebarre | Chebyshev interpolation polynomial-based tools for rigorous computing[END_REF][START_REF] Joldes | Rigorous Polynomial Approximations and Applications[END_REF], based on Taylor and Chebyshev expansions, respectively.

Numerous techniques have been developed over past decades to compute rigorous enclosures for solutions of ODEs [START_REF] Kedem | A posteriori error bounds for two-point boundary value problems[END_REF][START_REF] Plum | Computer-assisted existence proofs for two-point boundary value problems[END_REF][START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF][START_REF] Nakao | Numerical verification methods for solutions of ordinary and partial differential equations[END_REF][START_REF] Zgliczyński | C 1 Lohner algorithm[END_REF][START_REF] Lessard | Rigorous numerics for nonlinear differential equations using Chebyshev series[END_REF][START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach[END_REF] and of wider classes of equations (e.g., partial differential equations, delay equations), thereby contributing to the rise of computer-assisted proofs in dynamics [START_REF] Van Den Berg | Rigorous numerics in dynamics[END_REF]. More recently, several works [START_REF] Bournez | On the complexity of solving initial value problems[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF][START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF][START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] paved the way towards an algorithmic and complexity approach to such methods, aiming at the development of general-purpose reliable and efficient libraries for symbolic-numeric computations. This article brings its own contributions to this approach. The counterpart is the restriction to more limited classes of ODEs, notably D-finite functions, that is, linear ODE [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] with polynomial or rational coefficients a i (x). However, we emphasize the fact that many of the ODE validation methods cited above proceed by linearization of the ODE after bounding the nonlinear terms. Also, analytic coefficients a i (x) in (1) can be rigorously represented with RPAs. Hence, the application scope of these algorithmic works, including this article in particular, goes far beyond D-finite functions: they provide an algorithmic core for the validation of differential equations in a broad perspective.

Setting and Contributions

A powerful representation of D-finite functions y in computer algebra is via differential operators L for which L{y} = 0 (together with sufficiently many initial conditions), where:

L = ∂ r + a r-1 (x)∂ r-1 + • • • + a 1 (x)∂ + a 0 (x),
with a i polynomials or rational functions, belongs to the Ore algebra of differential operators [START_REF] Chyzak | The ABC of Creative Telescoping: Algorithms, Bounds, Complexity. Memoir of accreditation to supervise research (HDR)[END_REF][START_REF] Koutschan | Advanced Applications of the Holonomic Systems Approach[END_REF], spanned by x (x{y} = x → xy(x)) and ∂ (∂{y} = y ), with the skew commutation rule ∂x -x∂ = 1. In the theoretical developments of this article, the a i may actually be analytic functions over the real compact interval [x l , x r ] of interest. For the algorithmic parts, they are supposed to be polynomials (with rational, floating-point or interval coefficients) for the sake of clarity, even if adapting the algorithms to RPA representations of the a i (x) is rather straightforward (see the example of Section 4.3). We thus propose a fully algorithmic validation method for Problem 1.1 with analytic [f ]n n-th Chebyshev coefficient of f r order of L пN : C 0 → RN [x] N-th Chebyshev truncation ai coefficients of L Ч 1 Banach space of absolutely summable s = max 0 i<r deg ai maximum degree Chebyshev series h right-hand side of differential equation

• Ч 1 norm associated to Ч 1 y unknown function in differential equation K [N ] truncated integral operator y (i) iterated derivative of y I identity operator y [i] y [0] = y and y [i+1] = ar-1-iy -∂{y [i] } IN order N identity matrix K(x, t) bivariate kernel K * L composition of kernels [x l , xr|x0] = {(x, t) s.t. x l x t x0 R(x, t) resolvent kernel or x0 t x xr} R • (x, t) polynomial approximation of R(x, t) K integral operator of kernel K(x, t)

Eρ Bernstein ellipse of parameter ρ > 1 g right-hand side of integral equation NG truncation order in Newton-Galerkin f unknown function in integral equation N R degree of R • (x, t) in Newton-Picard Table 1. Notation.

coefficients a i representable by RPAs. In particular, iterating this construction allows for the rigorous numerical treatment of D n -finite functions [START_REF] Jiménez-Pastor | Some structural results on D n -finite functions[END_REF]. Also, in the case of non-homogeneous linear ODEs, the same remarks hold for the right-hand side h in [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. Our validation method is based on an automatic reformulation of IVP (1) into an equivalent Volterra linear integral equation of the second kind [START_REF] Lalescu | Introduction à la Théorie des Équations Intégrales (Introduction to the Theory of Integral Equations[END_REF] (simply referred to as integral equation in this article):

f (x) + K{f }(x) = g(x),
where K{f }(x)

= x x0 K(x, t)f (t)dt. (2)
The bivariate kernel K(x, t) associated to the integral operator K is analytic, and even polynomial if the a i are. Therefore, representations through integral equations are equally algorithmic. Moreover, such integral equations are particularly suited for numerical resolution via Galerkin spectral methods [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF][START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF][START_REF] Olver | A fast and well-conditioned spectral method[END_REF], which compute a polynomial approximation y • for Problem 1.1 using truncated expansions (e.g., Taylor, Fourier, Chebyshev), by transforming functional equation ( 2) into a linear system on the unknown coefficients. Several works, e.g. [START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF][START_REF] Lessard | Rigorous numerics for nonlinear differential equations using Chebyshev series[END_REF][START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach[END_REF], proposed ODE validation methods based on a posteriori Newtonlike methods in well-chosen coefficient spaces. The main idea is to approximate the (infinite-dimensional) inverse Jacobian by finite-dimensional truncations, as done by spectral methods in numerical analysis, in order to construct a contracting Newton-like operator and apply the Banach fixed-point theorem or a variant of it. Such methods, which we propose to call Newton-Galerkin validation methods, have been extensively used in the mathematical community of computer-assisted proofs. In [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], we developed a detailed analysis of a Newton-Galerkin validation algorithm in the Chebyshev basis, by enhancing this general scheme with the linear algebra algorithms of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] for so-called almost-banded linear systems. In the present article, we propose an alternative validation algorithm, which we call Newton-Picard, to overcome the complexity shortcomings of Newton-Galerkin validation methods, acknowledged in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. The leitmotif is to approximate the inverse Jacobian by computing polynomial approximations of the so-called resolvent kernel associated to (1) which arises from the well-known Picard iterations. More or less related ideas have already been used in previous works, notably [START_REF] Kedem | A posteriori error bounds for two-point boundary value problems[END_REF][START_REF] Plum | Computer-assisted existence proofs for two-point boundary value problems[END_REF][START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF]. However, the algorithmic description and complexity analysis we propose here lead to a robust and efficient validation algorithm, together with a fair comparison with Newton-Galerkin. Our conclusions are the following.

1. We prove that Newton-Picard runs in polynomial time with respect to the magnitude of the coefficients in the input equation, while Newton-Galerkin's complexity can be exponential in the same parameters, even, surprisingly, in the analytic case we consider here. The result is that moderate to hard instances of Problem 1.1 (e.g., exponentially increasing or highly oscillating functions) may be intractable with Newton-Galerkin, but solved within a few seconds with Newton-Picard on a modern computer. 2. The Newton-Picard validation method is widely independent of the underlying approximation tools, although presented in a Chebyshev approximation theory flavor in this article. Indeed, the algorithms could be parametrized by a generic type of RPAs to rigorously enclose functions without major modifications. Instead of Chebyshev models, one could use Fourier approximations (for trigonometric functions), splines (for C k functions), depending on the context. Validation of solutions with singularities or on unbounded intervals also belongs to future applications of Newton-Picard method. 3. A challenging goal we want to address in the future is the formalization in a proof assistant, for example Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development: Coq'Art: the Calculus of Inductive Constructions[END_REF], of an ODE a posteriori validation algorithm. To this purpose, Newton-Picard turns out to be a better candidate than Newton-Galerkin. Indeed, the correctness of the algorithms presented in this article requires very limited functional analysis background, contrary to Newton-Galerkin which involves more Chebyshev approximation theory. Moreover, the resolvent kernel approximation can be safely computed outside the proof assistant without compromising the final reliability.

Throughout this article, we consider an arithmetic model for the complexity, where each arithmetic operation (+, -, ×, ÷), typically on floating-point numbers or intervals, is assumed to have a constant cost O(1). Therefore, the claimed complexity estimates do not take into account the fact that particularly stiff ODEs may require a larger floating-point bit precision to provide sufficiently accurate results.

Outline

First, Section 2 provides the necessary background (a posteriori Newton-like validation principle, reformulation of IVP (1) into integral equation (2), Chebyshev approximation theory and spectral methods) and summarizes the main lines of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] about the Newton-Galerkin validation method. Then, Section 3 introduces Newton-Picard validation method, with the main algorithm NewtonPicardValid and the approximation routine ResolventKernel for the resolvent kernel approximation, together with detailed complexity results. After that, Section 4 presents how the method is implemented in the C library ChebValid2 and illustrates its performances on practical examples, in comparison with other methods. We finally discuss the obtained results and the future perspectives of this work.

A Posteriori Validation and Newton-Galerkin Method

This section lays out the main prerequisites to this article: a posteriori validation paradigm using Newton-like fixed-point operators (Section 2.1), transformation(s) of linear IVPs (1) into integral equations (2) (Section 2.2), some highlights of Chebyshev approximation theory and in particular spectral methods (Section 2.3). Finally, Section 2.4 gives an overview of the principles and limitations of Newton-Galerkin validation method, thoroughly investigated in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. This serves as a basis for comparison with the novel Newton-Picard algorithm presented in Section 3.

Principle of A Posteriori Validation using Newton-like Fixed-Point Operators

In an a posteriori validation scheme, an approximation x • to the exact solution x * is first constructed by whatever numerical routine of choice, and then a rigorous error bound between x • and x * is computed in a totally independent second step. Some advantages of this paradigm include the possibility to use external non-trusted numerical routines for the approximation step without compromising the final reliability; the need for rigorous numerics (e.g., interval arithmetics) in the validation step only; and a rather elegant mathematical formulation that eases the verification by a formal proof assistant.

A wide range of such methods involve fixed-point based validation, in which x * is expressed as the fixed point of a well chosen equation T(x) = x, with T a contracting operator. Then a rigorous error bound is obtained using a suitable fixed-point theorem, for example this formulation of the Banach fixed-point theorem [3, Thm. 2.1].

Theorem 2.1. Let T : X → X with (X, • ) a Banach space and D ⊆ X closed, such that:

• D is T-stable: T(D) ⊆ D;
• T is contracting over D: there exists λ ∈ [0, 1) such that:

T(x) -T(y) λ x -y , for all x, y ∈ D.
Then:

• T has a unique fixed point x * in D;

• For any x • ∈ D, d 1 + λ x • -x * d 1 -λ , where d = T(x • ) -x • .
Most of the problems encountered in scientific computing are certainly not of this form. In this context, Newton-like validation methods (also called Krawczyk methods) reuse the ideas behind Newton-Raphson iterations in numerical analysis to transform the initial problem into an instance on which Theorem 2.1 applies. Assume that the original equation is of the form F(x) = 0, with F : X → Y a Fréchet differentiable map between Banach spaces X and Y . The exact solution x * is supposed to be simple root of F, so that the Fréchet derivative DF(x) is invertible in a neighborhood of x * . We are given an approximation x • of x * .

1. Compute a bounded injective linear operator A approximating DF(x • ) -1 , the inverse of the Fréchet derivative of F at x • . Then F(x) = 0 is equivalent to the fixed-point equation associated to the Newton-like operator T:

T(x) = x, where T(x) = x -A(F(x)).
2. Compute rigorously a closed subset D ⊆ X containing x • (e.g. a ball centered at x • with suitable radius r) and a contraction ratio λ < 1 satisfying the assumptions of Theorem 2.1.

Return the error bound

ε = A(F(x • )) /(1 -λ) as in Theorem 2.1. If F has several roots in X,
also ensure that the unique fixed-point of T inside D is necessarily x * .

Remark 2.2. The IVP (1) (or the equivalent integral equation (2)) we consider in this article are linear, which somewhat simplifies the scheme above. The Newton-like operator T is indeed affine, and hence it is locally contracting if and only if it is globally contracting (i.e. D = X), and λ is given by the operator norm of its linear part. Moreover, the solution being unique by the Picard-Lindelöf theorem, no additional verification is required in step 3 above.

Reformulation of Linear ODE into Integral Equation

Several reasons may lead us to prefer integral equations rather than differential ones in numerical analysis. From a theoretical point of view, the differentiation operator ∂ can rarely be seen as a bounded linear endomorphism in one function space, since it typically reduces the regularity. Manipulating differential equations thus requires to consider a separate function space for each iterated application of ∂ (which is for example the approach of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] using ultraspherical polynomials). On the contrary, the operation of integration increases the regularity, which allows us to work in a single function space.

Moreover, the compactness property, necessary for the projection approach of spectral methods, easily follows (see Section 2.3). On the practical side also, integral equations often lead to better conditioned numerical problems [START_REF] Greengard | Spectral integration and two-point boundary value problems[END_REF].

There are several possibilities to transform a linear ODE into an equivalent integral equation, acting on the same function y or one of its iterated derivatives. In this article, we focus on two such transforms, giving integral equations acting on y and y (r) , respectively. The obvious "duality" between them will play an important role later in Section 3.

Integral Equation on

Highest-Order Derivative. This transformation, considered for example in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], is the simplest one, acting on f = y (r) . We here consider a differential equation L{y} = h with a linear differential operator L given in a "standard form", which we call x -∂ form:

L = ∂ r + a r-1 (x)∂ r-1 + • • • + a 1 (x)∂ + a 0 (x). (3)
Proposition 2.3. The linear initial value problem:

L{y}(x) = h(x), x ∈ [x l , x r ],
y (i) (x 0 ) = v i , 0 i < r, with L in x -∂ form (3),
is equivalent to the integral equation on f = y (r) :

f (x) + K{f }(x) = g(x), where K{f }(x) = x x0 K(x, t)f (t)dt, with K(x, t) = r-1 i=0 a i (x) (x -t) r-1-i (r -1 -i)! , and
g(x) = h(x) - r-1 i=0 a i (x) r-1 j=i (x -x 0 ) j-i (j -i)! v j .
Proof. By a repeated use of the integration by part formula, one has for 0 i < r:

y (i) (x) = x x0 t1 x0 • • • tr-1-i x0 f (t r-i )dt r-1 . . . dt 1 + r-1 j=i (x -x 0 ) j-i (j -i)! v j = x x0 (x -t) r-1-i (r -1 -i)! f (t)dt + r-1 j=i (x -x 0 ) j-i (j -i)! v j .
Doing the substitutions in L{y} = h, regrouping the integral terms and moving the initial condition terms to the right, one obtains the expected kernel K(x, t) and the right-hand side g(x).

2.2.2. Integral Equation on the Same Function. Another possible reformulation, used for example in [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF], is to define an integral equation directly on f = y. To do so, we need to have L in a ∂ -x form, with the polynomial (or analytic) coefficients a i (x) located to the right of the ∂ symbols:

L = (-∂) r + (-∂) r-1 a r-1 (x) + • • • -∂a 1 (x) + a 0 (x). (4)
We systematically use -∂ instead of simply ∂ in ∂ -x forms to be consistent with the adjunction L → L * that will play a crucial role in Section 3.3. While a x -∂ form linear ODE relates together the iterated derivatives y (i) where y (i+1) = ∂{y (i) }, it is more natural with a ∂ -x form to consider functions y [i] for 0 i r, with y [0] = y, y [i+1] = a r-i-1 y -∂{y [i] }, and then equating y [r] (x) = h(x). By doing so, we have y [i] = L [i] {y} where:

L [i] = (-∂) i + (-∂) i-1 a r-1 (x) + • • • -∂a r-i+1 (x) + a r-i (x),
so that the "natural" ∂ -x initial conditions at x 0 are y

[i] (x 0 ) = L [i] {y}(x 0 ) = w i for 0 i < r.
Remark 2.4. It is clear that, as long as one can differentiate the a i (x) symbolically sufficiently many times, one easily converts between x -∂ and ∂ -x operators L and initial conditions. Note however that for large r, this may result in a differential operator with enormous coefficients, potentially leading to numerical stability issues.

The integral equation on f = y follows by merely integrating the linear ODE r times.

Proposition 2.5. The linear initial value problem:

L{y}(x) = h(x), x ∈ [x l , x r ], L [i] {y}(x 0 ) = w i , 0 i < r, with L in ∂ -x form (4),
is equivalent to the integral equation on y:

y(x) + K{y}(x) = g(x), where K{y}(x) = x x0 K(x, t)y(t)dt, with K(x, t) = - r-1 i=0 (t -x) r-1-i (r -1 -i)! a i (t), and
g(x) = - x x0 (t -x) r-1 (r -1)! h(t)dt + r-1 i=0 w i (x 0 -x) i i! .
Proof. We first prove by induction on 0 k r: . . .

t k-1 x0 y [k] (t k )dt k . . . dt 1 = (-1) k y - x x0 r-1 i=r-k (t -x) r-1-i (r -1 -i)! a i (t)y(t)dt - k-1 i=0 (x 0 -x) i i! w i . ( 5
)
This is trivially true for k = 0. If it holds for k < r, since y

[k+1] = a r-k-1 y -∂{y [k] }, then: x x0 . . . t k x0 y [k+1] (t k+1 )dt k+1 . . . dt 1 = x x0 . . . t k x0 a r-k-1 (t)y(t)dt k+1 . . . dt 1 - x x0
. . .

t k-1 x0 (y [k] (t k ) -w k )dt k . . . dt 1 = x x0 (x -t) k k! a r-k-1 (t)y(t)dt - x x0
. . .

t k-1 x0 y [k] (t k )dt k . . . dt 1 + (x -x 0 ) k k! w k = (-1) k+1 y - x x0 r-1 i=r-k-1 (t -x) r-1-i (r -1 -i)! a i (t)y(t)dt - k i=0 (x 0 -x) i i! w i .
This completes the proof of [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF]. Taking this equality with k = r allows us to rewrite:

(-1) r x x0
. . .

t k-1 x0 y [r] (t k )dt k . . . dt 1 = (-1) r x x0
. . .

t k-1 x0 h(t k )dt k . . . dt 1 ,
into the desired integral equation with expected K(x, t) and g(x).

Chebyshev Approximation Theory and Spectral Methods

When dealing with analytic functions over a compact real interval [x l , x r ], Chebyshev approximation theory turns out to be an efficient generic tool to provide global, smooth and remarkably accurate approximants. Here we simply recall some fundamentals of this beautiful and rich theory, to which numerous introductory books are dedicated [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF][START_REF] Fox | Chebyshev Polynomials in Numerical Analysis[END_REF][START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF].

The Chebyshev polynomials form a graded family of univariate polynomials, commonly defined over [-1, 1] by the three-term recurrence:

T 0 (x) = 1, T 1 (x) = x, T n+2 (x) = 2xT n+1 (x) -T n (x) for n 0.
The mantra from a complexity point of view is that manipulating polynomials in the Chebyshev basis is asymptotically equivalent to the same operations in the monomial basis. Indeed, multiplication and differentiation/integration are given by the formulas:

T n (x)T m (x) = 1 2 (T n+m (x) + T n-m (x)), for 0 m n, (6)
T n (x) = 1 2

T n+1 (x) n + 1 - T n-1 (x) n -1 , for n 2. ( 7
)
Remark 2.6. Multiplication of two degree n polynomials in the Chebyshev basis can be computed in O(n log n) floating-point arithmetic operations thanks to fast algorithms for Discrete Cosine Transform (DCT), which is the analog to the Fast Fourier Transform (FFT). However, designing a numerically stable algorithm for fast multiplication of polynomials with interval coefficients (in the monomial or Chebyshev basis) is still an open and debated topic3 . This is why in Section 3 we will assume a quadratic complexity when interval arithmetics is needed.

From an analysis point of view, the relation T n (cos ϑ) = cos(nϑ) strongly connects Chebyshev polynomials with Fourier approximation theory. In particular, |T n (x)| 1 for -1 x 1 and they form an orthogonal family w.r.t. the scalar product:

f, g = 1 -1 f (x)g(x) √ 1 -x 2 dx = π 0 f (cos ϑ)g(cos ϑ)dϑ,
which motivates the definition of Chebyshev coefficients for a continuous function f over [-1, 1] by orthogonal projection:

[f ] n = f, T n T n , T n =      1 π f, 1 for n = 0, 2 π f, T n for n 1.
Remark 2.7. By an affine change of variable, the above definitions are transposed from [-1, 1] to any compact interval. We assume in this article that the Chebyshev polynomials T n are defined over [x l , x r ].

Chebyshev approximations are defined by truncated Chebyshev series using the truncation (or projection) operator:

п N : C 0 → R N [x], f → п N {f } = N n=0 [f ] n T n .
The key point of Chebyshev approximation theory is that the more regular f is, the faster п N {f } converges to f as N → ∞. In the analytic case we consider -and contrary to Taylor approximations -convergence on [x l , x r ] always holds, at exponential rate with ratio ρ determined by the location of singularities in the complex plane (Lemma 3.15 in Section 3 gives a quantitative estimate).

When working with infinite sequences of Chebyshev coefficients, a relevant subspace of C 0 is the space Ч 1 of absolutely summable Chebyshev series [8, Sec. 2.2], which is the analog of the Wiener algebra A(T) for Fourier series.

Definition 2.8. The Ч 1 space is the space of continuous functions f over [x l , x r ] for which the quantity,

f Ч 1 = ∞ n=0 |[f ] n |, is finite. This makes (Ч 1 , • Ч 1 ) a Banach space, with the property f Ч 1 f ∞ .
Obviously, Chebyshev coefficients [f ] n cannot be directly computed most of the time. Turning back to the integral equation ( 2), a spectral method approach [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF][START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF] is to cast this problem into an infinite-dimensional linear system in the Chebyshev space of coefficients. More specifically, the integral operator K of polynomial kernel K(x, t) becomes a bounded endomorphism of Ч 1 , acting on the vector of Chebyshev coefficients of a function f as an infinite-dimensional matrix (depicted in Figure 1a). The multiplication [START_REF] Bournez | On the complexity of solving initial value problems[END_REF] and integration (7) formulas imply that K has a sparse structure called almostbanded (following the terminology of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF]), that is, nonzero coefficients are located in a finite number of rows and diagonals. Moreover, K is a compact endomorphism of Ч 1 , where coefficients in the matrix tend to zero. Therefore, spectral methods relax this infinite-dimensional problem by projecting it onto a finite-dimensional subspace, using the truncated integral operator K [N ] = п N K п N (see Figure 1b). The resulting almost-banded linear system of dimension N + 1 is eventually solved in linear time w.r.t. the approximation degree N using the QR-like algorithms of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF].

Proposition 2.9 (see [START_REF] Olver | A fast and well-conditioned spectral method[END_REF]Sec. 5]). The solution of the integral equation (2) can be approximated by a degree N polynomial in the Chebyshev basis using a spectral method requiring O(N (r + s) 2 ) arithmetic operations. More precisely:

• A first algorithm AlmostBandedQRFactor computes a (structured) QR factorization of I N +1 + K [N ] in O(N (r + s) 2 ) arithmetic operations. • Given this QR factorization, a second algorithm AlmostBandedBackSubs performs a backsubstitution to solve:

(I N +1 + K [N ] ){y • } = п N {g}, of unknown y • = N n=0 c n T n , in O(N (r + s)) arithmetic operations.

A Brief Summary of the Newton-Galerkin Validation Method

Problem 1.1 typically falls within the class of function space problems that can be solved using the fixed-point based a posteriori validation scheme presented in Section 2.1. One needs to compute a rigorous error bound between an approximation f • and the exact solution f * of the affine equation:

F{f } = 0, where F{f } = f + K{f } -g : Ч 1 → Ч 1 .
We recall the main ideas of Newton-Galerkin validation method, in the way it is presented in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], with complexity analysis results and discussion about its shortcomings.

Newton-Galerkin Fixed-Point

Operator and Validation Algorithm. The principle of the Newton-Galerkin validation algorithm is to construct the linear operator A ≈ (I + K) -1 (required by step 1 of Newton-like validation) using finite dimensional projections. Indeed, K [N ] converges to K, so the algorithm chooses a suitable truncation order N G (which is totally independent of the approximation f • to be validated) and computes a numerical inverse A of the (N G + 1)-square matrix

I N G +1 + K [N G ] .
A is defined by extending the matrix A over Ч 1 with identity (see Figure 2a).

(a) A (b) K -K [N G ] Figure 2. Bounding the truncation error A(K -K [N G ] ) Ч 1 in Newton-Galerkin method.
After that, in step 2, a contraction ratio λ for the resulting Newton-Galerkin fixed-point operator T is automatically computed by rigorously bounding the operator norm (w.r.t. the Ч 1 -norm) of the linear part, simply denoted by DT:

DT Ч 1 = I -A(I + K) Ч 1 I -A(I + K [N G ] ) Ч 1 + A(K [N G ] -K) Ч 1 . (8)
• The first quantity

I-A(I+K [N G ]
) Ч 1 is due to the fact that A is computed as a numerical inverse (either dense, or even approximated by an almost-banded matrix, see [8, sec. 4.2]). Bounding this quantity only requires straightforward operations on finite-dimensional matrices.

• The second quantity A(K [N G ] -K) Ч 1
, where the two matrices of this product are depicted in Figure 2, is a consequence of the fact that K [N G ] is only an approximation of K. The difficulty, addressed in details in [8, Sec. 5.1], consists in rigorously computing a tight bound for the resulting infinite-dimensional tail, in reasonable complexity. Indeed, gross overestimations would force us to use a larger truncation index N G to ensure a contraction ratio λ < 1, thus strongly impacting the overall efficiency.

Finally, the defect

d = T{f • } -f • Ч 1 = A{f • + K{f • } -g} Ч 1 (step 3
) is rigorously computed in a straightforward manner, yielding the final bound d/(1 -λ).

Theorem 2.10 (see [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]Prop. 5.3]). Let f • be an approximate solution for the integral equation (2). Then, given a validation truncation index N G , algorithms in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], if they do not fail, realize the following4 :

• A Newton-Galerkin fixed-point validation operator T with truncation index N G is computed and bounded in O N 2 G (r + s) arithmetic operations,
where r is the order of L and s = max i deg a i in Problem 1.1. • Having this Newton-Galerkin operator and if the rigorously computed contraction ratio λ < 1, then a candidate approximation f • (with right-hand side g) is validated in

O N G (n + m + r + s) + n(r + s) 2 arithmetic operations,
where n = deg f • and m = deg g.

Obviously, the overall efficiency of the procedure depends on how large N G must be to ensure a contracting operator T. This crucial question is discussed in the next section.

Limitations of the Newton-Galerkin Validation Method.

As pointed out in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], the Newton-Galerkin validation method, although widely used in the mathematical community of computerassisted proofs, suffers from several limitations, especially when targeting certified implementation in a proof assistant.

The first shortcoming is the difficulty to process "harder" instances, e.g. LODEs with highly oscillating or rapidly increasing solutions, due to the value of the truncation index N G . Indeed, K [N G ] converges to K in O(1/N G ) only, while the Ч 1 -norm of (I + K) -1 (and hence the approximation A) is exponentially bounded by the magnitude a i Ч 1 of the input LODE coefficients. This explains why the minimum N G to ensure a contracting operator T can grow exponentially fast with the a i Ч 1 [8, Sec. 5.2]. Although other norms for the space of Chebyshev coefficients have been investigated, e.g. in [START_REF] Lessard | Rigorous numerics for nonlinear differential equations using Chebyshev series[END_REF][START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach[END_REF], it seems that none of them avoid this exponential bound.

Another limitation of Newton-Galerkin validation methods is that implementations can hardly be generic w.r.t. the approximation setting. Indeed, changing the approximation basis from (T n ) to another (ϕ n ) (e.g., Fourier, Legendre) or the norm implies substantial modifications in the technical formulas used to bound the operator norm of T. Moreover, not all norms are suitable when working with coefficient spaces. This is for example why we use

• Ч 1 to overapproximate • ∞ here.
Finally, the need to consider infinite-dimensional matrix representations of the involved operators, as well as the rather technical formulas that are used to bound them, makes Newton-Galerkin not the best candidate for a certified implementation in a proof assistant. The Newton-Picard validation algorithm presented in the next section not only offers a far more satisfying solution to the abovementioned issues of complexity and genericness, it is also in essence more "algebraic" and adapted to formal proof.

Newton-Picard Validation Algorithm

The conclusion to draw from Newton-Galerkin's pessimistic complexity estimates is that (I + K) -1 is poorly approximated by finite-dimensional truncations (I + K [N G ] ) -1 . By contrast, numerical experiments (see for example Figure 3c related to the Airy function example in Section 4.2) show that (I + K) -1 is "asymptotically" almost-banded, meaning that it is very well approximated by almostbanded operators of moderate band width. This phenomenon is investigated in Section 3.1 from the point of view of integral operators. More specifically, we prove that (I + K) -1 = I + R, where R is an integral operator of kernel R(x, t) called resolvent kernel. This elegant framework is the cornerstone of the Newton-Picard validation algorithm presented in Section 3.2.

Contrary to the validation method [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF] which performs explicit Picard iterations to construct a contracting operator, our novel approach relies on the Chebyshev spectral method to compute a polynomial approximation of the resolvent kernel R(x, t) (Section 3.3). The resulting complexity is carefully investigated in Section 3.4 and compared with Newton-Galerkin.

Kernel Composition, Picard Iterations and the Resolvent Kernel

The essence of Newton-Picard validation is to represent an integral operator K not as an infinitedimensional matrix in some coefficient space, but using its kernel K(x, t) defined over [x l , x r |x 0 ]:

K{f }(x) = x x0 K(x, t)f (t)dt, f ∈ C 0 , [x l , x r |x 0 ] = {(x, t) ∈ R 2 s.t. x l x t x 0 or x 0 t x x r }. (9)
The following lemma shows the bijectivity of this correspondence.

Lemma 3.1. Let K be an integral operator with a continuous kernel K(x, t) as in [START_REF] Bréhard | A certificate-based approach to formally verified approximations[END_REF].

Then K{f } = 0 for all f ∈ C 0 if and only if K(x, t) = 0 for all (x, t) ∈ [x l , x r |x 0 ].
Proof. The "if" part of the proof is straightforward. For the "only-if" part, let x ∈ [x l , x r ] and suppose, without loss of generality, that x x 0 . Define f ∈ C 0 by f (t) = K(x, t) for t ∈ [x 0 , x] and extend it constantwise over [x l , x 0 ] and [x, x r ]. Then:

K{f }(x) = x x0 K(x, t)f (t)dt = x x0 K(x, t) 2 dt = 0, which implies K(x, t) = 0 for all t ∈ [x 0 , x]. Therefore, K(x, t) = 0 for all (x, t) ∈ [x l , x r |x 0 ].
Another nice property of the class of integral operators is its closedness under composition, where the resulting kernel is given by a "convolution-like" operation, denoted by * throughout this article. Lemma 3.2. Let K and L be integral operators with respective continuous kernels K(x, t) and L(x, t). Then the composition K L is again an integral operator, with kernel K * L given by the associative law:

(K * L)(x, t) = x t K(x, s)L(s, t)ds, (x, t) ∈ [x l , x r |x 0 ].
In particular, if K and L are polynomials, then so is K * L, with:

deg K * L deg K + deg L + 1.
Proof. The kernel composition formula simply follows from an application of Fubini's theorem:

K{L{f }} = x x0 K(x, s) s x0 L(s, t)f (t)dtds = x x0 x t K(x, s)L(s, t)f (t)dsdt = x x0 (K * L)(x, t)f (t)dt.
Moreover, since the composition • of operators is associative and by injectivity given by Lemma 3.1, we deduce that the operation * of kernel composition is associative too.

Historically, the expression of (I + K) -1 involving an integral operator arose from the Picard iterations [START_REF] Lalescu | Introduction à la Théorie des Équations Intégrales (Introduction to the Theory of Integral Equations[END_REF][START_REF] Rall | Resolvent kernels of Green's function kernels and other finite-rank modifications of Fredholm and Volterra kernels[END_REF], which build a sequence (f n ) of functions converging to the solution f * of (2):

f 0 = g, f n+1 = g -K{f n }, n 0.
This indeed motivates to write (I + K) -1 (at least symbolically) as the Neumann series:

(I + K) -1 = I -K + K 2 -• • • + (-1) n K n + . . . (10)
This matches with the observation that (I + K) -1 is "asymptotically" almost-banded, since K n is almost-banded, of band width n times larger than K. The classic Theorem 3.3 below, for which we provide a proof for the sake of completeness, makes it rigorous by establishing the convergence of this series. This gives rise to the resolvent kernel R(x, t) that plays a central role in Newton-Picard validation method.

Theorem 3.3. Let K of continuous kernel K(x, t) act as an endomorphism of C 0 . Then (I + K) -1 = I + R, where the resolvent kernel R(x, t) is defined using the iterated kernels

K * n : R = ∞ n=1 (-1) n K * n , where K * 1 = K, K * n+1 = K * K * n , n 1. ( 11
)
We have the following estimates:

|R(x, t)| M e M |x-t| , for (x, t) ∈ [x l , x r |x 0 ], (I + K) -1 ∞ = I + R ∞ e M T ,
where M = sup

(x,t)∈[x l ,xr|x0] |K(x, t)|.
Proof. First, we prove by induction on n 1 that:

|K * n (x, t)| M n |x -t| n-1 (n -1)! for (x, t) ∈ [x l , x r |x 0 ].
This clearly holds for n = 1. Suppose that this is true for n. Take w.l.o.g. x 0 t x x r :

|K * n+1 (x, t)| x t |K(x, s)||K * n (s, t)|ds x t M n+1 (s -t) n-1 (n -1)! ds = M n+1 (x -t) n n! .
It follows therefrom that the series (11

) defining R(x, t) converges in norm, that K n ∞ (M T ) n
n! , and that the Neumann series [START_REF] Brisebarre | Chebyshev interpolation polynomial-based tools for rigorous computing[END_REF] makes sense, with the expected estimates for |R(x, t)| and I+R ∞ .

The notion of resolvent kernel is closely related to that of Green's function associated to a linear differential operator [START_REF] Rall | Resolvent kernels of Green's function kernels and other finite-rank modifications of Fredholm and Volterra kernels[END_REF]. The validation methods proposed in [START_REF] Kedem | A posteriori error bounds for two-point boundary value problems[END_REF] and [START_REF] Plum | Computer-assisted existence proofs for two-point boundary value problems[END_REF] rigorously bound the Green's function to provide a bound on the exact inverse operator associated to L. Our validation algorithm, presented in next section, proceeds differently: it approximates the resolvent kernel and rigorously bounds the operator norm of the resulting Newton-like operator T. This allows for sharper bounds and a more convenient complexity analysis.

The Validation Algorithm

Algorithm NewtonPicardValid below implements the three steps of the Newton-like a posteriori validation scheme presented in Section 2.1, by approximating the inverse derivative (I + K) -1 as I + R • , where R • is an integral operator of polynomial kernel R • (x, t) approximating the analytic resolvent kernel R(x, t). Similarly to the truncation index N G in Newton-Galerkin validation procedure, a parameter N R (the approximation degree for R • (x, t)) is passed to the validation algorithm. In general, N R is at least r + s, where r is the order of L and s = max 0 i<r deg a i , so that we make the assumption r + s = O(N R ) in the complexity estimates.

1. First, the subroutine ResolventKernel (postponed in Section 3.3) computes degree N R polynomials α • i , β • i and sets:

R • (x, t) = r-1 i=0 α • i (x)β • i (t) ≈ R(x, t), deg α • i , deg β • i N R . (12)
This defines the Newton-Picard fixed-point operator

T{f } = f -(I + R • ){f + K{f } -g}.
2. The linear part of affine operator T:

E = I -(I + R • )(I + K) = -K -R • -R • K, is an integral operator whose polynomial kernel E = -K -R • -R • * K is explicitly computed as i,j ε ij T i (x)T j (t)
. This yields the bound λ:

E ∞ = sup x∈I sgn(x -x 0 ) x x0 |E(x, t)|dt T i,j |ε ij | = λ.
3. If λ 1, then go back to step 1 using a larger approximation degree N R . Otherwise, compute and bound the defect as d, using, for instance, the Ч 1 norm (or any other more refined, yet rigorous overestimation of the uniform norm):

T{f • } -f • ∞ (I + R • ){f • + K{f • } -g} Ч 1 = d,
and return d/(1 -λ).

For the sake of efficiency, kernels are represented using low-rank decompositions.

Definition 3.4. A bivariate kernel K(x, t) is said to be of rank k if there are functions α 0 (x), . . . , α k-1 (x) and β 0 (t), . . . , β k-1 (t) over [x l , x r ] such that:

K(x, t) = k-1 i=0 α i (x)β i (t), (x, t) ∈ [x l , x r |x 0 ].
The kernel K(x, t) in Propositions 2.3 and 2.5 is clearly of rank r. Also, we shall see in the next section that the resolvent kernel R(x, t) is of rank r too (which is not obvious from its definition using iterated kernels in Theorem 3.3). This explains why we compute a R • (x, t) of rank r to approximate R(x, t) in [START_REF] Epstein | Ultra-arithmetic. I. Function data types[END_REF]. Therefore, elementary operations on low rank kernels (evaluation on a polynomial and composition with * ) are necessary to implement these three steps in Algorithm NewtonPicardValid. Their complexity is investigated in the subsequent Lemmas 3.7 and 3.8.

Algorithm 1 NewtonPicardValid(K, g, f • , N R )
Input: The polynomial kernel K(x, t) and polynomial right-hand side g(x) coming from L{y} = h using Proposition 2.3 or 2.5, a polynomial approximation f • (x) and a validation degree N R .

Output: A rigorous error bound ε f • -f * ∞ .
Step 1: approximating the resolvent kernel R(x, t) using floating-point arithmetics

1: Write K(x, t) as            r-1 i=0 κi(x)Ti(t) if L is in x -∂ form or - r-1 i=0 Ti(x)κi(t) if L is in ∂ -x form 2: R • (x, t) = r-1 i=0 α • i (x)β • i (t) ← ResolventKernel(K, N R )
Step 2: Bounding the contraction ratio of T using interval arithmetics

3: M(x, t) ← KernelComp(R • , K) 4: E(x, t) = i,j ε ij T i (x)T j (t) ← -R • (x, t) -K(x, t) -M(x, t) 5: λ ← T i,j |ε ij |
Step 3: Compute the final bound using interval arithmetics 6: if λ < 1 then

7: δ(x) ← f • (x) + K{f • }(x) -g(x) 8: δ(x) ← δ(x) + R • {δ}(x) 9: d ← δ Ч 1 10: return d 1-λ 11: else 12:
return "FAIL" 13: end if Theorem 3.5. Let f • be a degree n polynomial approximation of the solution f * of the integral equation of Proposition 2.3 or 2.5, and N R the degree used for R • (x, t). If NewtonPicardValid(K, g, f • , N R ) does not fail, then it outputs a rigorous error bound for f • -f * . More precisely:

•

Input: Polynomial kernels K(x, t) = k-1 i=0 α i (x)β i (t) and L(x, t) = l-1 j=0 γ j (x)δ j (t). Output: Polynomial kernel M(x, t) = k+l-1 i=0 µ i (x)ν i (t) equal to K * L.
1: for i = 0 to k -1 and j = 0 to l -1 do ζ ij (s) ← β i (s)γ j (s)ds end for 2: for i = 0 to k -1 do 3:

µ i (x) ← α i (x) 4: ν i (t) ← - l-1 j=0 ζ ij (t)δ j (t)
5: end for 6: for j = 0 to l -1 do 7:

µ k+j (x) ← k-1 i=0 α i (x)ζ ij (x) 8:
ν k+j (t) ← δ j (t) 9: end for

10: return M(x, t) ← k+l-1 i=0 µ i (x)ν i (t) Lemma 3.8. Let K(x, t) = k-1 i=0 α i (x)β i (t) and L(x, t) = l-1 j=0
γ j (x)δ j (t) be polynomial kernels of respective ranks k and l, and respective total degrees n and m. Then KernelComp(K, L) computes a rank k + l decomposition of K * L with degree at most n + m + 1, in O(kl(n + m) 2 ) arithmetic operations.

Proof. Algorithm KernelComp first computes primitives in line 1:

ζ ij (s) = β i (s)γ j (s)ds, 0 i < k, 0 j < l,
in O(klnm) arithmetic operations, giving polynomials of degree at most n + m + 1.

The computation of M(x, t) then follows:

(K * L)(x, t) = x t K(x, s)L(s, t)ds = k-1 i=0 l-1 j=0 α i (x)(ζ ij (x) -ζ ij (t))δ j (t) = l-1 j=0 k-1 i=0 α i (x)ζ ij (x) µ k+j (x) δ j (t) ν k+j (t) - k-1 i=0 α i (x) µi(x)   l-1 j=0 ζ ij (t)δ j (t)   νi(t)
.

More precisely, ν i (t) (line 4) requires O(l(n + m)m) arithmetic operations, and similarly µ k+j (x) (line 7) requires O(kn(n + m)) ones. This concludes the proof of the claimed complexity.

Proof of Theorem 3.5. First, for step

d δ Ч 1 δ ∞ computed line 9 requires deg δ(x) = O(2N R + r + s + n + m) arithmetic operations. Thus, O((N R + n + m)N R r) is a compact majorant for the complexity, assuming r + s = O(N R ).
Remark 3.9. The total complexity can be made quasi-linear w.r.t. N R by using a fast DCT-based algorithm for polynomial multiplication in the kernel composition (line 3), and bounding the lowrank kernel E(x, t) (lines 4-5) using an adequate orthogonalization process (to reduce overestimation). However, as explained in Remark 2.6, we choose to keep the quadratic estimate in this article.

Remark 3.10. The use of Chebyshev approximations for R • (x, t) is not crucial to the principle of the validation algorithm. The reason of this choice is that they are asymptotically excellent candidates and allow for an elegant complexity analysis in Section 3.4. However, for intermediate precision, other approximants (e.g., splines) may perform better.

Resolvent Kernel Approximation

Algorithm NewtonPicardValid above requires a polynomial approximation R • (x, t) of R(x, t). We first show in Section 3.3.1 that R admits an analytic rank r expression, involving the solutions of the homogeneous linear ODE L{y} = 0, together with those of the adjoint equation L * {y} = 0. The adjunction L → L * , defined as:

x * = x, ∂ * = -∂, (L 1 L 2 ) * = L * 2 L * 1 , (L 1 + L 2 ) * = L * 1 + L * 2
, plays a crucial role to this purpose. Clearly, the x -∂ form differential operator in (3) and the one in ∂ -x form in (4) are mutually adjoint. For the sake of clarity, we opt for the convention to refer to the former as L, and the latter as L * . We also call K and K(x, t) the integral operator and kernel associated to L by Proposition 2.3, and similarly K * and K * (x, t) for L * in Proposition 2.5. Note the following adjunction relation on integral operators5 :

K * (x, t) = -K(t, x), (x, t) ∈ [x l , x r ] 2 .
The formulas given in Proposition 3.12 below for the resolvent kernels R and R * associated to K and K * also satisfy R * (x, t) = -R(t, x).

Based on this, Section 3.3.2 presents the approximation routine ResolventKernel which numerically solves the homogeneous equations using the Chebyshev spectral method presented in Section 2.3. This completes the algorithmic description of the Newton-Picard validation procedure.

Y (x) =        y(x) y (x)
. . . . . .

y (r-1) (x)        , A(x) =        0 1 0 • • • 0 0 0 1 • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • 1 -a 0 (x) -a 1 (x) -a 2 (x) • • • -a r-1 (x)        , H(x) =        0 . . . . . . 0 h(x)        .
It is well known that the solution with initial conditions v = (v 0 , . . . , v r-1 ) T can be expressed using an integral operator with kernel Φ(x)Φ(t) -1 :

Y (x) = Φ(x) v + x x0 Φ(t) -1 H(t)dt ,
where Φ(x) is the fundamental matrix, satisfying the matrix ODE Φ (x) = A(x)Φ(x) with Φ(x 0 ) = I r . In other words, Φ(x) contains the canonical basis of solutions for the homogeneous linear ODE L{y} = 0 and their derivatives up to order r -1. In particular, its columns are linearly independent for all x. Hence, Φ(x) is invertible and we have:

(Φ(x) -1 ) = -Φ(x) -1 Φ (x)Φ(x) -1 = -Φ(x) -1 A(x)Φ(x)Φ(x) -1 = -Φ(x) -1 A(x),
which means that the inverse transpose x → Φ(x) -T satisfies the matrix ODE M (x) = -A(x) T M (x).

On the other hand, the scalar linear ODE L * {y} = h with L * given in ∂ -x form as in ( 4) is equivalent to the vectorial ODE Y (x) = B(x)Y (x) + H(x), where:

Y (x) =        y(x)
y [1] (x) . . . . . .

y [r-1] (x)        , B(x) =        a r-1 (x) -1 0 . . . 0 a r-2 (x) 0 -1 . . . 0 . . . . . . . . . . . . . . . a 1 (x) 0 0 . . . -1 a 0 (x) 0 0 . . . 0        , H(x) =        0 . . . . . . 0 -h(x)        .
Similarly, the fundamental matrix Ψ(x) is defined as the solution of Ψ (x) = B(x)Ψ(x) with Ψ(x 0 ) = I r . A key observation is that:

J r B(x)J r = -A(x) T , where J r =    0 . . . 1 . . . . . . . . . 1 . . . 0    , (13)
leading to the following crucial orthogonality relations.

Lemma 3.11. The canonical basis ϕ 0 , . . . , ϕ r-1 for the solution space of L:

L{ϕ i } = 0, ϕ (j) i (x 0 ) = 1 if i = j 0 otherwise,
and the dual basis ψ 0 , . . . , ψ r-1 for the solution space of the adjoint L * :

L * {ψ i } = 0, ψ [j] i (x 0 ) = 1 if i = j 0 otherwise,
satisfy the orthogonality relations for 0 j, k < r:

r-1 i=0 ϕ (j) i (x) ψ [k] r-1-i (x) = 1 if k + j = r -1, 0 otherwise, for all x ∈ [x l , x r ]. (14)
Proof. The matrices Φ(x) = (ϕ (i) j (x)) 0 i,j<r and Ψ(x) = (ψ

[i] j (x)) 0 i,j<r satisfy Φ(x 0 ) = Ψ(x 0 ) = I r , and are solutions to the matrix linear ODEs Φ (x) = A(x)Φ(x) and Ψ (x) = B(x)Ψ(x), respectively. From identity (13), we have that:

(J r Ψ(x)J r ) = J r B(x)Ψ(x)J r = -A(x) T (J r Ψ(x)J r ).
Hence, x → J r Ψ(x)J r satisfies the same ODE M (x) = -A(x) T M (x) as x → Φ(x) -T , and both are equal to I r at x 0 . It follows therefrom that Φ(x) -T = J r Ψ(x)J r , and in particular:

Φ(x)J r Ψ(x) T J r = I r for all x ∈ [x l , x r ],
which unfolds as the orthogonality relations [START_REF] Fox | Chebyshev Polynomials in Numerical Analysis[END_REF].

Proposition 3.12. The resolvent kernels R and R * are given by:

R(x, t) = r-1 i=0 ϕ (r) i (x)ψ r-1-i (t), (15)
R * (x, t) = - r-1 i=0 ψ i (x)ϕ (r) r-1-i (t), (16)
with the ϕ i , ψ i defined in Lemma 3.11.

Proof. • For R, call R (x, t) the right-hand side of ( 15) and R the integral operator defined by this kernel. Since, by Theorem 3.3, I + K is already known to be invertible, of inverse I + R, it is sufficient to prove that

(I + K)(I + R ) = I, that is K * R + R + K = 0. For 0 i < r, L{ϕ i } = 0 implies from Proposition 2.3 by renaming x 0 into t: ϕ (r) i (x) + x t K(x, s)ϕ (r) i (s)ds = - r-1 j=0 a j (x) r-1 k=j (x -t) k-j (k -j)! ϕ (k) i (t).
Multiplying the left-hand side by ψ r-1-i (t) and summing over i gives R (x, t) + (K * R )(x, t). Doing the same operations on the right-hand side and commuting the summation symbols allows us to isolate the term

r-1 i=0 ϕ (k)
i (t)ψ r-1-i (t), which by Lemma 3.11 is nonzero (and equal to 1) if and only if k = r-1. We thus obtain the expected equality:

R (x, t) + (K * R )(x, t) = - r-1 j=0 a j (x) (x -t) r-1-j (r -1 -j)! = -K(x, t).
• For R * , call R (x, t) the right-hand side of [START_REF] Greengard | Spectral integration and two-point boundary value problems[END_REF]. Clearly, R (x, t) = -R(t, x). We have:

(K * * R )(x, t) + R (x, t) + K * (x, t) = x t K * (x, s)R (s, t)ds + R (x, t) + K * (x, t) = x t K(s, x)R(t, s)ds -R(t, x) -K(t, x) = -(R * K)(t, x) -R(t, x) -K(t, x) = 0, since -R K -R -K = I -(I + R)(I + K) = 0.
Hence the corresponding operator R satisfies (I + K * )(I + R ) = I, which implies R = R * and hence R = R * .

Resolvent Kernel Approximation Algorithm.

A major advantage of the formulas in Proposition 3.12 from the computational point of view, is that in both x -∂ and ∂ -x cases, one requires the ψ i , solutions to an ODE in ∂ -x form, and the ϕ (r) i

where the ϕ i are solutions to an ODE in x -∂ form. Hence, no conversion from x -∂ and ∂ -x forms is needed here.

Algorithm ResolventKernel below computes degree N R polynomial approximations of these analytic functions. To do so, the Chebyshev spectral method of Section 2.3 is applied on the following integral equations, equivalent by Propositions 2.3 and 2.5 to the IVPs defining ϕ i and ψ i in Lemma 3.11:

(I + K){ϕ (r) i } = g i , with g i (x) = - i j=0 a j (x) (x -x 0 ) i-j (i -j)! , (I + K * ){ψ i } = g * i , with g * i (x) = (x 0 -x) i i! . ( 17
)
Algorithm 3 ResolventKernel(K, N R )

Input: Kernel K(x, t) obtained from L using Proposition 2.3 or 2.5, and approximation degree N R .

Output:

Polynomial kernel R • (x, t) = r-1 i=0 α • i (x)β • i (t)
approximating the resolvent kernel R(x, t).

1: K * (x, t) ← -K(t, x) 2: if L is in ∂ -x form then 3: K(x, t) ↔ K * (x, t)
Exchange the role of L and its adjoint

L * 4: end if Compute approximations ϕ (r)• i to the ϕ (r) i associated to L{ϕ} = 0 5: M ← I N R +1 + K [N R ] 6: (Q, R) ← AlmostBandedQRFactor(M ) 7: for i = 0 to r -1 do 8: g ii ← -a i 9: g ij ← g i-1,j x-x0 i-j for j = 0 to i -1 10: g i ← g i0 + • • • + g ii 11: ϕ (r)• i ← AlmostBandedBackSubs(Q, R, g i ) 12: end for Compute approximations ψ • i to the ψ i associated to L * {ψ} = 0 13: M * ← I N R +1 + K * [NR] 14: (Q * , R * ) ← AlmostBandedQRFactor(M * ) 15: g * 0 ← 1 16: for i = 0 to r -1 do 17: ψ • i ← AlmostBandedBackSubs(Q * , R * , g * i )
18:

g * i+1 ← g * i x0-x i+1 19: end for Return polynomial kernel R • (x, t) 20: if L is in x -∂ form then 21: (α • i , β • i ) ← (ϕ (r)• i , ψ • r-1-i ) for i = 0 to r -1 22: else if L is in ∂ -x form then 23: (α • i , β • i ) ← (-ψ • i , ϕ (r)• r-1-i ) for i = 0 to r -1 24: end if 25: return R • (x, t) ← r-1 i=0 α • i (x)β • i (t) Proposition 3.13. Algorithm ResolventKernel computes a rank r, degree 2N R polynomial kernel R • (x, t) approximating the resolvent kernel R(x, t) in O(N R (r + s) 2 ) arithmetic operations.
Proof. First, the computation of M (line 5) and M * (line 13) only requires O(N R r(r + s)) arithmetic operations, due to the specific form of K(x, t) built in line 1 of Algorithm NewtonPicardValid. Indeed, for 0 i < r and 0 j N R , the computation of κ i (x)

x x0 T i (t)T j (t)dt or T i (t)
x x0 κ i (t)T j (t)dt involves two multiplications where at most one operand has O(r + s) nonzero coefficients (O(1) for the other one), and one integral where the integrand has at most O(r + s) nonzero coefficients. Hence, for 0 j N R , the nonzero coefficients of K{T j } and K * {T j } (forming the horizontal and diagonal bands of M and M * ) are computed in O(r(r + s)) arithmetic operations each.

Then, according to Proposition 2.9, the computation of the QR factorizations in lines 6 and 14 using Algorithm AlmostBandedQRFactor requires O(N R (r + s) 2 ) arithmetic operations.

After that, for each 0 i < r, the computation of the right-hand side g i (lines 8-10) requires O(r(r + s)) arithmetic operations, while the updating of the right-hand side g * i of the adjoint equation (line 18) needs O(r + s) ones. Finally, solving the resulting almost-banded linear systems using Al-mostBandedBackSubs in lines 11 and 17 with the previously computed QR factorizations requires O(N R (r+s)) arithmetic operations. Iterating over 0 i < r, this yields a complexity of O(N R r(r+s)) for lines 7-12 and 16-19, assuming r + s = O(N R ).

Hence, the overall complexity of ResolventKernel is O(N R (r + s) 2 ).

Remark 3.14. Instead of computing a fundamental matrix for the adjoint equation L * {y} = 0 (lines 16-19 in ResolventKernel), we can also directly approximate the inverse of the fundamental matrix of L{y} = 0 (already computed in lines 7-12). To do so, sample the fundamental matrix at N R + 1 Chebyshev nodes (with complexity O(N R log N R r 2 ) using fast DCT), invert the matrix at each node (with a total of O(N R r ω ) arithmetic operations), and finally re-interpolate in complexity O(N R log N R r 2 ).

Newton-Galerkin vs Newton-Picard Validation Complexity

The overall complexity of NewtonPicardValid depends on how large the degree N R of the approximate resolvent kernel R • (x, t) has to be to ensure a contracting Newton-Picard operator T. This, in turn, is determined by the approximation errors ε

(N R ) i = ϕ (r)• i -ϕ (r) i ∞ and η (N R ) i = ψ • i -ψ i ∞ of the ϕ (r)• i and ψ •
i computed by ResolventKernel. The use of spectral methods is motivated by both efficiency and accuracy, and in practice the resulting error is close to the one from the mere truncation of the exact Chebyshev series. In the analytic framework we consider (where the coefficients of L are polynomials or analytic functions over a complex neighborhood of [x l , x r ], represented using RPAs), Chebyshev series enjoy far better convergence properties than the simple summability condition defining the Ч 1 space. Let E ρ denote the Bernstein ellipse [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Chap. 8] of parameter ρ > 1 rescaled over [x l , x r ], i.e. the ellipse in the complex plane of foci x l and x r and eccentricity 2/(ρ + ρ -1 ). It is the typical domain of convergence of Chebyshev series, which play the same role as the convergence disks do for Taylor series. The following lemma (see e.g. [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Thm. 8.1]) asserts the exponential decay of Chebyshev coefficients for analytic functions.

Lemma 3.15. Let f be an analytic function over E ρ . Then its Chebyshev coefficients [f ] n satisfy:

|[f ] n | f Eρ for n = 0, 2 f Eρ ρ -n for n 1,
where

f Eρ = sup z∈Eρ |f (z)| is the uniform norm over E ρ . In particular, f -п n {f } ∞ f -п n {f } Ч 1 2 f Eρ ρ -n ρ -1 , n 0.
Let us fix ρ > 1 such that all coefficients of L (and consequently the IVP solutions) are analytic over E ρ . Theorem 3.16 provides an estimate for N R under the simplifying assumption that exact truncated Chebyshev series of degree N R for ϕ (r) i and ψ i can be computed. We discuss later why including the approximation error of the Chebyshev spectral method used in ResolventKernel does not significantly affect it. Theorem 3.16. Under the hypothesis that R(x, t) is approximated using exact truncated Chebyshev series for ϕ (r) i and ψ i :

ϕ (r)• i = п N R {ϕ (r) i }, ψ • i = п N R {ψ i }, 0 i < r, (18)
the minimum value for N R to ensure a contracting Newton-Picard fixed-point operator is bounded by the asymptotic estimate:

N R = O     log 1 + r-1 i=0 a i ∞ T r-i (r-i)! + log T r-1 i=0 ϕ (r) i Eρ ψ r-1-i Eρ -log(ρ -1) log ρ     .
Proof. • We first bound the operator norm of the linear part

E = I-(I+R • )(I+K) = (R-R • )(I+K)
of the Newton-Picard fixed-point operator T as:

E ∞ T r-1 i=0 ( ϕ (r) i ∞ η (N R ) r-1-i + ψ r-1-i ∞ ε (N R ) i + ε (N R ) i η (N R ) r-1-i 1 + r-1 i=0 a i ∞ T r-i (r -i)! . (19)
Indeed, the definition of K(x, t) in Proposition 2.3 or 2.5 yields:

K{f } ∞ r-1 i=0 a i ∞ sup x∈[x l ,xr] sgn(x -x 0 ) x x0 |x -t| r-1-i (r -1 -i)! dt f ∞ r-1 i=0 a i ∞ T r-i (r -1)! f ∞ ,
from which we have the expected bound for K ∞ . Moreover, in the x -∂ case:

R • (x, t) -R(x, t) = r-1 i=0 ϕ (r) i (x)(ψ • r-1-i -ψ r-1-i )(t) + (ϕ (r)• i -ϕ (r) i )(x)ψ r-1-i (t) +(ϕ (r)• i -ϕ (r)
i )(x)(ψ • r-1-i -ψ r-1-i )(t) ,
and a similar decomposition holds in the ∂ -x case. This gives the bound:

R • -R ∞ T r-1 i=0 ( ϕ (r) i ∞ η (N R ) r-1-i + ε (N R ) i ψ r-1-i ∞ + ε (N R ) i η (N R )
r-1-i .

• Instantiating ε

(N R ) i = 2 ϕ (r) i Eρ ρ -N R /(ρ -1) and η (N R ) i = 2 ψ i Eρ ρ -N R /(ρ -1) in (19) yields: E 4T ρ -N R ρ -1 + ρ -2N R (ρ -1) 2 r-1 i=0 ϕ (r) i Eρ ψ r-1-i Eρ 1 + r-1 i=0 a i ∞ T r-i (r -i)! .
Taking the logarithm and reordering the terms gives the claimed asymptotic lower bound for N R .

This estimate depends on the magnitude of the transition matrix associated to the linear ODE, which is consistent with the intuition that stiff ODEs with large amplitude solutions are harder to validate. We can however provide a more conservative estimate depending on the initial data only.

Corollary 3.17. Let T ρ = (x r -x l )(ρ+ρ -1 )/2 T be the major axis of E ρ and

M ρ = r-1 i=0 a i Eρ T r-1-i ρ (r-1-i)! .
Then, under the same simplifying assumption (18), the minimum N R is given by:

N R = O T ρ M ρ -log(ρ -1) log ρ .
Proof. Reusing the notations of ( 17), the formulas ϕ (r) i = (I + R){g i } and ψ i = (I + R * ){g * i } extend over the Bernstein ellipse E ρ in the complex plane. An adaptation of the inverse bound estimate of Theorem 3.3 yields I + R Eρ , I + R * Eρ e TρMρ . Also, for 0 i < r:

g i Eρ g * r-1-i Eρ i j=0 a j Eρ T r-1-j ρ (i -j)!(r -1 -i)! M ρ .
Hence, in the estimate of Theorem 3.16:

log T ρ r-1 i=0 ϕ (r) i Eρ ψ r-1-i Eρ log(rT ρ M ρ e 2TρMρ ) = O(T ρ M ρ ),
from which we obtain the claimed asymptotic estimate.

This estimate confirms the significant efficiency gap between the Newton-Picard and the Newton-Galerkin validation methods for moderate to hard instances. Both algorithms run in quadratic time w.r.t. the parameter N G (truncation index in Newton-Galerkin) or N R (resolvent kernel approximation degree in Newton-Picard). However, contrary to the linear estimate for N R given above, N G can grow exponentially fast w.r.t. the magnitude of the input ODE, as explained in Section 2.4.2.

This conclusion must be tempered by noticing that the quantity M ρ characterizes the magnitude of the input ODE over the Bernstein ellipse E ρ . In fact, a i Eρ may be in theory way larger than a i Ч 1 (and consequently than a i ∞ ). Indeed, for n 0,

T n ∞ = T n Ч 1 = 1, while T n Eρ = O(ρ n
). However, most of the time in practice, the a i are polynomials of reasonable degrees or analytic functions with fast decaying Chebyshev coefficients, so that this ratio still makes N R way smaller than N G . If not, then both Newton-Galerkin and Newton-Picard perform poorly for the mere reason that even the numerical spectral method computing a candidate Chebyshev approximation will require a very high degree, since the Chebyshev coefficients of the solution may be large and decrease slowly.

Including the Spectral Method Error. The estimates for N R given in Theorem 3.16 and Corollary 3.17 above rely on the simplifying assumption (18) that R • (x, t) is built from the exact truncated Chebyshev series for the ϕ (r) i and ψ i defining the resolvent kernel R(x, t). In practice, such truncations cannot be directly obtained, and Algorithm ResolventKernel resorts to spectral methods to compute the desired polynomial approximations ϕ (r)• i and ψ • i . Therefore, we prove in Theorem 3.18 below that the same asymptotic estimate for N R holds in this more realistic setting where method errors of spectral methods are taken into account.

Chebyshev spectral methods perform roughly as well as merely truncating the Chebyshev series (see e.g. [START_REF] Olver | A fast and well-conditioned spectral method[END_REF]Thm. 4.5]), with a convergence in O(ρ -N ) with N the approximation degree in the analytic case. However, the proof is asymptotic, and valid as long as ϕ + K [N ] {ϕ} = g has a solution ϕ with deg ϕ N . Obviously, since I + K is itself invertible, so is I + K [N ] for N sufficiently large, e.g. if

(I + K) -1 (K -K [N ] ) Ч 1
1. But this quantity, already encountered in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], precisely leads to an exponential bound for N = N G in the Newton-Galerkin validation method. To avoid this phenomenon for the degree N = N R used in Newton-Picard, we propose a slightly modified procedure ResolventKernel, by replacing the linear system with a least-squares problem to approximate the exact solution ϕ * = (I + K) -1 {g}. It is important to note that the following development is essentially theoretical and aims at justifying the linear asymptotic estimate for N R given in Corollary 3.17. In practice however, ResolventKernel is still implemented as presented in Section 3.3.2.

Let m = N -r -s, and call K [N,m] = п N K п m the (N + 1) × (m + 1) rectangular projection of K. The key observation is that K [N,m] {ϕ} = K [N ] {ϕ} = K{ϕ} for all ϕ with deg ϕ m. To solve the following least-squares problem, we replace • Ч 1 by an Euclidean norm, namely

h 2 = N i=0 a 2 i for h = N i=0 a i T i , keeping in mind the classical norm equivalence h 2 h Ч 1 √ N + 1 h 2 . min ϕ (I + K [N,m] ){ϕ} -g 2 , s.t. deg ϕ m. (20)
Let ϕ be the optimal solution to [START_REF] Kaucher | Validating computation in a function space[END_REF]. By comparing it with suboptimal ϕ = п m {ϕ * }, we get:

(I + K [N,m] ){ ϕ} -g 2 (I + K [N,m] ){ϕ} -g 2 (I + K){ϕ} -g Ч 1 I + K Ч1 ϕ -ϕ * Ч 1 I + K Ч1 2 ϕ * Eρ ρ -N +r+s ρ -1 .
We therefore deduce the bound:

ϕ -ϕ * ∞ (I + K) -1 ∞ (I + K [N,m] ){ ϕ} -g ∞ (I + K) -1 ∞ √ N + 1 (I + K [N,m] ){ ϕ} -g 2 (I + K) -1 ∞ I + K Ч1 2 ϕ * Eρ √ N + 1 ρ -N +r+s ρ -1 . (21)
Theorem 3.18. With the polynomial approximations ϕ (r)• i and ψ • i computed by ResolventKernel using the modified spectral method presented above, the minimum value for N R to ensure a contracting Newton-Picard fixed-point operator is bounded by the asymptotic estimate:

N R = O T ρ M ρ -log(ρ -1) log ρ .
Proof. Combining the method error for the (modified) Chebyshev spectral method [START_REF] Kedem | A posteriori error bounds for two-point boundary value problems[END_REF] with

(I + K) -1 ∞ e M T e MρTρ and I + K Ч 1 1 + 2MρTρ ρ-1 , we get for ϕ = ϕ (r) i or ψ i and N = N R : ϕ • -ϕ ∞ e MρTρ 1 + 2M ρ T ρ ρ -1 2 ϕ Eρ N R + 1 ρ -N R +r+s ρ -1 .
Instantiating the error bounds ε

(N R ) i
and η

(N R ) i
with these values instead of the mere truncation bounds in the proof of Theorem 3.16 yields the same asymptotic estimate for N = N R in Corollary 3.17.

To conclude, we also note that this modified procedure still runs in O(N R (r + s) 2 ) arithmetic operations. Let M = I N +1,m+1 + K [N,m] . The solutions ϕ to the least squares problem [START_REF] Kaucher | Validating computation in a function space[END_REF] are the solutions of:

M T M ϕ = M T g. ( 22
)
Since M is (rectangular) almost-banded, adapting the algorithms of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] leads to a QR decomposition of M in complexity O(N R (r + s) 2 ) and a solution ϕ of ( 22) in complexity O(N R (r + s)).

Remark 3.19. The overall complexity analysis in this article does not include the effect of rounding errors of floating-point arithmetics and the growth of interval width in the validation part. This could be avoided by using rational numbers instead. However, in most examples we tried, including those presented in the following section, rounding errors do not seem to play an important role. Also, the Ч 1 norm used in line 9 of NewtonPicardValid could be replaced by a polynomial extrema isolation technique to get sharper bounds.

Implementation, Examples and Discussion

We briefly present our implementation of the Newton-Picard validation algorithm presented in this paper, and draw up more detailed observations based on practical examples.

Implementation in the ChebValid library

Algorithm NewtonPicardValid presented in this article is implemented in the C library ChebValid6 dedicated to rigorous numerics with RPAs in the Chebyshev basis (a.k.a. Chebyshev models). This library, which relies on MPFR7 for floating-point arithmetics in arbitrary precision and its extension MPFI8 for interval arithmetics, was initially designed to implement the Newton-Galerkin validation method of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. It implements data types as polynomials with double floating-point numbers, arbitrary precision floating-point numbers or interval coefficients, under the name double_chebpoly_t, mpfr_chebpoly_t and mpfi_chebpoly_t. The type chebmodel_t implements Chebyshev models, that is a pair consisting of a mpfi_chebpoly_t and a floating-point bound for the remainder.

For the purpose of this article, we enriched ChebValid with low-rank integral operators, implemented as the type prefix_lrintop_t (where prefix_ = double_chebpoly_, mpfr_chebpoly_ or chebmodel_), together with relevant operations, in particular evaluation on a polynomial or RPA, composition (that is, operation * on the kernels), and rigorous bound. The validation method is implemented in mpfi_chebpoly_newtonpicard/. In Section 4.2 below, we test it on a D-finite example, namely the Airy function of the first kind, and discuss the performances compared to other a-posteriori or self-validating methods. Also, Section 4.3 exemplifies the extension of NewtonPicard-Valid to non-polynomial coefficients, represented using RPAs (i.e., chebmodel_t), and implemented in chebmodel_newtonpicard/.

A D-finite Example: the Airy Function

The Airy function of the first kind Ai is a special function, notably used in different branches of physics (e.g., quantum mechanics and optics). It is a D-finite function, satisfying the following second order linear ODE:

y (x) -xy(x) = 0, (23)
with transcendental initial conditions v 0 = Ai(0) = 3 -2/3 /Γ( 23 ) and v 1 = Ai (0) = -3 -1/3 /Γ( 13 ). While Ai converges exponentially fast as x → +∞, namely Ai(x) ∼ x -1 4 e -2 3 x 3/2 /2 √ π, the Airy function of the second kind Bi, which together with Ai forms a basis of the solution space of [START_REF] Lalescu | Introduction à la Théorie des Équations Intégrales (Introduction to the Theory of Integral Equations[END_REF], is exponentially divergent, with Bi(x) ∼ x -1 4 e 2 3 x 3/2 / √ π as x → +∞ (see Figure 3a). This means that even slight perturbations of the (non-rational) initial conditions v 0 and v 1 result in exponentially growing errors when evaluating Ai(x) using the D-finite equation [START_REF] Lalescu | Introduction à la Théorie des Équations Intégrales (Introduction to the Theory of Integral Equations[END_REF]. This justifies the use of rigorous validation methods.

First, this IVP for Ai is transformed into an equivalent integral equation on f = Ai using Proposition 2.3:

f (x) + x 0 (xt -x 2 )f (t)dt = v 0 x + v 1 x 2 . (24)
A polynomial approximation y • is computed using the Chebyshev spectral method described in Section 2.3. The magnitude of the Chebyshev coefficients of Ai over the symmetric interval [-a, a] = [-15, 15] is depicted in Figure 3b. The observed super-linear convergence in log-scale (that is, in o(ρ -n ) for any ρ > 1) is typical for entire functions. According to the plot, a Chebyshev polynomial approximation y • with double-precision accuracy (≈ 10 -16 ) needs to be of degree roughly 90. However, such a heuristic truncation, as done for example in the Chebfun library, is not a rigorous proof. As clearly evidenced by Figure 3c, the coefficients in the matrix representation in the Chebyshev basis of (I + K) -1 converge exponentially fast to 0 when moving off a few diagonals or upper rows, in accordance with the logarithmic bound for N R claimed in Theorem 3.16. On the contrary, it also appears that finite-dimensional truncations are, by far, less efficient approximations. This phenomenon is reflected by the minimum value for Newton-Galerkin's N G and Newton-Picard's N R needed to ensure a contracting fixed-point operator. More precisely, Figure 3d compares the evolution in function of [x l , x r ] = [-a, a] of N G (Newton-Galerkin), N R (Newton-Picard), and N P , the number of Picard iterations needed in the algorithm of [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF]. In addition, we provide the a priori conservative estimate for N R from Corollary 3.17, computed over a Bernstein ellipse with ρ = 2 (the estimate can be made effective rather than asymptotic by using the bounds in the proof). The first observation is that Newton-Galerkin validation method can hardly run with a 10, due to the exponential blow-up of the truncation index N G . This indeed requires to perform rigorous computations with matrices with several billions entries... The second observation is that although N R and N P should be similar in theory (which is reflected by the similar behavior of N P and the a priori estimate for N R ), in practice N P is significantly larger, because it is computed as a conservative, purely a priori estimate, contrary to N R where the actual contraction ratio of T is explicitly computed. Timing Comparison. To conclude the Airy example, we compare timings of computing RPAs in the Chebyshev basis using several a posteriori validation methods:

• The Newton-Picard validation algorithm NewtonPicardValid presented in this article and implemented in the C library ChebValid. • The Newton-Galerkin validation method detailed in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF], for which the ChebValid library was initially designed. • The Picard iteration based a posteriori validation algorithm of [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF], which we implemented (without any optimization) in ChebValid.

Although the application scope of self-validating methods slightly differs (the validation method constructs its own approximant, which is no longer freely provided by the user), we also compared for the sake of completeness with two implementations compatible with multiprecision:

• The CAPD library9 computes rigorous Taylor forms with remainders using a higher-order Lohner method [START_REF] Wilczak | C r -Lohner algorithm[END_REF]. The result is a piecewise analytic representation of the solution, as opposed to uniform polynomial approximants considered in this article. • The algorithm developed in [START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF] and implemented in the SageMath ore_algebra package 10combines a majorant series approach with a binary splitting algorithm to evaluate pointwise D-finite functions. The scope is even more restricted, but we fairly acknowledge the remarkable performances of this implementation.

We first note that in the category of Newton-like a posteriori validation methods, the Newton-Picard algorithm presented in this paper outperforms Newton-Galerkin and Picard iterations, which are both unable to deal with the instances where a 10 in reasonable time. For the former, the main reason is that the required truncation index N G becomes very large and rapidly exceeds the memory available on a standard machine. For the latter, the number of iterations remains reasonable (although larger than the degree N R used in Newton-Picard), but the iteration process using interval arithmetics suffers from numerical instability, which in the simple case a = 5 forced us to use MPFR floating-point numbers with 20, 000 bits of precision.

The comparison with self-validating methods used in CAPD or SageMath is more delicate, since those implementations are by far more mature than our prototypes. However, the timings highlight the "cost to pay" to construct a contracting fixed-point operator as a common feature of a posteriori validation methods. For very high precision, this constant cost is amortized by the linear asymptotic complexity w.r.t. the degree of the candidate approximation. In particular, NewtonPicardValid performs better than CAPD when high precision is required. Note however that we did not manage to use CAPD with Taylor models of degree more than 70, which may explain the poorer performance of this reputed library in high precision.

Finally, when it comes to the simple pointwise rigorous evaluation of Ai, the algorithm implemented in SageMath outperforms all the other methods. One reason for this is that it uses a binary splitting technique to evaluate the truncated Taylor series without computing explicitly the Taylor coefficients up to the truncation order. This is quite different from Problem 1.1 where we want to compute an RPA with an explicit polynomial approximation that can be manipulated. However, some techniques used in this implementation could certainly be reused to improve the existing algorithms for RPAs.

A Boundary Value Problem with Non-Polynomial Coefficients

As mentioned in the introduction, Algorithm NewtonPicardValid easily extends to the case of linear ODEs with non-plynomial coefficients, represented using RPAs. Indeed, one just needs to provide routines to manipulate low-rank kernels defined with univariate RPAs rather than pure polynomials, namely evaluation on another RPA, composition and rigorous bound. This is implemented in the ChebValid library, with the prefix chebmodel_ replacing mpfi_chebpoly_ in the names of the corresponding functions.

To illustrate this, let us consider the following boundary layer problem considered in [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] and reused in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]:

   y (x) - 2x ε cos x - 8 10 y (x) + 1 ε cos x - 8 10 y(x) = 0, x ∈ [-1, 1], y(-1) = y(1) = 1. (25)
This is a Boundary Value Problem (BVP) rather than an Initial Value Problem (IVP). However, since the problem is linear, this amounts to solving two IVPs on the same ODE, namely computing y = y 0 + λy 1 with y 0 , y 1 solutions to the IVP with initial conditions (1, 0) and (0, 1), respectively, and with λ = (1 -y 0 (1))/y 1 (1). Table 2. Timings comparison for the rigorous computation of Ai over [-a, a] with prescribed target accuracy, between Newton-Picard (N-P), Newton-Galerkin (N-G), Picard iterations (Pit) and Taylor forms in the CAPD library. Timings for the rigorous evaluation of Ai(a) using the SageMath ore_algebra package are also provided. "mpfr prec" denotes the MPFR precision in bits and "deg" the degree of the Chebyshev polynomial approximation used in (N-P), (N-G) and (Pit).

We first use Algorithm NewtonPicardValid to compute an RPA for cos over [-1, 1], using the differential equation y + y = 0. This allows us to provide RPAs of arbitrary accuracy for the two coefficients in [START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF]. Then we call this algorithm again on the resulting IVP, overloaded with RPAs. The smaller ε is, the steeper the slope of the solution y is, as it can be observed in Figure 4a. This results in an increasing degree N R to approximate the resolvent kernel (see Figure 4b). Note however that N R is way smaller than the truncation order N G required by Newton-Galerkin validation method (see [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]Sec. 7.3]). As a consequence, we can rigorously deal with instances with really small ε (up to 10 -6 within a few seconds), while Newton-Galerkin is limited to ε = 0.005 roughly.

On a last note, one remarks that y 0 and y 1 have very high amplitude compared to the BVP solution y. With the algorithm of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF], one can in fact directly numerically solve this BVP by replacing the first two rows of the almost-banded system by the corresponding equations representing the boundary conditions. An interesting future work is to extend Newton-Picard validation algorithm to BVPs, by investigating the resolvent kernel of the associated integral equation.

High Order D-finite Equations

In order to check that the proposed method successfully handles linear ODEs of order greater than 2, 21 random D-finite equations L 10 , L 11 , . . . , L 30 were generated, with order r i = i and degree 2 polynomial coefficients with integer coefficients uniformly chosen between -100 and 100, like the following one: Despite the high order and dense coefficients of the L i , Algorithm NewtonPicardValid requires quite moderate degrees N R to construct the Newton-Picard fixed-point operator (always less than 700 for those 21 randomly generated instances). Hence, we see in practice that the remarkable approximation properties of Chebyshev spectral methods are reflected in the validation procedure. Examples of "natural" high order ODEs are hard to find in the literature and applications. Yet, they frequently appear in the context of D-finite functions, e.g. when performing closure operations like addition, multiplication or algebraic substitution on simpler functions, and their properties may be quite different from those of random equations. Consider for example the D-finite function f (x) = cos x + cos 5x 2 -sin 3x 4 over [0, 2], plotted in Figure 5a. The holexprtodiffeq command of Maple's Gfun package [START_REF] Salvy | Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable[END_REF] automatically computes an order 6, degree 37 differential equation L{f } = 0, with a 6-dimensional solution space spanned by {cos x, sin x, cos 5x 2 , sin 5x 2 , cos 3x 4 , sin 3x 4 }. Unfortunately, the leading coefficient, a 6 (x) = 268738560000x 29 -5374771200x 27 -365158969344x 25 + 3732480000x 23 + 463694284800x 21 -6366325248x 19 -66390306048x 17 + 2217859200x 15 + 27057529296x 13 -2477301120x 11 + 3908641720x 9 -134485400x 7 + 1100233x 5 + 15403262x 3 , is non-constant, and in particular singular at x = 0. To circumvent this difficulty, the Desingularize command can be used to compute a non-singular, higher order equation. The result is an operator of order 9 and degree 35, with significantly larger coefficients, which is generated in approximately 10 seconds in Maple. Our experiments on this equation show the following:

• The Chebyshev spectral method provides remarkably accurate numerical approximations of this function, despite the considerable size of the equation. Indeed, the Chebyshev coefficients (see Figure 5b) converge at a super-linear rate, since f is entire. • Unfortunately, NewtonPicardValid fails to validate this example in reasonable time, due to a too large degree N R needed to approximate the resolvent kernel so as to get a contracting operator. Indeed, the functions ϕ i and ψ i defining R(x, t) are way larger than f , and cannot be accurately approximated with a moderate degree. The reason is that the desingularization step produces a much bigger and ill-conditioned equation (in this example, max 0 i<9 a i Ч 1 = 1.26 10 22 ). Hence, the extra 3 functions in the new 9-dimensional solution space most probably have Chebyshev coefficients of size exponential in this bound, and are therefore not even representable in MPFR.

The key lesson of this example is that while the Chebyshev spectral method can perform well on specific initial conditions associated to a well-behaved particular function of the solution space of the differential equation, the Newton-Picard validation method needs to inspect the whole solution space (of L and L * ) to guarantee the maximum error "in every direction". Still, future improvements could help in taming the contribution of the ill-conditioned but non-relevant solutions. Also, in this example and more generally for D-finite functions, handling directly singular equations is an interesting future challenge.

Discussion and Future Work

Newton-Galerkin a posteriori validation methods have been widely used over recent years to provide successful computer-assisted proofs for the existence (and the effective approximation) of solutions to various kinds of functional equations, notably (linear or nonlinear) ODEs with initial or boundary value conditions, partial differential equations, delay-differential equations, etc. Elaborating on the simpler case of linear ODEs with initial values, we presented a new validation algorithm (which we propose to call Newton-Picard) and provided a thorough complexity analysis to explain why it performs better than other a posteriori validation algorithms in the literature. This algorithm is implemented in the prototype C library ChebValid.

This work is to be seen as part of a long term effort to combine the field of rigorous numerics towards computer-assisted proofs in mathematics with the algorithmic and complexity approach of computer algebra, in order to design efficient and fully trusted symbolic-numeric libraries for function space computations. In particular, several future goals can be identified:

• Combining the Newton-Picard fixed-point operator with the radii polynomial approach of [START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach[END_REF] to bound nonlinear terms will allow for the efficient treatment of nonlinear ODEs, where even proving the existence of a solution on a given interval is a nontrivial problem. In particular, investigating the complexity and comparing with the results of [START_REF] Bournez | On the complexity of solving initial value problems[END_REF] is an interesting challenge. • Deriving formulas for the resolvent kernel in the case of boundary conditions is the key to adapting this algorithm to BVPs. • Even in the case of D-finite functions, dealing with singular points (that is, when the leading coefficient, instead of being 1, vanishes at some point of the interval) is a difficult problem. While rigorous evaluation of D-finite functions near regular singular points has already been investigated, e.g. in [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions D-finies[END_REF], rigorous uniform approximation (with a suitable norm) is, up to our knowledge, still an open problem. A similar question is the case of unbounded intervals. • Computing rigorous approximants for PDE solutions is a major challenge due to the (time and space) complexity issues related to the several variable case. Hence, proposing efficient and reliable algorithms is crucial in the development of computer-assisted proofs for such problems.

• Finally, we propose to implement this validation algorithm in the Coq proof assistant in order to provide a fully reliable tool to engineers and mathematicians. More specifically, we want to integrate it into our experimental development of certified RPAs in Coq using an a posteriori (or certificate-based) approach [START_REF] Bréhard | A certificate-based approach to formally verified approximations[END_REF]. The algebraic framework of Newton-Picard with integral operators and polynomial kernels makes this new validation algorithm a relevant choice for that task.

Figure 1 .

 1 Figure 1. Almost-banded matrix representations in the Chebyshev basis of an integral operator K and its truncation K [N ] .

Remark 3 . 6 .Lemma 3 . 7 .

 3637 It computes a Newton-Picard fixed-point operator T and a rigorous contraction ratio λ in O(N 2 R r 2 ) arithmetic operations. • Given T and λ, it computes for any right-hand side g and any approximation f • a rigorous error bound for f • -f * in O((N R + n + m)N R r) arithmetic operations, where n = deg f • and m = deg g. Theorem 3.5 gives a separate complexity estimate for Step 3, since Steps 1 and 2 are completely independent of f • and g. This means that the same Newton-Picard fixed-point validation operator T can be reused to validate another approximation f • in Problem 1.1, possibly with a different right-hand side h and/or different initial conditions v i . If K is an integral operator with polynomial kernel K(x, t) = k-1 i=0 α i (x)β i (t) of rank k and total degree n, and f (x) a univariate polynomial of degree m, then K{f } = x x0 K(x, t)f (t)dt is a polynomial of degree at most n + m + 1 that can be computed in O(kn(n + m)) arithmetic operations.Proof. For each 0 i < k, α i (x)x x0 β i (t)f (t)dt is a univariate polynomial of degree at most deg α i + deg β i + deg f + 1 n + m + 1,that requires O(deg β i deg f ) arithmetic operations for the inner multiplication, O(deg β i + deg f ) ones for the integration, and finally O(deg α i (deg β i + deg f )) ones for the outer multiplication, giving a total of O(n(n + m)) arithmetic operations. The claimed overall complexity follows from the iteration for each 0 i < k and the additions. Algorithm 2 KernelComp(K, L)

3. 3 . 1 .

 31 Formulas for the Resolvent Kernel. The scalar linear ODE (3) can be transformed into the equivalent r-dimensional first-order vectorial linear ODE Y (x) = A(x)Y (x) + H(x), where:

( a )

 a Airy functions Ai and Bi. (b) Chebyshev coefficients of y = Ai for a = 15. (I + K) -1 with K in (24) and a = 15. (d) Comparison of Newton-Picard's N R (real value and a priori estimate) with Newton-Galerkin's NG and the number NP of Picard iterations in [2].

Figure 3 .

 3 Figure 3. Validation of the Airy function Ai over [x l , x r ] = [-a, a].

L 10 =

 10 ∂ 10 + (-93x 2 + 5x + 73)∂ 9 + (-51x 2 + 22x + 74)∂ 8 + (-8x 2 -41x + 41)∂ 7 + (-18x 2 -76x -11)∂ 6 + (79x 2 -45x + 70)∂ 5 + (-11x 2 + 61x -83)∂ 4 + (-74x 2 -5x -84)∂ 3 + (-x 2 -63x + 96)∂ 2 + (-95x 2 + 86x -57)∂ + (-8x 2 + 21x -56),together with random integer initial conditions between -100 and 100.

( a )

 a Solution for different values of (b) Required degree N R in function of ε.

Figure 4 .

 4 Figure 4. Validation of the boundary layer example using NewtonPicardValid.

Figure 5 .

 5 Figure 5. An example of a D-finite function with high-order, ill-conditioned differential equation.

  independent variable • ∞ uniform norm over [x l , xr] [x l , xr] real compact interval of interest C 0 space of continuous functions over [x l , xr] x0 initial time in [x l , xr] RN [x] space of univariate polynomials with real T = max(xr -x0, x0 -x l ) maximum time coefficients of degree N ∂ differentiation operator w.r.t. x Tn n-th Chebyshev polynomial of the first L linear differential operator kind, rescaled over [x l , xr] L

	x

* adjoint differential operator

  1, rewriting K(x, t) given in Propositions 2.3 or 2.5 as a rank r expression (line 1) requires O(r 2 s) arithmetic operations. Then, anticipating the complexity estimate for ResolventKernel given in Proposition 3.13, line 2 requires O(N R (r+s) 2 ) arithmetic operations. Moving on to step 2, the kernel composition in line 3 runs in complexity O((N R + r + s) 2 r 2 ), according to Lemma 3.8. After that, the computation of a bound λ for E(x, t) (lines 4 and 5) runs in O((N and O(2N R (2N R + r + s + n + m)r) arithmetic operations, using Lemma 3.7. Finally, the bound

R + r + s) 2 r). This justifies the overall complexity for steps 1 and 2, assuming r + s = O(N R ).

Concerning step 3, lines 7 and 8 computing δ(x) respectively require O(r(r + s)(r + s + n) + m)

https://www.chebfun.org/

https://gitlab.inria.fr/brisebar/tchebyapprox/-/releases/v0.0.1

see, e.g., https://www.texmacs.org/joris/stablemult/stablemult-abs.html

Refined complexity estimates are given in[START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] Prop.

5.3] when using an almost-banded A instead of a dense one.

The adjunction relation requires the kernels to be defined over [x l , xr] 2 rather than [x l , xr|x 0 ] only, which is always the case in this article.

https://gitlab.inria.fr/brisebar/tchebyapprox/-/releases/v0.0.1

https://www.mpfr.org/

https://gforge.inria.fr/projects/mpfi/

http://capd.ii.uj.edu.pl/

https://github.com/mkauers/ore_algebra/

Acknowledgments. I would like to thank the anonymous referees for their careful reviewing work and their useful suggestions, which helped improve the clarity of this article.