Polynomial interpolation and a priori bootstrap for computer-assisted proofs in nonlinear ODEs

TL;DR

This paper introduces a rigorous computational method using polynomial interpolation and a priori bootstrap to prove existence of solutions and special orbits in nonlinear ODEs, demonstrated on Lorenz and ABC flow systems.

Contribution

It develops a novel a priori bootstrap technique combined with a fixed point approach for computer-assisted proofs in nonlinear ODEs.

Findings

Proved existence of solutions for Lorenz system at classical parameters.

Established existence of ballistic spiral orbits in ABC flows.

Validated the method on complex nonlinear systems.

Abstract

In this work, we introduce a method based on piecewise polynomial interpolation to enclose rigorously solutions of nonlinear ODEs. Using a technique which we call a priori bootstrap, we transform the problem of solving the ODE into one of looking for a fixed point of a high order smoothing Picard-like operator. We then develop a rigorous computational method based on a Newton-Kantorovich type argument (the radii polynomial approach) to prove existence of a fixed point of the Picard-like operator. We present all necessary estimates in full generality and for any nonlinearities. With our approach, we study two systems of nonlinear equations: the Lorenz system and the ABC flow. For the Lorenz system, we solve Cauchy problems and prove existence of periodic and connecting orbits at the classical parameters, and for ABC flows, we prove existence of ballistic spiral orbits.

Figures

Tables

Table 1. Table 1: Comparisons for p = 1 𝑝 1 p=1 . τ m ​ a ​ x subscript 𝜏 𝑚 𝑎 𝑥 \tau_{max} is the longest integration time for which the proof succeeds, and r 𝑟 r is the associated validation radius, that is a bound of the distance (in 𝒞 0 superscript 𝒞 0 \mathcal{C}^{0} norm) between the numerical data used and the true solution.

$k$	$m$	$nm\left(k+1\right)$	${\tau }_{max}$	$r$
1	2333	13998	0.69	$2.3233\times {10}^{-5}$
2	1556	14004	0.64	$1.1524\times {10}^{-7}$
3	1167	14004	0.58	$8.6805\times {10}^{-9}$

Table 2. Table 2: Comparisons for p = 2 𝑝 2 p=2 . τ m ​ a ​ x subscript 𝜏 𝑚 𝑎 𝑥 \tau_{max} is the longest integration time for which the proof succeeds, and r 𝑟 r is the associated validation radius, that is a bound of the distance (in 𝒞 0 superscript 𝒞 0 \mathcal{C}^{0} norm) between the numerical data used and the true solution.

$k$	$m$	$nm\left(k+1\right)$	${\tau }_{max}$	$r$
1	2333	13998	0.97	$1.5718\times {10}^{-3}$
2	1556	14004	5.6	$8.4373\times {10}^{-5}$
3	1167	14004	5.5	$8.4184\times {10}^{-8}$
4	933	13995	4.9	$7.9190\times {10}^{-9}$

Table 3. Table 3: Comparisons for p = 3 𝑝 3 p=3 . τ m ​ a ​ x subscript 𝜏 𝑚 𝑎 𝑥 \tau_{max} is the longest integration time for which the proof succeeds, and r 𝑟 r is the associated validation radius, that is a bound of the distance (in 𝒞 0 superscript 𝒞 0 \mathcal{C}^{0} norm) between the numerical data used and the true solution.

$k$	$m$	$nm\left(k+1\right)$	${\tau }_{max}$	$r$
2	1556	14004	5.6	$7.6716\times {10}^{-4}$
3	1167	14004	8.1	$9.3043\times {10}^{-6}$
4	933	13995	8.1	$8.8204\times {10}^{-8}$
5	778	14004	8.0	$1.6175\times {10}^{-8}$
9	467	14010	7.9	$1.3748\times {10}^{-8}$
19	233	13980	6.9	$2.2998\times {10}^{-8}$

Table 4. Table 4: Minimal number of coefficients needed to validate the orbit of length τ = 2 𝜏 2 \tau=2 , starting from u 0 subscript 𝑢 0 u_{0} given in ( 22 ).

	$k=1$	$k=2$	$k=3$	$k=4$
$p=2$	no proof	$m=416$	$m=415$	$m=377$
	no proof	$nm\left(k+1\right)=3744$	$nm\left(k+1\right)=4980$	$nm\left(k+1\right)=5655$
	$k=2$	$k=3$	$k=4$	$k=5$
$p=3$	$m=470$	$m=125$	$m=110$	$m=99$
	$nm\left(k+1\right)=4230$	$nm\left(k+1\right)=1500$	$nm\left(k+1\right)=1650$	$nm\left(k+1\right)=1782$

Table 5. Table 5: The intervals [ τ A − , τ A + ] subscript superscript 𝜏 𝐴 subscript superscript 𝜏 𝐴 [\tau^{-}_{A},\tau^{+}_{A}] , for A = 0.1 , 0.2 , … , 1 𝐴 0.1 0.2 … 1 A=0.1,0.2,\ldots,1 in which the period τ A subscript 𝜏 𝐴 \tau_{A} of the solution described in Theorem  7.1 is proved to be.

$A$	$1$	$0.9$	$0.8$	$0.7$	$0.6$
${\tau }_A^{-}$	$3.23527736$	$3.41779635$	$3.62512508$	$3.86405419$	$4.14464726$
${\tau }_A^{+}$	$3.23527746$	$3.41779647$	$3.62512521$	$3.86405436$	$4.14464749$

Equations

$$\{\begin{array}{rl}& \frac{du}{dt}(t)=\varphi (u(t)),t\in [0,\tau ],\\ & BV(u(0),u(\tau ))=0.\end{array}$$

$$\{\begin{array}{rl}& \frac{du}{dt}(t)=\varphi (u(t)),t\in [0,\tau ],\\ & u(0)=u_0,\end{array}$$

$$f:\{\begin{array}{rl}{\mathcal{C}}^0([0,1],{\mathbb{R}}^n)& \to {\mathcal{C}}^1([0,1],{\mathbb{R}}^n)\\ u& \mapsto f(u),\end{array}$$

$$f(u)(t)\stackrel{\text{mbox}[}{=}{u}_0+\tau {\int }_0^t\varphi (u(s))ds,\text{for all}\text{leavevmode}t\in [0,1].$$

$${\varphi }^{[0]}(u)\stackrel{\text{mbox}[}{=}u\text{and}{\varphi }^{[p+1]}(u)\stackrel{\text{mbox}[}{=}D{\varphi }^{[p]}(u)\varphi (u)\text{for all }u\in {\mathbb{R}}^n\text{ and for all }p=0,\dots ,\mathbf{r}.$$

$$\{\begin{array}{rl}& \frac{d^{p}u}{d{t}^p}(t)={\varphi }^{[p]}(u(t)),t\in [0,\tau ],\\ & \frac{d^{q}u}{d{t}^q}(0)={\varphi }^{[q]}(u_0),\text{for all }q=0,\dots ,p-1.\end{array}$$

$$\tilde{g}:\{\begin{array}{rl}{\mathcal{C}}^0([0,1],{\mathbb{R}}^n)& \to {\mathcal{C}}^p([0,1],{\mathbb{R}}^n)\\ u& \mapsto \tilde{g}(u),\end{array}$$

$$\tilde{g}(u)(t)\stackrel{\text{mbox}[}{=}\sum _{q=0}^{p-1}{\tau }^q\frac{t^q}{q!}{\varphi }^{[q]}(u_0)+{\tau }^p{\int }_0^t\frac{(t-s{)}^{p-1}}{(p-1)!}{\varphi }^{[p]}(u(s))ds,\text{for all}\text{leavevmode}t\in [0,1].$$

$${\Delta }_m\stackrel{\text{mbox}[}{=}\{t_0,t_1,\dots ,t_m\},$$

$$g:\{\begin{array}{rl}{\mathcal{C}}_{{\Delta }_m}^0([0,1],{\mathbb{R}}^n)& \to {\mathcal{C}}_{{\Delta }_m}^p([0,1],{\mathbb{R}}^n)\\ u& \mapsto g(u),\end{array}$$

$$g(u)(t)\stackrel{\text{mbox}[}{=}\sum _{q=0}^{p-1}{\tau }^q\frac{(t-t_j{)}^q}{q!}{\varphi }^{[q]}(u(t_j^{-}))+{\tau }^p{\int }_{t_j}^t\frac{(t-s{)}^{p-1}}{(p-1)!}{\varphi }^{[p]}(u(s))ds,$$

$$G(u)\stackrel{\text{mbox}[}{=}g(u)-u.$$

$$u(t_j^{-})=\underset{t\to t_j^{+}}{\lim }g(u)(t)=\underset{t\to t_j^{+}}{\lim }u(t)=u(t_j^{+}).\text{qed}$$

$$t_{j,l}\stackrel{\text{mbox}[}{=}{t}_j+\frac{x_l^k+1}{2}(t_{j+1}-t_j),\text{for }l=1,\dots ,k-1,$$

$$x_l^k\stackrel{\text{mbox}[}{=}\cos {\theta }_l^k,{\theta }_l^k\stackrel{\text{mbox}[}{=}\frac{k-l}{l}\pi ,\text{for }l=0,\dots ,k.$$

$$S_{m,k}^n\stackrel{\text{mbox}[}{=}\left\{u\in {\mathcal{C}}_{{\Delta }_m}^0([0,1],{\mathbb{R}}^n):u_{{|}_{(t_j,t_{j+1})}}\text{ is a polynomial of degree }k\text{ for all }j=0,\dots ,m-1\right\}.$$

$${\Pi }_{m,k}^n:\{\begin{array}{rl}{\mathcal{C}}_{{\Delta }_m}^0([0,1],{\mathbb{R}}^n)& \to S_{m,k}^n\\ u& \mapsto \bar{u}={\Pi }_{m,k}^n(u),\end{array}$$

$$u(t_j^{-})=\bar{u}(t_j^{-})\text{and}u(t_j^{+})=\bar{u}(t_j^{+}).$$

$${\Vert (Id-{\Pi }_{m,k}^1)u\Vert }_{\infty }={\Vert u-\bar{u}\Vert }_{\infty }\le C_k\underset{0\le j<m}{\max }\left((t_{j+1}-t_j{)}^{k+1}\underset{t\in [t_j,t_{j+1}]}{\max }\left|\frac{d^{k+1}u}{d{t}^{k+1}}(t)\right|\right),$$

$$C_k\stackrel{\text{mbox}[}{=}\frac{1}{(k+1)!2^{2k}}.$$

$${\Vert (Id-{\Pi }_{m,k}^1)u\Vert }_{\infty }={\Vert u-\bar{u}\Vert }_{\infty }\le {\tilde{C}}_{k,l}\underset{0\le j<m}{\max }\left((t_{j+1}-t_j{)}^l\underset{t\in [t_j,t_{j+1}]}{\max }\left|\frac{d^{l}u}{d{t}^l}(t)\right|\right),$$

$${\tilde{C}}_{k,l}\stackrel{\text{mbox}[}{=}\min \left[\left(1+{\Lambda }_k\right){\left(\frac{\pi }{4}\right)}^l\frac{(k+1-l)!}{(k+1)!},\frac{1}{l!2^l}\sum _{q=0}^{\left[\frac{l-1}{2}\right]}\frac{1}{4^q}\left(\genfrac{}{}{0ex}{}{l-1}{2q}\right)\left(\genfrac{}{}{0ex}{}{2q}{q}\right)\right],$$

$${\mathcal{E}}_{m,k}^n\stackrel{\text{mbox}[}{=}\left\{(i,j,l)\in {\mathbb{N}}^3,1\le i\le n,0\le j\le m-1,0\le l\le k\right\}.$$

$${\bar{u}}^{(i)}(t)=\sum _{l=0}^k{\bar{u}}_{j,l}^{(i)}{T}_l\left(\frac{t-t_j}{t_{j+1}-t_j}-\frac{t_{j+1}-t}{t_{j+1}-t_j}\right),\text{for all }j=0,\dots ,m-1\text{ and }t\in (t_j,t_{j+1}),$$

$$G(u)\left(t_{j,l}\right)=\sum _{q=0}^{p-1}{\tau }^q\frac{(t_{j,l}-t_j{)}^q}{q!}{\varphi }^{[q]}(u(t_j^{-}))+{\tau }^p{\int }_{t_j}^{t_{j,l}}\frac{(t_{j,l}-s{)}^{p-1}}{(p-1)!}{\varphi }^{[p]}(u(s))ds-u(t_{j,l}),$$

$$\text{for all}\text{leavevmode}j=0,\dots ,m-1\text{ and}\text{leavevmode}l=0,\dots ,k.$$

$${\mathcal{X}}_{m,k}^n\stackrel{\text{mbox}[}{=}{S}_{m,k}^n\text{and}{\mathcal{X}}_{\infty }^n\stackrel{\text{mbox}[}{=}(Id-{\Pi }_{m,k}^n){\mathcal{C}}_{{\Delta }_m}^0([0,1],{\mathbb{R}}^n).$$

$${\Pi }_{m,k}^n:\{\begin{array}{rl}{\mathcal{X}}^n& \to {\mathcal{X}}_{m,k}^n\\ u& \mapsto \bar{u}={\Pi }_{m,k}^n(u),\end{array}$$

$${\Pi }_{\infty }^n:\{\begin{array}{rl}{\mathcal{X}}^n& \to {\mathcal{X}}_{\infty }^n\\ u& \mapsto {\Pi }_{\infty }^n(u)=u-\bar{u}=(Id-{\Pi }_{m,k}^n)(u).\end{array}$$

$${\Vert {\Pi }_{m,k}^n(u)\Vert }_{{\mathcal{X}}_{m,k}^n}\stackrel{\text{mbox}[}{=}\underset{(i,j,l)\in {\mathcal{E}}_{m,k}^n}{\max }\left|{\bar{u}}^{(i)}(t_{j,l})\right|\text{and}{\Vert {\Pi }_{\infty }^n(u)\Vert }_{{\mathcal{X}}_{\infty }^n}\stackrel{\text{mbox}[}{=}\underset{1\le i\le n}{\max }{\Vert u^{(i)}-{\bar{u}}^{(i)}\Vert }_{\infty }.$$

$${\Vert u\Vert }_{{\mathcal{X}}^n}\stackrel{\text{mbox}[}{=}\max \left({\Vert {\Pi }_{m,k}^n(u)\Vert }_{{\mathcal{X}}_{m,k}^n},\frac{1}{r_{\infty }}{\Vert {\Pi }_{\infty }^n(u)\Vert }_{{\mathcal{X}}_{\infty }^n}\right),$$

Full text

Key words Computer-assisted proof $\cdot$ Nonlinear ODEs $\cdot$ Picard-like operator

Piecewise polynomial interpolation $\cdot$ Newton-Kantorovich $\cdot$ ABC flows

2010 AMS Subject Classification 65G20 $\cdot$ 65P99 $\cdot$ 65D30 $\cdot$ 37M99 $\cdot$ 37C27

1 Introduction

This paper introduces an approach based on polynomial interpolation to obtain mathematically rigorous results about existence of solutions of nonlinear ordinary differential equations (ODEs). Our motivation for the present work is threefold. First, we believe that polynomial interpolation techniques are versatile and can lead to efficient and general computational methods to approximate solutions of ODEs with complicated (non polynomial) nonlinearities. Second, while polynomial interpolation techniques have be used to produce computer-assisted proofs in ODEs, their applicability to produce proofs is sometimes limited by the formulation of the problem itself. More precisely, a standard way to prove (both theoretically and computationally) existence of solutions of systems of ODEs is to reformulate the problem into an integral equation (often in the form of a Picard operator) and then to apply the contraction mapping theorem to get existence. If one is interested to produce computer-assisted proofs using that approach, the analytic estimates required to perform the proofs depend on the amount of regularity gained by applying the integral operator. This observation motivated developing what we call the a priori bootstrap, which consists of reformulating the original ODE problem into one of looking for the fixed point of a higher order smoothing Picard-like operator. Third, we believe (and hope) that our proposed method can be adapted to study infinite dimensional continuous dynamical systems (e.g. partial differential equations and delay differential equations) for which spectral methods may sometimes be difficult to apply (for instance because of the shape of the spatial domains or because the differential operators are difficult to invert in a given spectral basis).

It is important to realize that computer-assisted arguments to study differential equations are by now standard, and that providing a complete overview of the literature would require a tremendous effort and is outside the scope of this paper. However, we encourage the reader to read the survey papers [1, 2, 3, 4, 5, 6], the books [7, 8] and to consult the webpage of the CAPD group [9] to get a flavour of the extraordinary recent advances made in the field.

More closely related to the present work are methods based on the contraction mapping theorem using the radii polynomial approach (first introduced in [10]), which has been developed in the last decade to study fixed points, periodic orbits, invariant manifolds and connecting orbits of ODEs, partial differential equations and delay differential equations (see for instance [11, 12, 13, 14, 15, 16]). The numerics and a posteriori analysis in those works mainly use spectral methods like Fourier and Chebyshev series, and Taylor methods. First order polynomial (piecewise linear) approximations were also used using the radii polynomial approach (see [17, 18, 19]), but more seldom, mainly because the numerical cost was higher and the accuracy was lower than for spectral methods. The computational cost of these low order methods is essentially due to the above mentioned low gain of regularity of the Picard operators chosen to perform the computer-assisted proofs.

In an attempt to address the low gain of regularity problem, we present here a new technique that we call a priori bootstrap which, when combined with the use of higher order interpolation, significantly improves the efficiency of computer-assisted proofs with polynomial interpolation methods. We stress that the limitations that affected the previous works using interpolation were not solely due to the use of first order methods, and that the a priori bootstrap is crucial (that is, just increasing the order of the interpolation does not significantly improve the results in those previous works). This point is illustrated in Section 6.

While we believe that one of the advantage of our proposed method is the versatility of the polynomial interpolations to tackle problems with complicated (non polynomial) nonlinearities, we hasten to mention the existence of previous powerful methods which have been developed in rigorous computing to study such problems. For instance, automatic differentiation (AD) techniques provide a beautiful and efficient means of computing solutions of nonlinear problems (e.g. see [7, 20, 21]) and are often combined with Taylor series expansions to prove existence of solutions of differential equations with non polynomial nonlinearities (e.g. see [1, 22, 23, 24, 25, 26, 27]). Also, in the recent work [28], the ideas of AD are combined with Fourier series to prove existence of periodic orbits in non polynomial vector fields. Independently of AD techniques, a method involving Chebyshev series to approximate the nonlinearities have been proposed recently in [29]. Finally, the fast Fourier transform (FFT) algorithm is used in [30] to control general nonlinearities in the context of computer-assisted proofs in KAM theory.

In this paper we consider $\phi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ a $C^{{{\boldsymbol{r}}}}$ vector field (not necessarily polynomial) with ${{{\boldsymbol{r}}}}\geq 1$, and we present a rigorous numerical procedure to study problems of the form

[TABLE]

We treat three special cases for $BV$, corresponding to an initial value problem, a problem of looking for periodic orbits and a problem of looking for connecting orbits. We also note that, as for the already existing spectral methods, the technique presented here extends easily to treat parameter dependent versions of (1) (e.g. using ideas from [31, 32, 33]). For the sake of simplicity, we expose all the general arguments for the initial value problem only, that is for

[TABLE]

given an integration time $\tau>0$ and an initial condition $u_{0}\in\mathbb{R}^{n}$. We explain how (2) needs to be modified for different problems as we introduce them in Sections 6 and 7.

Our paper is organized as follows. In Section 2, we start by presenting our a priori bootstrap technique, together with a piecewise reformulation of the operator that we use throughout this work. We then recall in Section 3 some definitions and error estimates about polynomial interpolation, and explain how to combine them with our a priori bootstrap formulation to get computer-assisted proofs. The precise estimates needed for the proofs are then derived in Section 4, and their dependency with respect to the a priori bootstrap and to the parameters of the polynomial interpolation is commented in Section 5. This discussion is complemented by several examples in Section 6, where we apply our technique to validate solutions for the Lorenz system. Finally we give another example of application in Section 7, where we prove the existence of some specific orbits for ABC flows.

2 Reformulations of the Cauchy problem

2.1 A priori bootstrap

One of the usual strategies used to study (2), both theoretically and numerically, is to recast it as a fixed point problem, as in the following lemma.

Lemma 2.1**.**

Consider the standard Picard operator

[TABLE]

where

[TABLE]

Then $u$ is a fixed point of $f$ if and only if $v:t\mapsto u(\frac{t}{\tau})$ is a solution of (2).

In previous works using this reformulation, the limiting factor in the estimates needed to apply the contraction mapping theorem was a consequence of the fact that $f$ only gains one order of regularity, that is maps $\mathcal{C}^{0}$ into $\mathcal{C}^{1}$. This fact will be made precise in Section 4 where we derive the estimates in question and in Section 5 where we discuss how those estimates affect the effectiveness of our technique.

To circumvent this limitation, we propose a different reformulation that we call a priori bootstrap. This approach provides operators which gain more regularity, and therefore lead to sharper analytic estimates. First we introduce some notations. The following definition allows concisely describing the higher order equations obtained by taking successive derivatives of (2).

Definition 2.2**.**

Consider the sequence of vector fields $\left(\phi^{[p]}\right)_{0\leq p\leq{{\boldsymbol{r}}}+1}$ with $\phi^{[p]}:\mathbb{R}^{n}\to\mathbb{R}^{n}$,

[TABLE]

Lemma 2.3**.**

For any $1\leq p\leq{{\boldsymbol{r}}}+1$, $u$ solves (2) if and only if $u$ solves the following Cauchy problem

[TABLE]

Proof.

The direct implication is trivial. To prove the converse application, we consider $e\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\frac{du}{dt}-\phi(u)$ and show that is solve a linear ODE of order $p-1$, with initial conditions $\frac{d^{q}e}{de^{q}}(0)=0$ for all $q=0,\dots,p-2$, which implies that $e\equiv 0$. ∎

Integrating the $p^{th}$ order Cauchy problem (4) $p$ times leads to a new fixed point operator which now maps $\mathcal{C}^{0}$ into $\mathcal{C}^{p}$.

Lemma 2.4**.**

Let $1\leq p\leq{{\boldsymbol{r}}}+1$ and consider the Picard-like operator

[TABLE]

where

[TABLE]

Then $u$ is a fixed point of $\tilde{g}$ if and only if $v:t\mapsto u(\frac{t}{\tau})$ is a solution of (4) (and thus of (2)).

Proof.

If $u$ is a fixed point of $\tilde{g}$, an elementary computation yields that $v$ solves (4). Conversely, if $v$ solves (4) then Taylor’s formula with integral reminder shows that $u$ is a fixed point of $\tilde{g}$. ∎

It is worth noting that, in the same framework of rigorous computation as the one used here, approximations using piecewise linear functions were used in [18] to prove existence of homoclinic orbits for the Gray-Scott equation. In that case the system of ODEs considered is of order $2$, and therefore the equivalent integral operator is very similar to $\tilde{g}$ in (5) for $p=2$. Similarly, piecewise linear functions were used in [19] to prove existence of connecting orbits in the Lorenz equations. In that case, the standard Picard operator (3) was used.

Now that we have an operator which provides a gain of several orders of regularity, it becomes interesting to consider polynomial interpolation of higher order, and again this will be detailed in Section 5 and applied in Section 6.

2.2 Piecewise reformulation of the Picard-like operator

We finish this section by a last equivalent formulation of the initial value problem (2), that will be the one used in the present paper to perform the computer-assisted proofs. Given $m\in\mathbb{N}$, we introduce the mesh of $[0,1]$

[TABLE]

where $t_{0}=0<t_{1}<\ldots<t_{m}=1$. Then we consider $\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ (respectively $\mathcal{C}^{k}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$) the space of piecewise continuous (respectively $\mathcal{C}^{k})$ functions having possible discontinuities only on the mesh $\Delta_{m}$. More precisely, we use the following definition.

Definition 2.5**.**

For $k\in\mathbb{N}$, we say that $u\in\mathcal{C}^{k}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ if $u_{|_{(t_{j},t_{j+1})}}\in\mathcal{C}^{k}((t_{j},t_{j+1}),\mathbb{R}^{n})$ and can be extended to a $\mathcal{C}^{k}$ function on $[t_{j},t_{j+1}]$, for all $j=0,\ldots,m-1$.

We then introduce

[TABLE]

defined on the interval $(t_{j},t_{j+1})$ ($j=0,\dots,m-1$) by

[TABLE]

where $u(t_{j}^{-})$ denotes the left limit of $u$ at $t_{j}$, and $u(t_{0}^{-})$ must be replaced by $u_{0}$ (this last convention will be used throughout the paper).

Remark 2.6**.**

We point out that our computer-assisted proof is based on the operator $g$ (defined in (6)), which differs slightly from the operator $\tilde{g}$ (defined in (5)), which was used in previous studies such as [17, 18, 19]. The only difference is that the integral in $g$ is in some sense reseted at each $t_{j}$. We introduce this piecewise reformulation because it allows for sharper estimates (see Remark 4.7). **

We finally introduce $G:\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})\to\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ as

[TABLE]

Lemma 2.7**.**

Let $u\in\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$. Then $G(u)=0$ if and only if $v:t\mapsto u(\frac{t}{\tau})$ solves (2).

Proof.

This result is similar to Lemma 2.4. The only additional property that we need to check is that, if $u\in\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ satisfies $G(u)=0$, then $u$ cannot be discontinuous. Indeed, if $G(u)=0$ then $g(u)=u$, and for all $j\in\{1,\ldots,m-1\}$ one has

[TABLE]

At this point, it might seems as if defining $G$ on $\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ brings unnecessary complications, and that we should simply define it on $\mathcal{C}^{0}([0,1],\mathbb{R}^{n})$. While this is indeed a possibility, it will quickly become apparent that the present choice is more convenient, both for theoretical and numerical considerations (see Remark 3.4).

Finding a zero of $G$ is the formulation of our initial problem (2) that we are going to use in the rest of this paper.

3 General framework for the polynomial interpolation

3.1 Preliminaries

Given a mesh $\Delta_{m}$ as defined in Section 2.2, we introduce the refined mesh $\Delta_{m,k}$ where, for all $j\in\{0,\ldots,m-1\}$ we add $k-1$ points between $t_{j}$ and $t_{j+1}$. More precisely we suppose that these points are the Chebyshev points (of the second kind) between $t_{j}$ and $t_{j+1}$, that is we add the following points:

[TABLE]

where

[TABLE]

Notice that the above definition extends to $t_{j,0}=t_{j}$ and $t_{j,k}=t_{j+1}$, and that $k=1$ corresponds to the mesh used in previous studies with first order interpolation (e.g. see [18, 19, 17]).

We then introduce the subspace $S_{m,k}^{n}\subset\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ of piecewise polynomial functions of degree $k$ on $\Delta_{m}$

[TABLE]

Next, we define the projection operator

[TABLE]

where $\bar{u}$ is the function in $S_{m,k}^{n}$ that matches the values of $u$ on the mesh $\Delta_{m,k}$. Notice that $u$ can have discontinuities at the points $t_{j}$, therefore the matching of $u$ and $\bar{u}$ at those points must be understood as

[TABLE]

In the sequel we will need to control the error between a function $u$ and its interpolation $\bar{u}$. This is the content of the following propositions, where $\left\|\cdot\right\|_{\infty}$ denotes the sup norm on $[0,1]$.

Proposition 3.1**.**

For all $u\in\mathcal{C}^{k+1}_{\Delta_{m}}([0,1],\mathbb{R})$,

[TABLE]

where

[TABLE]

Proposition 3.2**.**

Fix $l\in\mathbb{N}$ such that $1\leq l\leq k$. For all $u\in\mathcal{C}^{l}_{\Delta_{m}}([0,1],\mathbb{R})$,

[TABLE]

where

[TABLE]

$\Lambda_{k}$* being the Lebesgue constant (see for instance [34]), and $\left[\frac{l-1}{2}\right]$ denoting the integer part of $\frac{l-1}{2}$.*

Remark 3.3**.**

More information about the Lebesgue constant, and in particular sharp computable upper bounds for it, can be found in the Appendix, together with references and proofs of the two above propositions. **

3.2 Finite dimensional projection

To get an approximate zero of $G$ (and thus an approximate solution of (2)), we are going to look for a zero of $\bar{G}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\Pi^{n}_{m,k}G_{|S^{n}_{m,k}}$. But first, we need a convenient way to represent the elements of $S^{n}_{m,k}$. Here and in the sequel, we use the exponent (i) to denote the $i$-th component of a vector in $\mathbb{R}^{n}$, but we will work with all the components at once as often as possible to avoid burdening the notations with this exponent (i). Let us introduce the set of indexes

[TABLE]

Perhaps the most natural way to characterize an element $\bar{u}$ of $S^{n}_{m,k}$ is to give all the values $\bar{u}^{(i)}(t_{j,l})$ for $(i,j,l)\in\mathcal{E}^{n}_{m,k}$. However, we will also use another representation, more suited to numerical computations, which consists of decomposing $\bar{u}$ on the Chebyshev basis. That is, we write

[TABLE]

where $T_{l}$ is the $l$-th Chebyshev polynomial of the first kind. We can thus also describe uniquely any function $u$ belonging to $S^{n}_{m,k}$ by the family of Chebyshev coefficients $\left(\bar{u}^{(i)}_{j,l}\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$.

Remark 3.4**.**

Let us mention how considering functions with possible discontinuities on the mesh points in $\Delta_{m}$ comes in handy. By restricting ourselves to functions in $\mathcal{C}^{0}([0,1],\mathbb{R}^{n})$, we would need additional constraints on the Chebyshev coefficients to impose the continuity at each of the mesh point $t_{j}$ ($j=1,\dots,m-1$) and keep track of them in all computations. Instead, the choice of working with $\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ allows avoiding these additional constraints. **

For $u\in\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ we have

[TABLE]

We recall that all the values $G^{(i)}(u)\left(t_{j,l}\right)$ for $(i,j,l)\in\mathcal{E}^{n}_{k,m}$ uniquely characterize $\Pi^{n}_{m,k}G(u)$.

Using the isomorphisms $\bar{u}\mapsto\left(\bar{u}^{(i)}(t_{j,l})\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$ and $\bar{u}\mapsto\left(\bar{u}^{(i)}_{j,l}\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$ to identify $S^{n}_{m,k}$ and $\mathbb{R}^{nm(k+1)}$, we can in fact see $\bar{G}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\Pi^{n}_{m,k}G_{|S^{n}_{m,k}}$ as a function from $\mathbb{R}^{nm(k+1)}$ to itself, that associates to the coefficients $\left(\bar{u}^{(i)}_{j,l}\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$ the values $\left(G^{(i)}(\bar{u})(t_{j,l})\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$. Thus we can numerically find a zero $\bar{u}$ of $\bar{G}$, which is going to be our approximate solution. We note that we use these identifications between $S^{n}_{m,k}$ and $\mathbb{R}^{nm(k+1)}$ throughout the present work. Our objective is now to validate this numerical solution $\bar{u}$, that is to prove that within a given neighbourhood of $\bar{u}$ lies a true zero $u$ of $G$.

3.3 Back to a fixed point formulation

We consider the space $\mathcal{X}^{n}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\mathcal{C}^{0}_{\Delta_{m}}([0,1],\mathbb{R}^{n})$ and its decomposition $\mathcal{X}^{n}=\mathcal{X}^{n}_{m,k}\oplus\mathcal{X}^{n}_{\infty}$, where

[TABLE]

We already have a projection onto $\mathcal{X}^{n}_{m,k}$

[TABLE]

and we also define its complementary

[TABLE]

We then define the norms

[TABLE]

On $\mathcal{X}^{n}$ we consider the norm

[TABLE]

where $r_{\infty}$ is a positive parameter. Notice that for all $r_{\infty}>0$, $\left(\mathcal{X}^{n},\|\cdot\|_{\mathcal{X}^{n}}\right)$ is a Banach space. For any $r,r_{\infty}>0$, we denote by $B_{\mathcal{X}^{n}}(r,r_{\infty})$ the closed neighbourhood of [math] defined as

[TABLE]

Suppose that we now have computed a numerical zero $\bar{u}$ of $\bar{G}$. We define $A_{m,k}^{{\dagger}}=D\bar{G}\left(\bar{u}\right)$ and consider $A_{m,k}$ an injective numerical inverse of $A_{m,k}^{{\dagger}}$. Finally, we introduce the Newton-like operator $T:\mathcal{X}^{n}\to\mathcal{X}^{n}$ defined by

[TABLE]

Notice that the fixed points of $T$ are in one-to-one correspondence with the zeros of $G$. We now give a finite set of sufficient conditions, that can be rigorously checked on a computer using interval arithmetic, to ensure that $T$ is a contraction on a given ball around $\bar{u}$. If those conditions are satisfied, the Banach fixed point theorem then yields the existence and local uniqueness of a zero of $G$. This is the content of the following statement (based on [35], see also [10] for a detailed proof).

Theorem 3.5**.**

Let

[TABLE]

and

[TABLE]

Assume that we have bounds $Y$ and $Z(r,r_{\infty})$ satisfying

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

If there exist $r,r_{\infty}>0$ such that

[TABLE]

and

[TABLE]

then there exists a unique zero of $G$ within the set $\bar{u}+B_{\mathcal{X}^{n}}(r,r_{\infty})\subset\mathcal{X}^{n}$.

The quantities $p^{(i)}_{j,l}(r,r_{\infty})$ and $p^{(i)}_{\infty}(r,r_{\infty})$ given respectively in (15) and (16) are called the radii polynomials.

In the next section, we show how to obtain bounds $Y$ and $Z$ satisfying (11)-(14). Before doing so, let us make a quick remark about the different representations and norms we can use on $\mathcal{X}^{n}_{m,k}$.

Remark 3.6**.**

As explained in Section 3.2, in practice we will mostly work with $\bar{u}\in\mathcal{X}^{n}_{m,k}$ represented by its Chebyshev coefficients as in (7). However, there are going to be instances where the values $\bar{u}^{(i)}(t_{j,l})$ are needed, for instance to compute $\left\|\bar{u}\right\|_{\mathcal{X}^{n}_{m,k}}$. We point out that numerically, passing from one representation to the other can be done easily by using the Fast Fourier Transform.

One other important point is that, at some point in the next section we are going to need upper bounds for $\left\|\bar{u}^{(i)}\right\|_{\infty}$. To get such a bound from our finite dimensional data, we have two options, namely

[TABLE]

or

[TABLE]

If $\bar{u}$ is given, then (17) is usually better, whereas (18) is better if $\bar{u}$ is any function in a given ball of $\mathcal{X}^{n}_{m,k}$. Notice that (17) simply follows from the fact that the Chebyshev polynomials satisfy $|T_{l}(t)|\leq 1$ for all $t\in[-1,1]$ and all $l\in\mathbb{N}$. For more information about the bound (18), see the Appendix and the references therein. **

4 Formula for the bounds

In this section, we give formulas for $Y^{(i)}_{j,l}$, $Y^{(i)}_{\infty}$, $Z^{(i)}_{j,l}$ and $Z^{(i)}_{\infty}$ satisfying the assumptions (11)-(14) of Theorem 3.5. To make the exposition clearer, we focus strictly on the derivation of the different bounds in this section. In particular, the discussion about the impact of the level of an priori bootstrap (that is the value of $p$) and the order of polynomial approximation (that is the value of $k$) is done in Section 5.

4.1 The $Y$ bounds

In this section we derive the $Y$ bounds, which measure the defect associated with a numerical solution $\bar{u}$, that is how close $G(\bar{u})$ is to [math]. We start by the finite dimensional part.

Proposition 4.1**.**

Let $y$ be defined as in (9) and consider

[TABLE]

where $\bar{G}(\bar{u})$ is here seen as the vector $\left(G^{(i)}(\bar{u})(t_{j,l})\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$. Then (11) holds.

Proof.

Simply notice that $\Pi^{n}_{m,k}y=-A_{m,k}\Pi^{n}_{m,k}G(\bar{u})$. ∎

Remark 4.2**.**

The above bound is not completely satisfactory, in the sense that is not directly implementable. Indeed, to compute $Y^{(i)}_{j,l}$ we need to evaluate (or at least to bound) the quantities $G^{(i)}(\bar{u})(t_{j,l})$. In particular (see (3.2)), we need to evaluate the integrals

[TABLE]

where

[TABLE]

If $\phi$ is a non polynomial vector field, we use a Taylor approximation of order $k_{0}$ of $\Psi$ to get an approximate value of the integral by computing

[TABLE]

Notice that this quantity can be evaluated explicitly. The error made in this approximation is then controlled as follows

[TABLE]

Notice that this error term is effective, since $\max_{s\in[-1,1]}\left\|\frac{d^{k_{0}+1}\Psi}{ds^{k_{0}+1}}(s)\right\|$ can be bounded using interval arithmetic. Therefore, the quantity $Y^{(i)}_{j,l}$ that we end up implementing is of the form

[TABLE]

where the vector $\hat{G}(\bar{u})$ contains the approximate integrals and the vector $G_{\epsilon}(\bar{u})$ contains the errors bounds for these approximations. Here and in the sequel, absolute values applied to a matrix, like $|A_{m,k}|$, must be understood component-wise. We point out that, in practice, if the mesh $\Delta_{m}$ is refined enough, then $\bar{u}$ is not going to be varying much on each subinterval $[t_{j},t_{j+1}]$, and thus we can get rather precise approximations even with a lower order $k_{0}$ for the Taylor expansion.

We mention that when the vector field $\phi$ is polynomial, $\Psi$ has a finite Taylor expansion, therefore up the integrals can in fact be computed exactly (i.e. we can get $G_{\epsilon}(\bar{u})=0$). **

We now turn our attention to the second part of the $Y$ bound.

Proposition 4.3**.**

Let $y$ be defined as in (9) and consider

[TABLE]

Then (12) holds.

Proof.

We have $\Pi^{n}_{\infty}y=\Pi^{n}_{\infty}(G(\bar{u})+\bar{u})=\Pi^{n}_{\infty}g(\bar{u})$. Since $\frac{d^{p}}{dt^{p}}g(\bar{u})=\tau^{p}\phi^{[p]}(\bar{u})$, we have that $\frac{d^{k+1}}{dt^{k+1}}g(\bar{u})=\tau^{p}\frac{d^{k+1-p}}{dt^{k+1-p}}\phi^{[p]}(\bar{u})$ and Proposition 3.1 yields

[TABLE]

Remark 4.4**.**

As comment similar to the one of Remark 4.2 applies here. Indeed, the bound given in Proposition 4.3 is not directly implementable because of the term

[TABLE]

but we can again get an explicit bound for this quantity by using a low order Taylor approximation and interval arithmetic. In the particular case where the vector field $\phi$ is polynomial, an explicit bound can also be obtained via the Chebyshev coefficients of the polynomial $\frac{d^{k+1-p}}{dt^{k+1-p}}(\phi^{[p]})^{(i)}(\bar{u})$, as in (17). **

4.2 The $Z$ bounds

In this section we derive the $Z$ bounds, which measure the contraction rate of $T$ on the ball of radius $r$ around $\bar{u}$. We begin with the finite dimensional part, that is the projection on $\mathcal{X}^{n}_{m,k}$. Let $z$ be defined as in (10). Then

[TABLE]

where $A_{m,k}$ and $A_{m,k}^{{\dagger}}$ are defined as in Section 3.3. We are going to bound each term separately as

[TABLE]

4.2.1 The bound $Z_{0}(r)$

The computation of the bounds $\left(Z_{0}(r)\right)^{(i)}_{j,l}$ estimating the first of the terms in the splitting (19) is rather straightforward and is simply a control on the precision of the numerical inverse.

Proposition 4.5**.**

Let $u_{2}\in B_{\mathcal{X}^{n}}(r,r_{\infty})$, define the vector $\mathbf{1}_{\mathcal{X}^{n}_{m,k}}\in\mathbb{R}^{nm(k+1)}$ by $\left(\mathbf{1}_{\mathcal{X}^{n}_{m,k}}\right)^{(i)}_{j,l}=1$ for all $(i,j,l)\in\mathcal{E}^{n}_{m,k}$ and let

[TABLE]

Then,

[TABLE]

4.2.2 The bound $Z_{1}(r,r_{\infty})$

We now construct the bounds $\left(Z_{1}(r,r_{\infty})\right)^{(i)}_{j,l}$ estimating the second term in the splitting (19).

Proposition 4.6**.**

Let $u_{2}\in B_{\mathcal{X}^{n}}(r,r_{\infty})$, consider $\rho=\left(\rho^{(i)}_{j,l}\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$ such that

[TABLE]

where $\mathbf{1}_{n}$ is the vector of size $n$ whose components all are equal to $1$. Let

[TABLE]

Then

[TABLE]

Proof.

By definition of $A_{m,k}^{{\dagger}}$ and $\bar{G}$, we have that

[TABLE]

Therefore, we can rewrite

[TABLE]

and we only need to prove that

[TABLE]

Remembering (3.2) and using that $\|\Pi^{n}_{\infty}u_{2}\|_{\mathcal{X}^{n}_{\infty}}\leq r_{\infty}r$, we estimate for all $(i,j,l)\in\mathcal{E}^{n}_{m,k}$,

[TABLE]

∎

Notice that Remark 4.4 also applies here.

Remark 4.7**.**

Had we used the operator $\tilde{g}$ (see (5)) instead of $g$ (see (6)), we would have gotten a bound like

[TABLE]

which is obviously worst because one has to consider the whole integral from [math] to $t_{j,l}$ instead of just from $t_{j}$ to $t_{j,l}$. **

4.2.3 The bound $Z_{2}$

We finally construct the bounds $\left(Z_{2}(r,r_{\infty})\right)^{(i)}_{j,l}$ estimating the last term in the splitting (19).

Proposition 4.8**.**

Let $u_{1},u_{2}\in B_{\mathcal{X}^{n}}(r,r_{\infty})$. consider $\\ varrho=\left(\\ varrho^{(i)}_{j,l}\right)_{(i,j,l)\in\mathcal{E}^{n}_{m,k}}$ such that

[TABLE]

Let

[TABLE]

Then

[TABLE]

Remark 4.9**.**

In the above proposition, $\left|D^{2+\delta}(\phi^{[p]})^{(i)}(\bar{u}(s))\right|\left(\mathbf{1}_{n}^{2+\delta}\right)$ must be understood as the evaluation of the $(2+\delta)$-linear form $\left|D^{2+\delta}(\phi^{[p]})^{(i)}(\bar{u}(s))\right|$ at the vectors $(\mathbf{1}_{n},\ldots,\mathbf{1}_{n})$, that is

[TABLE]

Proof.

(of Proposition 4.8) We only have to prove that

[TABLE]

Using (18) we have that

[TABLE]

Then we estimate for all $(i,j,l)\in\mathcal{E}^{n}_{m,k}$,

[TABLE]

Notice that Remark 4.4 also applies here.

4.2.4 The $Z_{\infty}$ bound

We are left with the remainder part of the $Z$ bound, which we treat in this section.

Proposition 4.10**.**

Let $u_{1},u_{2}\in B_{\mathcal{X}^{n}}(r,r_{\infty})$ and $z$ as in (10). Define for all $i\in\{1,\ldots,n\}$

[TABLE]

where $C^{opt}_{k,p}$ is one of the two constants given by Propositions 3.1 and 3.2, namely

[TABLE]

Then (14) holds.

Proof.

We need to estimate

[TABLE]

For any continuous function $\gamma$, one has

[TABLE]

thus we get, for all $1\leq i\leq n$

[TABLE]

Notice that Remark 4.4 also applies here.

4.3 The radii polynomials and interval arithmetics

The following proposition sums up what has been proven up to now in this section, namely that we have derived bounds that satisfy the requirements (11) to (14) from Theorem 3.5.

Proposition 4.11**.**

Let $y$ and $z$ defined as in (9) and (10). Then, the bound defined in Proposition 4.1 satisfies (11) and the one from Proposition 4.3 satisfies (12). Also, consider the bounds defined in Propositions 4.5 to 4.8. Then

[TABLE]

satisfies (13) and finally the bound from Proposition 4.10 satisfies (14).

Notice that, the way these bounds are defined, they are polynomials in $r$ and $r_{\infty}$, whose coefficients are all positive and can be computed explicitly with the help of the computer, since they depend on the numerical data of an approximate solution $\bar{u}$. Also, we make sure to control possible round-off errors by using interval arithmetic (in our case INTLAB [27]).

In practice, we first consider $r_{\infty}$ so that it satisfies the constraint (20) introduced in the next section. If there does not exist such positive $r_{\infty}$, we increase $m$ and/or $k$ and/or $p$ and try again. Once $r_{\infty}$ is fixed, we try to find a positive $r$ such that the last conditions (15) and (16) of Theorem 3.5 hold. If there is no such positive $r$, we increase $m$ and/or $k$ and/or $p$ and try again. If we finally find a positive $r$ satisfying (15) and (16), then we have proven that Theorem 3.5 applies, that is there exists a unique zero of $G$ in $B_{\mathcal{X}^{n}}(r,r_{\infty})$.

In Sections 6 and 7, we give several examples where the procedure described just above is successfully used to validate solutions of an initial value problem, as well as periodic solutions and heteroclinic orbits. But before doing so, we discuss in the next section the role of the parameters $k$, $m$ and $p$, and how they influence the bounds.

5 About the choice of the parameters

In this section, we explain how the parameters $k$, $m$ and $p$ should be chosen, and in particular we highlight how the a priori bootstrap (that is taking $p\geq 2$) helps improving the efficiency of the computer-assisted procedure we propose. The discussion will be rather informal, but we hope it helps the reader understand the results of the various comparisons presented in Section 6. Also, to make things slightly simpler we assume here that the grid $\Delta_{m}$ is uniform, therefore in the estimates each instance of $t_{j+1}-t_{j}$ can be replaced by $\frac{1}{m}$.

Our main constraint is that we want the method to be successful while minimizing the size of our numerical data, that is the dimension of our finite dimensional space $\mathcal{X}^{n}_{m,k}$, which is $nm(k+1)$. Since $n$ is fixed by the dimension of the vector field $\phi$, we want to minimize the product $m(k+1)$. As we see in the examples of Section 6, the usual limiting factor when trying to satisfy the radii polynomial inequalities (15) and (16) is to get the order one term (in $r$) to be negative. For the finite part (that is (15)), that means basically having that

[TABLE]

(see Proposition 4.6), and for the remainder part (that is (16)) we get a condition like

[TABLE]

where $\alpha$ and $\beta$ are constants depending on the numerical solution $\bar{u}$ and on the vector field $\phi$, but not on the parameters $k$, $m$ and $p$ that we can tune. This leads to

[TABLE]

We want to be able to chose a $r_{\infty}$ satisfying the above inequalities, and a necessary and sufficient condition for that is

[TABLE]

which we can rewrite

[TABLE]

Remember that we want (21) to be satisfied, while minimizing the product $m(k+1)$. When $p$ is fixed, and $k$ becomes large, notice that $C^{opt}_{k,p}$ is decreasing like $\frac{\ln(k)}{k^{p}}$. However, satisfying (21) requires, roughly speaking, to decrease $\left(\frac{\tau}{m}\right)^{p}C^{opt}_{k,p}$ as much as possible. This suggests two things, which we confirm in our explicit examples of Section 6. First, that it is slightly better to increase $m$ than $k$ (because of the $ln(k)$ factor) and second, that if we take $p$ equal to 2 or more (that is if we use a priori bootstrap) then we can satisfy (21) while taking $m(k+1)$ much smaller than if we had $p$ equal to 1.

Finally, we point out that taking $k=p-1$ seems optimal for the conditon (21) given by the order one term. Indeed, increasing $k$ from $p-1$ to $p$ increases the total number of coefficients, but brings no gain with respect to (21) since

[TABLE]

However, for the proof to succeed (that is for (15) and (16) to be satisfied) we also need small enough $Y$ and $Y_{\infty}$ bounds. Looking more precisely at $Y_{\infty}$, we see that it is of the form

[TABLE]

where $\gamma$ depends on the numerical data $\bar{u}$ and also on $k$, but the dependency on $k$ is way less important than in the $C_{k}\frac{1}{m^{k+1}}$ term, so we neglect it here. Looking back to the definition of $C_{k}$ in Proposition 3.1, we see that the term that we want to be small is of the form

[TABLE]

Therefore, if we really need to decrease the $Y_{\infty}$ bound, increasing $k$ is drastically better than increasing $m$. That is why, in practice we often take $k=p$, even though $k=p-1$ would be enough to satisfy (21). Finally we point out that, if we are not simply focused on getting an existence result, but also care about having sharp error bound, then we should definitively take care of having small $Y$ and $Y_{\infty}$ bounds, which, as we will show in the next section, can be achieved by slightly increasing $k$ (that is taking $k>p$).

In the next section, we present several comparisons for different choices of parameters, that confirm the heuristic presented in this section.

6 Examples of applications for the Lorenz system

In this section, we consider the Lorenz system, that is

[TABLE]

with standard parameter values $(\sigma,\beta,\rho)=(10,\frac{8}{3},28)$. Here, we first consider the initial value problem  (2), and use those bounds to try and validate orbits of various length with different parameters, to highlight the significant improvement made possible by the a priori bootstrap technique (that is taking $p\geq 2$). Then, we show that the a priori bootstrap also allows to validate more interesting solutions (from a dynamical point of view), namely periodic orbits and connecting orbits.

6.1 Comparisons for the initial value problem

The aim of this section is to showcase the improvements allowed by the use of a priori bootstrap, and to validate the heuristics made in Section 5. To do so, we fix an initial data (chosen close to the attractor of the Lorenz system)

[TABLE]

and do two kinds of comparisons. First, we try to validate the longest possible orbits for $p=1,2,3$ at various values of $m$ and $k$. We recall that by validating, we mean getting the existence of a true solution near a numerical one, by checking that the hypotheses of Theorem 3.5 hold. To make the comparison fair, we fix the total number of coefficients used for the numerical approximation, that is the dimension of $\mathcal{X}^{n}_{m,k}$, given by $nm(k+1)$. This quantity is usually the bottleneck of our approach, since we need to store and invert the matrix $A^{{\dagger}}_{m,k}$ which is of size $nm(k+1)\times nm(k+1)$. Here, we take $nm(k+1)=14000$ (or as close as possible to $14000$). The computations were made on a laptop with 8GB of RAM, and of course $nm(k+1)$ could be taken larger on a computer with more memory.

The first set of results are given in Table 1 (we recall that we work with the Lorenz system, therefore $n=3$).

In all cases, the proof fails for longer time $\tau$, because (21) is no longer satisfied. We see here that, as announced in Section 5, it is better to take $k$ as small as possible to get the longest possible orbit, but that increasing $k$ helps reducing the $Y_{\infty}$ bound, and thus the validation radius $r$. We see that simply increasing the order of the polynomial interpolation (given by $k$), allows to get better accuracy but does not really help to prove longer orbits. However, we are going to show on the next examples (see Table 2) that combining a priori bootstrap (that is taking $p\geq 2$) with higher order polynomial interpolation does allow to get much longer orbits.

First, comparing the $k=1$ case when $p=1$ and $p=2$, we see that using a priori bootstrap allows to get a slightly longer orbit, even for linear interpolation. Also, even for the longest possible orbit in that case ($\tau=0.97$), we still have much room to satisfy (21) (the quantity given by (21) is $\ll 1$). However, we cannot get a longer orbit in that case even with $p=2$, because the $Y_{\infty}$ bound becomes too large. This can be dealt with by increasing $k$, and we see that we can then get much longer orbits. To finish this set of comparisons, we show that doing one more iteration of the a priori bootstrap process (that is taking $p=3$ instead of $p=2$) still improves the results and allows to get longer orbits (see Table 3).

We sum up this set of comparisons by displaying the longest orbit obtained with $p=1$, $p=2$ and $p=3$ (see Figure 1).

We then finish this section with another set of comparisons, where we now fix the length of the orbit, here $\tau=2$, and instead look for the minimal total number of coefficients for which we can validate this orbit (for different values of $p$). The aim of this experiment is to show that using a priori bootstrap enables to validate solutions that one would not be able to validate without using it. Indeed, we are going to see that taking $p$ greater than one allows us to use way less coefficients to validate the solutions. Thus, if for a given solution, the proof without a priori bootstrap requires more coefficients than what our computer can handle, one can reduce this number by using a priori bootstrap and then possibly validate the orbit. For instance, still with the initial condition given by (22), we cannot validate the orbit of length $\tau=2$ without a priori bootstrap (that is with $p=1$), at least not with less that $14000$ coefficients. However, the next table of results shows that we can validate it with $p=2$, and also using even less coefficients with $p=3$.

6.2 Validation of a periodic orbit.

To study periodic orbits, instead of an initial value problem, the system (1) has to be slightly modified into a boundary value problem

[TABLE]

where $\tau$ is now an unknown of the problem, and where $u_{0},v_{0}\in\mathbb{R}^{n}$. The last equation is sometimes called Poincaré phase condition and is here to isolate to periodic orbit.

As for the initial value problem, we can then consider an equivalent integral formulation (possibly with a priori bootstrap) and define an equivalent fixed point operator $T$ very similar to the one introduced in Section 3. The additional phase condition and the fact that $\tau$ is now a variable only require minor modifications of $T$ and of the bounds derived in Section 4 (see for intance [12]).

Using a priori bootsrtap, we are able to validate fairly complicated periodic orbits (see Figure 2).

6.3 Validation of a connecting orbit

In this section, we present a computer-assisted proof of existence of a connecting orbit in the Lorenz system for the standard parameter values $(\sigma,\beta,\rho)=(10,\frac{8}{3},28)$. It is well know that at these parameter values, the Lorenz system admits a transverse connecting orbit between $\left(\sqrt{\beta(\rho-1)},\sqrt{\beta(\rho-1)},\rho-1\right)$ and the origin.

While computer-assisted proofs of connecting orbits were already investigated several times using topological and analytical approaches (see [18, 19, 25, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]), a paper of particular relevance to the present work is [19], where a particular case of our method is developed, with only linear interpolation and no a priori bootstrap (that is $k=1$ and $p=1$). While the authors in [19] were able to validate several connecting orbits for the Lorenz system, they could not validate the aforementioned connecting orbit for the standard parameter values. In fact, one of the main motivations for the present work was to improve the setting of [19] to be able to validate more complicated orbits.

As we showcased in Section 6.1, using a priori bootstrap enables us to validate significantly more complicated orbits for the initial value problem, and this is also true for connecting orbits. Indeed we are able to validate the standard connecting orbit for the Lorenz system, with parameter values $(\sigma,\beta,\rho)=(10,\frac{8}{3},28)$. Before exposing the results, we briefly describe how to modify (2) to be able to handle connecting orbits.

Compared to an initial value problem on a given time interval, or to a periodic orbit, connecting orbits present an aditionnal difficulty which is that they are defined on an infinite time interval (from $-\infty$ to $+\infty$). To circumvent this difficulty and get back to a time interval of finite length, which is more suited to numerical computations (and to computer-assisted proofs), we are going to use local stable and unstable manifolds of the fixed points. By a computer-assisted method very similar to the one presented here, we first compute and validate local parameterization of the unstable manifold at $\left(\sqrt{\beta(\rho-1)},\sqrt{\beta(\rho-1)},\rho-1\right)$ and of the stable manifold at the origin. Since the main object of this work is not the rigorous computations of those manifolds, we simply assume that they are given (with validation radius) and do not give more details here. The interested reader can find more information about the computations and validations of these parameterizations in [13, 48] and the references therein, and also more detailed examples of their usage to get connecting orbits in [18, 19, 11, 12].

We denote by $P$ a local parameterization of the stable manifold of the origin, and by $Q$ a local parameterization of the unstable manifold of a local parameterization of the stable manifold of $\left(\sqrt{\beta(\rho-1)},\sqrt{\beta(\rho-1)},\rho-1\right)$. We point out that both manifolds are two-dimensional. We then want to solve

[TABLE]

where $\varphi$ and $\theta$ each is a one dimensional parameter, the parameter in the other dimension being fixed to isolate the solution. Notice that $\tau$ is now an unknown of the system. As for the initial value problem, we can then consider an equivalent integral formulation (possibly with a priori bootstrap) and define an equivalent fixed point operator $T$ very similar to the one introduced in Section 3. The additional equation $u(\tau)=P(\theta)$ and the fact that we have three additional variables $\tau$, $\theta$ and $\varphi$ only requires minor modifications of $T$ and of the bounds derived in Section 4 (see for intance [12, 19]).

Using $p=3$, $k=3$ and $m=1150$ (that is a total number of 13803 coefficients) we are then able to rigorously compute a solution of (24) (see Figure 3).

7 Examples of applications for ABC flows

In this section, we apply our method to the non polynomial vector field

[TABLE]

The map $\phi_{A,B,C}$ is usually referred to as the Arnold-Beltrami-Childress (ABC) vector field, and gives a prime example of complex steady incompressible periodic flows in 3D (see [49, 50, 51] and the references therein).

The main point of this section is to briefly illustrate the applicability of our technique to non polynomial vector fields. We plan on studying more thoroughly ABC flows with the help of our a posteriori validation method in a future work. Recently, the existence of orbits, that are periodic up to a shift of $2\pi$ in one coordinates, have been proven in the cases $A=B=C=1$ and $0<A\ll 1,\ B=C=1$ [52, 53]. Applying the method developed in this paper, we were able to complete these results by proving the following statements.

Theorem 7.1**.**

For all $A=0.1,0.2,\ldots,1$, with $B=C=1$, there exists $\tau_{A}\in[\tau^{-}_{A},\tau^{+}_{A}]$ (see Table 5) and a solution $(x,y,z)$ of the ABC flow such that

[TABLE]

Proof.

The proof is done by running script_proofs_A11.m (available at [56]), which for each $A=0.1,0.2,\ldots,1$ computes an approximate solution, then computes bounds satisfying (11)-(14) as described in Section 4, and finally finds positive $r_{\infty}$ and $r$ such that (15)-(16) holds. ∎

Theorem 7.2**.**

For all $A=B=C=1$, there exists $\tau\in[7.797656,7.797666]$ and a solution $(x,y,z)$ of the ABC flow such that

[TABLE]

and $x(\cdot+\tau)\neq x(\cdot)+2\pi$.

Proof.

The proof is done by running script_proofs_111.m (available at [56]), which computes an approximate solution, then computes bounds satisfying (11)-(14) as described in Section 4, and finally finds positive $r_{\infty}$ and $r$ such that (15)-(16) holds. ∎

The solutions given by Theorem 7.1 are represented in Figure 4, and the solution given by Theorem 7.2 is represented in Figure 5.

Appendix

For the sake of completeness, we give here some properties of the Lebesgue constant $\Lambda_{k}$, as well as proofs of Proposition 3.1 and Proposition 3.2. We will assume here that $u$ is defined and smooth on $[-1,1]$. The corresponding estimates on $[t_{j},t_{j+1}]$ can the easily be deduced by rescaling.

We recall that $\Lambda_{k}$ is defined as the norm of the interpolation operator mapping $\mathcal{C}^{0}([-1,1],\mathbb{R})$ to itself, and associating a continuous function $u$ to its interpolation polynomial $P_{k}(u)$ of order $k$. Of course this operator (and its norm) depend on the interpolation points, which in this work are the Chebyshev interpolation points of the second kind

[TABLE]

Introducing the basis consisting of the Lagrange functions

[TABLE]

we have that the Lagrange interpolation polynomial of order $k$ is given by

[TABLE]

One can then show (see for instance [34]), that

[TABLE]

and therefore we get

[TABLE]

which is exactly (18).

Since we used several times (18) and Proposition 3.2 in Section 4, the bounds that we obtained there depend on the Lebesgue constant $\Lambda_{k}$. Therefore we need computable (and as sharp as possible) upper bounds for $\Lambda_{k}$. One possibility is to use the well known bound (again see for instance [34])

[TABLE]

However, we can do better, at least when $k$ is odd. In that case, it has been shown (see [54]) that

[TABLE]

and this formula can be evaluated rigorously using interval arithmetic. Unfortunately, there is no such formula for $k$ even. For small values ($k=2$ and $k=4$) we computed $\Lambda_{k}$ by hand using (26), and for $k\geq 6$ we used (27) (it is also know that $\Lambda_{k}\sim\frac{2}{\pi}\ln(k)$, therefore (27) is sharp for large $k$).

We now turn our attention to the interpolation estimates of Section 3. We point out that the analogue of Proposition 3.1 for the Chebyshev interpolation points of the first kind is very standard, and can be found in many textbooks. However, the case of the Chebyshev interpolation points of the second kind is seldom discussed, therefore we include a proof here for the sake of completeness (which is nothing but a slight adaptation of the standard proof).

Proof of Proposition 3.1. We consider the polynomial $W_{k}(x)=\prod_{l=0}^{k}(x-x_{l}^{k})$ and use the standard interpolation error estimate for a function $u\in\mathcal{C}^{k+1}$ (see for instance [55]),

[TABLE]

To prove Proposition 3.1, we only have to show that $\left\|W_{k}\right\|_{\infty}=\frac{1}{2^{k-1}}$ (the remaining factor $\frac{1}{2^{k+1}}$ coming from the rescaling). Introducing, for $k\in\mathbb{N}$, the $k$-th Chebyshev polynomial of the second kind $U_{k}$, defined by

[TABLE]

we have that

[TABLE]

Indeed, the right hand side of (28) is a unitary polynomial of degree $k+1$, that has the same zeros as $W_{k}$. We can then rewrite

[TABLE]

and we end up with

[TABLE]

so Proposition 3.1 is proven. ∎

Proof of Proposition 3.2. The first part of the bound, namely

[TABLE]

comes from a combination of the Lebesgue constant and Jackson’s Theorem, and can be found in [55]. However, it does not give a very sharp interpolation error estimate for small values of $k$ and $l$, therefore we derive here the second part of the bound, namely

[TABLE]

that can be used in those cases.

Letting $u(x)=x^{p}$ (with $p\in\{0,\ldots,k\}$) in (25) leads to

[TABLE]

We now fix a function $u\in\mathcal{C}^{l}$. Using (29) with $p=0$ (that is $1=\sum_{i=0}^{k}L^{k}_{i}(x)$), we get

[TABLE]

Using Taylor’s formula, we then get

[TABLE]

for some $y_{i}$ in $[-1,1]$. Then, expanding the $(x^{k}_{i}-x)^{p}$ terms and using again (29), we get that

[TABLE]

and thus

[TABLE]

Letting

[TABLE]

we can easily observe that $L_{i}^{k}(x)=\lambda^{k}_{i}W_{k}(x)/(x-x_{i}^{k})$, and therefore

[TABLE]

According to [34], in case the points $x^{k}_{i}$ are the Chebyshev interpolation points of the second kind, we have

[TABLE]

Remembering that $\left\|W_{k}\right\|_{\infty}=\frac{1}{2^{k-1}}$, we get

[TABLE]

The function

[TABLE]

is even and increasing on $[0,1]$, therefore its maximum is reached at $x=1$ and we get

[TABLE]

Then, we compute

[TABLE]

We end up with

[TABLE]

and Proposition 3.2 is proven (the lacking $\frac{1}{2^{l}}$ factor coming from the time rescaling). ∎