Validated computations for connecting orbits in polynomial vector fields

TL;DR

This paper introduces a computer-assisted method to rigorously prove the existence of connecting orbits in polynomial vector fields, combining Taylor approximations, Chebyshev series, and fixed-point analysis.

Contribution

It develops a novel, rigorous computational framework for establishing the existence and transversality of heteroclinic orbits in polynomial vector fields.

Findings

Successfully proves existence of heteroclinic orbits

Provides explicit bounds on approximation errors

Establishes transversality of stable and unstable manifolds

Abstract

In this paper we present a computer-assisted procedure for proving the existence of transverse heteroclinic orbits connecting hyperbolic equilibria of polynomial vector fields. The idea is to compute high-order Taylor approximations of local charts on the (un)stable manifolds by using the Parameterization Method and to use Chebyshev series to parameterize the orbit in between, which solves a boundary value problem. The existence of a heteroclinic orbit can then be established by setting up an appropriate fixed-point problem amenable to computer-assisted analysis. The fixed point problem simultaneously solves for the local (un)stable manifolds and the orbit which connects these. We obtain explicit rigorous control on the distance between the numerical approximation and the heteroclinic orbit. Transversality of the stable and unstable manifolds is also proven.

Tables

Table 1. Table 6.1: The (approximate) computational parameters used to validate connecting orbits for γ ∈ [ 0.4557 , 10.50 ] 𝛾 0.4557 10.50 \gamma\in\left[0.4557,10.50\right] . Each row in the table corresponds to an interval on which we performed (non-rigorous) pseudo-arclength continuation and validated rigorously a finite number of connecting orbits with the same computational parameters. In particular, we used K u = [ 15 15 ] superscript 𝐾 𝑢 matrix 15 15 K^{u}=\begin{bmatrix}15&15\end{bmatrix} on each interval (though validation with less Taylor coefficients is feasible). We were not able to validate connecting orbits past the right endpoints of the intervals without modifying the computational parameters. In each case, we have indicated the obstruction for validating connecting orbits near the (right) endpoint of the interval.

$\gamma$	$L$	$N$	$K^s$	${ϵ}_{u,k}$	${ϵ}_{s,k}$	${\nu }_i$	Obstruction
$[0.4557,0.7046]$	$30$	$\left[6261\right]$	$12$	$5.4\cdot {10}^{-2}$	$5.3\cdot {10}^{-3}$	$1.15$	${\hat{\lambda }}_1^u\approx 2{\hat{\lambda }}_2^u$
$[0.7146,1.233]$	$30$	$\left[6261\right]$	$9$	$5.4\cdot {10}^{-2}$	$5.3\cdot {10}^{-3}$	$1.15$	${\hat{\lambda }}_1^u\approx {\hat{\lambda }}_2^u$
$[1.243,4.089]$	$8$	$63$	$9$	$0.12$	$0.17$	$1.11$	${\hat{\lambda }}_1^s+{\hat{\lambda }}_2^s\approx {\hat{\lambda }}_3^s$
$[4.202,10.50]$	$4$	$51$	$7$	$0.17$	$0.12$	$1.15$	$\hat{r}>r^{\ast }$

Equations

$$\{\begin{array}{ll}\frac{du}{dt}=g\left(u\right),& t\in [0,L],\\ u(0)\in W_{\text{loc}}^u\left(p_0\right),& \\ u(L)\in W_{\text{loc}}^s\left(q_0\right),& \end{array}$$

$$\{\begin{array}{l}\frac{\partial v}{\partial t}=D\frac{{\partial }^{2}v}{\partial x^2}+v\left(1-v-w\right),\\ \frac{\partial w}{\partial t}=\frac{{\partial }^{2}w}{\partial x^2}+aw\left(v-b\right),\end{array}$$

$$\underset{\tau \to -\infty }{\lim }\left({\zeta }_1\left(\tau \right),{\zeta }_2\left(\tau \right)\right)=(1,0),\underset{\tau \to \infty }{\lim }\left({\zeta }_1\left(\tau \right),{\zeta }_2\left(\tau \right)\right)=(b,1-b).$$

$$\frac{du}{dt}=\left[\begin{array}{c}-u_2\\ D^{-1}\left(\kappa u_2+u_1\left(1-u_1-u_3\right)\right)\\ -u_4\\ \left(\kappa u_4+a{u}_3\left(u_1-b\right)\right)\end{array}\right]$$

$$\displaystyle\frac{\partial v}{\partial t}=-\gamma\frac{\partial^{4}v}{\partial x^{4}}+\frac{\partial^{2}v}{\partial x^{2}}+\bigl{(}v-a\bigr{)}\bigl{(}1-v^{2}\bigr{)}$$

$$\frac{du}{dt}=-\left[\begin{array}{c}\gamma u_2\\ \gamma u_3\\ \gamma u_4\\ \kappa u_2+u_3+\left(u_1-a\right)\left(1-u_1^2\right)\end{array}\right],$$

$$u=a_0+2\sum _{k=1}^{\infty }{a}_k{T}_k$$

$$u=a_0+2\sum _{k=1}^{\infty }{a}_k{T}_k,$$

$$u\cdot v={\left(a\ast b\right)}_0+2\sum _{k=1}^{\infty }{\left(a\ast b\right)}_k{T}_k,\text{where}{\left(a\ast b\right)}_k:={\sum _{k_1+k_2=k}}_{k_1,k_2\in \mathbb{Z}}{a}_{|k_1|}{b}_{|k_2|}.$$

$${\ell }_{\nu ,n}^1:=\left\{{\left(a_k\right)}_{k\in {\mathbb{N}}_0}:a_k\in {\mathbb{C}}^n,\sum _{k=0}^{\infty }\left|{\left[a_k\right]}_j\right|{\nu }^k<\infty ,1\le j\le n\right\},$$

$${\Vert a\Vert }_{\nu ,n}:=\underset{1\le j\le n}{\max }\left\{\left|{\left[a_0\right]}_j\right|+2\sum _{k=1}^{\infty }\left|{\left[a_k\right]}_j\right|{\nu }^k\right\}.$$

$$a=\sum _{k=0}^{\infty }{a}_k{e}_k,$$

$${\epsilon }_{k,\nu }=\{\begin{array}{ll}1& k=0,\\ \frac{1}{2}{\nu }^{-k},& k\in \mathbb{N}.\end{array}$$

$${\Vert \mathcal{L}\Vert }_{\mathcal{B}\left({\ell }_{\nu }^1,X\right)}=\underset{k\in {\mathbb{N}}_0}{\sup }{\Vert \mathcal{L}\left({\xi }_k\right)\Vert }_X.$$

$${\Vert \mathcal{L}\Vert }_{\mathcal{B}\left({\ell }_{\nu }^1,X\right)}\ge \underset{k\in {\mathbb{N}}_0}{\sup }{\Vert \mathcal{L}\left({\xi }_k\right)\Vert }_X,$$

$$a=a_0{\xi }_0+2\sum _{k=1}^{\infty }{a}_k{\xi }_k{\nu }^k.$$

$$\mathcal{L}\left(a\right)=a_0\mathcal{L}\left({\xi }_0\right)+2\sum _{k=1}^{\infty }{a}_k\mathcal{L}\left({\xi }_k\right){\nu }^k.$$

$${\Vert \mathcal{L}\left(a\right)\Vert }_X\le \underset{k\in {\mathbb{N}}_0}{\sup }{\Vert \mathcal{L}\left({\xi }_k\right)\Vert }_X{\Vert a\Vert }_{\nu }$$

$$\mathcal{L}\left(e_j\right)=\sum _{i=0}^{\infty }{\mathcal{L}}_{ij}{e}_i\simeq {\left[\begin{array}{cccc}{\mathcal{L}}_{0j}& {\mathcal{L}}_{1j}& \dots & \end{array}\right]}^T,j\in {\mathbb{N}}_0.$$

$$\mathcal{L}(a)=\left[\begin{array}{ccc}{\mathcal{L}}_{00}& {\mathcal{L}}_{01}& \dots \\ {\mathcal{L}}_{10}& {\mathcal{L}}_{11}& \dots \\ ⋮& & \end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\end{array}\right],$$

$${\Vert \mathcal{L}\Vert }_{\mathcal{B}\left({\ell }_{{\nu }_1}^1,{\ell }_{{\nu }_2}^1\right)}=\underset{l\in {\mathbb{N}}_0}{\sup }{\epsilon }_{l,{\nu }_1}{\Vert {\mathcal{L}}_{\left(\cdot ,l\right)}\Vert }_{{\nu }_2}.$$

$${\Vert \mathcal{L}\left(a\right)\Vert }_{{\nu }_2}$$

$$=\sum _{l=0}^{\infty }\left|a_l\right|{\Vert {\mathcal{L}}_{\left(\cdot ,l\right)}\Vert }_{{\nu }_2}$$

$$=\left|a_0\right|{\Vert {\mathcal{L}}_{\left(\cdot ,0\right)}\Vert }_{{\nu }_2}+2\sum _{l=1}^{\infty }{\epsilon }_{{\nu }_1,l}\left|a_l\right|{\Vert {\mathcal{L}}_{\left(\cdot ,l\right)}\Vert }_{{\nu }_2}{\nu }_1^l$$

$$\le \left(\underset{l\in {\mathbb{N}}_0}{\sup }{\epsilon }_{l,{\nu }_1}{\Vert {\mathcal{L}}_{\left(\cdot ,l\right)}\Vert }_{{\nu }_2}\right){\Vert a\Vert }_{{\nu }_1}.$$

$$p_k:=p_{k_1\dots k_d}:=p\left(k\right),k\in {\mathbb{N}}_0^d.$$

$$k⪯l{\iff }^{def}{k}_j\le l_j,\forall 1\le j\le d.$$

$$W_{\nu ,n,d}^1:=\left\{p:{\mathbb{N}}_0^d\to {\mathbb{C}}^n|\sum _{k\in {\mathbb{N}}_0^d}\left|{\left[p_k\right]}_j\right|{\nu }^{|k|}<\infty ,1\le j\le n\right\}$$

$${\Vert p\Vert }_{W_{\nu ,n,d}^1}:=\underset{1\le j\le n}{\max }\sum _{k\in {\mathbb{N}}_0^d}\left|{\left[p_k\right]}_j\right|{\nu }^{|k|}.$$

$$f(z)=\sum _{k\in {\mathbb{N}}_0^d}{f}_k{z}^k,g(z)=\sum _{k\in {\mathbb{N}}_0^d}{g}_k{z}^k,$$

$$(fg)(z)=\sum _{k\in {\mathbb{N}}_0^d}{\left(\tilde{f}\ast \tilde{g}\right)}_k{z}^k,{\left(\tilde{f}\ast \tilde{g}\right)}_k:=\sum _{\begin{array}{c}\alpha +\beta =k,\\ \alpha ,\beta \in {\mathbb{N}}_0^d\end{array}}{\tilde{f}}_{\alpha }{\tilde{g}}_{\beta },$$

$${\Vert \mathcal{L}\Vert }_{\mathcal{B}\left(W_{\nu }^1,X\right)}=\underset{k\in {\mathbb{N}}_0^d}{\sup }{\Vert \mathcal{L}\left({\mathbf{\xi }}_k\right)\Vert }_X.$$

Full text

1 Introduction

Connecting orbits play a central role in the study of dynamical systems. They provide a detailed picture of how a dynamical system can evolve from one “state” (e.g. an equilibrium, a periodic orbit or another type of recurrent set) into another. Furthermore, their existence can often be used to establish more complicated dynamical phenomena through forcing theorems. However, proving the existence of a connecting orbit for a given nonlinear ODE is in general a difficult (if not impossible) task to accomplish by hand. For this reason, one often resorts to numerical methods. While numerical methods can provide valuable insight into quantitive properties of a connecting orbit, which would otherwise be out of reach with merely a pen and paper analysis, the results are usually non-rigorous and cannot be used in mathematical arguments. In particular, a standard numerical method does not yield a proof for the existence of a connecting orbit.

In this paper we present a general computer-assisted method for proving the existence of transverse connecting orbits between hyperbolic equilibria for nonlinear ODEs. The method is based on solving the finite time boundary value problem

[TABLE]

where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a general polynomial vector field, $p_{0},q_{0}\in\mathbb{R}^{n}$ are hyperbolic equilibria and $L>0$ is the time needed to travel between the local (un)stable manifolds. We assume that $\dim\left(W^{u}\left(p_{0}\right)\right)+\dim\left(W^{s}\left(q_{0}\right)\right)=n+1$, which is a necessary condition for a transverse connecting orbit to exist. The idea is to solve (1.1) by computing Taylor expansions for charts on the local (un)stable manifolds via the parameterization method [ManifoldTaylor, MR3437754], which are used to supplant the boundary conditions in (1.1) with explicit equations, and to use Chebyshev series and domain decomposition techniques [domaindecomposition] to parameterize the orbit in between. We remark that the assumption that $g$ is polynomial is not as restrictive as it initially might seem, since many nonlinearities which consist of elementary functions can be brought into polynomial form by using automatic differentiation techniques, see [MR3545977] for instance.

Before we proceed with a more detailed description of our method, a few remarks concerning the development of numerical methods for connecting orbits are in order. Many numerical methods (both rigorous and non-rigorous) are based on approximating (un)stable manifolds and solving finite time boundary value problems. We mention the numerical (non-rigorous) methods implemented in the continuation packages Matcont [MR2000880] and AUTO [AUTO] in particular. Furthermore, many validated numerical methods have been developed over the last decade. It is beyond the scope of this text to give an overview. Nevertheless, we mention the functional analytic methods developed in [Maxime, domaindecomposition, MR3281845, MR3022075, suspensionbridge, MR3353132, MR3207723, MR3148084, MR2821596, Kepley, JayAMS], which are based on solving fixed point problems, the topological methods developed in [MR3463691, MR2173545, MR2505658, MR1961956], which are based on covering relations, cone conditions [MR2494688, MR2060531, MR2060532] and rigorous integration of the flow via Lohner-type algorithms [MR1930946], and the methods in [MR2302059, MR2339601] based on shadowing techniques.

The first step in the development of our validated numerical method is to recast (1.1) into an equivalent zero finding problem $F(x)=0$. The unknowns in this problem are the Taylor coefficients of the parameterizations of the local (un)stable manifolds, which include the equilibria and the associated eigendata, the coordinates of the endpoints $u(0)$ and $u(L)$ on the associated charts, and the Chebyshev coefficients of the orbit. Next, we use the computer to determine an approximate zero of $F$. The numerical computations are then combined with analysis on paper to construct a Newton-like map $T$ whose fixed points correspond to zeros of $F$. Finally, we use pen and paper estimates to derive a finite number of inequalities, which can be used to determine a neighborhood around the approximate solution on which $T$ is a contraction. An essential property of these inequalities is that they can be rigorously verified with the aid of a computer. This approach is in the literature often referred to as a parameterized Newton-Kantorovich method or the radii-polynomial approach (see [MR1639986, MR2338393]).

Our construction builds on the work in [suspensionbridge, MR3353132, MR3148084]. In those papers heteroclinic orbits for particular vector fields were studied using the parameterization method, Chebyshev series and the radii-polynomial approach to validate solutions of the boundary value problem (1.1). Whereas in [suspensionbridge, MR3353132, MR3148084] the charts on the local (un)stable manifolds were validated separately from the proof of the connecting orbit between them, the main contribution of the current paper is that we solve all ingredients of the problem simultaneously, and in great generality. While this means we have to introduce a fair amount of notation, the advantage is that our framework is very flexibly. Indeed, the mathematical analysis and the code are developed in full detail for general polynomial vector fields. Moreover, setting up this framework opens the door for various extensions. In particular, the formulation as a zero finding problem in an appropriately chosen (product) Banach space provides a foundation to study continuation and bifurcation problems, see also Section 1.2. Furthermore, the presented method deals in a unified and systematic manner with the cases of real and complex eigenvalues associated to the parameterizations of the (un)stable manifolds (see also the examples in Section 1.1). We note that in part of the paper we impose a technical “non-resonance” condition on the eigenvalues of the linearized problems at the equilibria, see Section 3.1. This condition is by no means fundamental, but it reduces the (already heavy) notational burden. As explained in Section 1.2, the general case follows by combining the current work with [ManifoldTaylor].

Finally, we note that it is a feature of the Newton-Kantorovich fixed point method that any heteroclinic orbit that we find using our methodology is the transverse intersection of the stable and unstable manifolds of the equilibria it connects, see Proposition 3.22. Such transversality information is often essential when one aims to use the connecting orbits as ingredients for forcing theorems.

In this paper we at times refer to previous work for certain proofs, bit we give sufficient details of the constructions to understand the algorithm. Hence the duo of paper and code is self-contained. The fully documented Matlab code is available at [heteroclinicscode].

1.1 Applications

To illustrate the efficacy of the method, we have chosen two travelling wave problems originating from partial differential equations Here we present some example results. Additional rigorously validated orbits and more details can be found in Section 6.

Application 1** (Traveling fronts in the Lotka-Volterra equations).**

The Lotka-Volterra equations are a system of reaction-diffusion equations given by

[TABLE]

where $D>0,a>0,b\in(0,1)$ and $\left(t,x\right)\in\mathbb{R}^{2}$. This system has three homogeneous equilibrium states: $(v,w)=(0,0)$, $(v,w)=(1,0)$ and $(v,w)=(b,1-b)$. We have used our method to prove the existence of solutions of (1.2) of the form $v(t,x)=\zeta_{1}(x-\kappa t)$ and $w(t,x)=\zeta_{2}(x-\kappa t)$, where $\zeta_{1},\zeta_{2}:\mathbb{R}\rightarrow\mathbb{R}$ and $\kappa<0$, which satisfy

[TABLE]

Such solutions are often referred to as traveling fronts with wave speed $\kappa$.

Substitution of the traveling wave Ansatz $\left(\zeta_{1},\zeta_{2}\right)$ into (1.2) shows that connecting orbits from $(b,0,1-b,0)$ to $(1,0,0,0)$ for the four dimensional system of ODEs

[TABLE]

correspond to traveling wave profiles $\left(\zeta_{1}(t),\zeta_{2}(t)\right)=\left(u_{1}(-t),u_{3}(-t)\right)$ and vice versa. We have successfully validated connecting orbits in (1.3) for various values of $\kappa\in\left[-1,-0.5938\right]$. For these parameter values, the equilibria $(b,0,1-b,0)$ and $(1,0,0,0)$ have a two dimensional unstable and three dimensional stable manifold, respectively. In particular, the stable eigenvalues of the linearization at $\left(1,0,0,0\right)$ consist of one complex conjugate pair of eigenvalues and one real eigenvalue. We have depicted a validated traveling wave profile and the corresponding connecting orbit for a particular wave speed in Figures 1.1 and 1.2, respectively. The reader is referred to Section 6.1 for the details.

Application 2** (Traveling fronts in a fourth order parabolic PDE).**

We have proven the existence of traveling fronts $v(t,x)=\zeta\left(x-\kappa t\right)$ for the following fourth order parabolic PDE:

[TABLE]

where $-1<a\leq 0$ and $\gamma>0$. Such travelling waves have been studied in [MR1629027] used geometric singular perturbation theory for small $\gamma$, and in [BHV] by Conley index techniques. In this paper we focus on traveling fronts between the homogeneous states $v\equiv-1$ and $v\equiv a$ for (relatively) large $\gamma$.

We have successfully established the existence of connecting orbits from $(-1,0,0,0)$ to $(a,0,0,0)$ for the four dimensional system of ODEs

[TABLE]

which correspond to traveling wave profiles $\zeta\left(t\right)=u_{1}\left(-\frac{t}{\gamma}\right)$, for fixed values of $a$ and the wave speed $\kappa$, and various values of $\gamma\in[0.4557,10.50]$. We rescaled time with a factor $\gamma$ so that the system in (1.5) is well-defined at $\gamma=0$. We have depicted a validated traveling wave profile and the corresponding connecting orbit in Figures 1.3 and 1.4, respectively. The reader is referred to Section 6.2 for the details.

1.2 Extensions and future work

We now discuss possible extensions. The first extension is to let the vector field explicitly depend on a parameter and to perform rigorous (pseudo-arclength) continuation of connecting orbits. This involves a relatively straightforward application of the uniform contraction principle and a slight modification of the estimates developed in this paper (see [Queirolo, MR2630003, MR3125637, suspensionbridge] for instance). Furthermore, in order to carry out continuation efficiently, we need to develop algorithms (heuristics) which automatically determine near-optimal parameter values for the validation of the charts on the local (un)stable manifolds and the connecting orbit. More specifically, during continuation it might become necessary to modify the number of Taylor coefficients, the size of the charts on the local (un)stable manifolds, the grid on which the connecting orbit is computed, the number of Chebyshev coefficients, or the integration time.

The second extension involves the incorporation of resonances. In this paper, we assume that the (un)stable eigenvalues associated to the (un)stable manifolds satisfy a so-called non-resonance condition. This condition is related to the regularity of the chart mappings obtained via the parameterization method. In short, the parameterization method is based on constructing a smooth conjugacy (analytic in our case) between the nonlinear flow on the (un)stable manifold and an “easier” fully understood model system. One can choose this model system to be linear, which we do in this paper, if the (un)stable eigenvalues satisfy a non-resonance condition. If there are resonant eigenvalues, however, one needs to use a nonlinear model system instead. This is explained in detail in [ManifoldTaylor]. Generically, one will encounter resonances during continuation. Therefore, in order to successfully perform validated continuation, we need to allow for the possibility of resonant eigenvalues and modify the current method accordingly as explained in [ManifoldTaylor]. In particular, we need to develop an algorithm which automatically detects when to “switch” between the linear and nonlinear model flow during continuation. Furthermore, a careful analysis of the case in which a pair of complex conjugate eigenvalues become real (or vice versa) is needed as well. After the above extensions have been implemented, one can start developing tools for the rigorous study of bifurcations of connecting orbits for nonlinear ODEs.

The computer-assisted method presented in this paper is implemented in an object oriented framework in Matlab using the Intlab package [Intlab] for interval arithmetic. A third and useful extension would be to incorporate an extra degree of freedom into the classes for the connecting orbit so that additional equations and variables can be added (or removed) in a convenient manner. This would facilitate the required modifications for dealing with non-polynomial vector fields via automatic differentiation techniques, analyzing connecting orbits in vector fields with symmetry, proving the existence of homoclinic instead of heteroclinic orbits, and performing bifurcation analysis.

Finally, we remark that the computational efficiency of the current implementation can be improved. For instance, the equations for the parameterizations of the local (un)stable manifolds and the connecting orbit in between are to a large extent uncoupled. As a consequence, the derivative of the zero finding map $F$ has a block structure, which can be exploited to reduce the computational costs of the computation of an approximate inverse (we need an explicit finite dimensional approximate inverse to construct a Newton-like map $T$). Furthermore, in applications it might not be necessary to resolve the full local stable manifold equally well in all directions, but a rather more focussed parameterization centered around the slow eigendirections of the equilibria is appropriate, since a connecting orbit generically tends to enter the stable manifold via these directions.

1.3 Outline of the paper

This paper is organized as follows. In Section 2 we review some basic facts about Chebyshev series, Taylor series and sequence spaces, which will be used extensively throughout this paper. In Section 3 we set up an equivalent zero finding problem for (1.1) by using domain decomposition, the parameterization method, Chebyshev series and Taylor series. In Section 4 we set up an equivalent fixed-point problem and explain how the existence of a zero can be established with the aid of a computer. This involves the construction of computable bounds which are developed in full detail in Section 5. Finally, in Section 6 we demonstrate the effectiveness of the method by proving the existence of traveling fronts in parabolic PDEs. We also discuss some algorithmic aspects.

2 Preliminaries

In this section we develop a functional analytic framework for analyzing maps which arise from the study of connecting orbits. We start in Section 2.1 by recalling basic results from Chebyshev approximation theory. In Sections 2.2 and 2.3 we introduce spaces of geometrically decaying sequences and multivariate arrays, respectively. In addition, we review methods for analyzing bounded linear operators on them.

2.1 Chebyshev series

In this section we recall basic notions and results from Chebyshev approximation theory. The reader is referred to [MR3012510] for the proofs and a more comprehensive introduction into the theory of Chebyshev approximations.

Definition 2.1**.**

The Chebyshev polynomials $T_{k}:[-1,1]\rightarrow\mathbb{R}$ are defined by the relation $T_{k}\left(\cos\left(\theta\right)\right)=\cos\left(k\theta\right)$, where $k\in\mathbb{N}_{0}$ and $\theta\in\left[0,\pi\right]$.

Chebyshev series constitute a non-periodic analog of Fourier cosine series and have similar convergence properties. For instance, any Lipschitz continuous function admits a unique Chebyshev expansion. In this paper, we will consider Chebyshev expansions of analytic functions. The Chebyshev coefficients of such regular functions decay (in analogy with Fourier series) at a geometric rate to zero. A more precise statement is given in the next proposition.

Proposition 2.2**.**

Suppose $u:\left[-1,1\right]\rightarrow\mathbb{R}$ is analytic and let

[TABLE]

be its Chebyshev expansion. Let $\mathcal{E}_{\nu}\subset\mathbb{C}$ denote an open ellipse with foci $\pm 1$ to which $u$ can be analytically extended, where $\nu>1$ is the sum of the semi-major and semi-minor axis of $\mathcal{E}_{\nu}$. If $u$ is bounded on $\mathcal{E}_{\nu}$, then $\left|a_{k}\right|\leq M\nu^{-k}$ for all $k\in\mathbb{N}_{0}$, where $M=\sup_{z\in\mathcal{E}_{\nu}}\left|u(z)\right|$.

The Chebyshev coefficients of the product of two Chebyshev series is (in direct analogy with Fourier cosine series) given by the symmetric discrete convolution:

Proposition 2.3**.**

Suppose $u,v:\left[-1,1\right]\rightarrow\mathbb{R}$ are Lipschitz continuous and let

[TABLE]

be the associated Chebyshev expansions. Then

[TABLE]

2.2 Geometrically decaying sequences

In this section we introduce a sequence space suitable for analyzing analytic functions, and elementary operations on them, via their Chebyshev coefficients. Recall that the Chebyshev coefficients of an analytic function $u:[a,b]\rightarrow\mathbb{R}^{n}$ decay exponentially fast to zero by Proposition 2.2. In light of this observation we define

[TABLE]

where $\left[a_{k}\right]_{j}$ denotes the $j$-th component of $a_{k}$ and $\nu>1$ is some prescribed weight, endowed with the norm

[TABLE]

In the special case that $n=1$ we shall write $\ell^{1}_{\nu}:=\ell^{1}_{\nu,1}$ and $\left\|\cdot\right\|_{\nu}:=\left\|\cdot\right\|_{\nu,1}$. It is a straightforward task to verify that $\ell^{1}_{\nu,n}$ equipped with this norm is a Banach space over $\mathbb{C}$.

Remark 2.4**.**

In this paper we are exclusively concerned with Chebyshev expansions of real-valued functions. From this perspective it is more natural to consider sequence spaces over $\mathbb{R}$ instead of $\mathbb{C}$ . The reason for using a space of complex valued sequences is that we wish to couple the Chebyshev expansions with chart maps for (un)stable manifolds, which might be complex-valued (see Section 3). We will proof a-posteriori that the Chebyshev coefficients are real by using arguments based on symmetry.

The operation of multiplying two Chebyshev series can be lifted to the level of sequences, giving rise to the symmetric discrete convolution $\ast$, as shown in Proposition 2.3. This additional product structure on $\ell^{1}_{\nu}$ yields a particularly nice space:

Proposition 2.5**.**

The space $\left(\ell^{1}_{\nu},\ast\right)$ is a commutative Banach algebra over $\mathbb{C}$.

Proof.

This follows directly from Proposition 2.3 and the triangle inequality. ∎

One of the reasons for using the space $\ell^{1}_{\nu}$ is to have a relatively simple and sharp convolution estimate. Another important reason is that it is easy to compute the norm of bounded linear operators. To explain how to compute the norm of a bounded linear operator on $\ell^{1}_{\nu}$ we introduce the notion of the corner points. Let $\left(e_{k}\right)_{k\in\mathbb{N}_{0}}$ denote the canonical Schauder basis for $\ell^{1}_{\nu}$, i.e. $\left(e_{k}\right)_{l}:=\delta_{kl}$ for $l\in\mathbb{N}_{0}$, so that

[TABLE]

for any $a\in\ell^{1}_{\nu}$.

Remark 2.6**.**

We shall frequently use the Schauder basis $\left(e_{k}\right)_{k\in\mathbb{N}_{0}}$ to identify an element $a\in\ell^{1}_{\nu}$ with the infinite column vector $\begin{bmatrix}a_{0}&a_{1}&\ldots&\end{bmatrix}^{T}$.

Definition 2.7**.**

The corner points $\left\{\xi_{k,\nu}\right\}_{k\in\mathbb{N}_{0}}\subset\ell^{1}_{\nu}$ of the unit ball in $\ell^{1}_{\nu}$ are defined by $\xi_{k,\nu}:=\varepsilon_{k,\nu}e_{k}$, where

[TABLE]

We shall write $\xi_{k,\nu}=\xi_{k}$ and $\varepsilon_{k,\nu}=\varepsilon_{k}$ whenever there is no chance of confusion.

The norm of a bounded linear operator on $\ell^{1}_{\nu}$ can be computed by simply evaluating it at the corner points as shown in the next proposition:

Proposition 2.8**.**

Let $\left(X,\left\|\cdot\right\|_{X}\right)$ be a normed vector space. If $\mathcal{L}\in\mathcal{B}\left(\ell^{1}_{\nu},X\right)$, then

[TABLE]

Proof.

It is clear that

[TABLE]

since $\left\|\xi_{k}\right\|_{\nu}=1$ for all $k\in\mathbb{N}_{0}$ by definition.

Conversely, let $a\in\ell^{1}_{\nu}$ be arbitrary and observe that

[TABLE]

Therefore, since $\mathcal{L}$ is bounded,

[TABLE]

Consequently,

[TABLE]

for any $a\in\ell^{1}_{\nu}$, which proves the claim. ∎

Now, suppose $\mathcal{L}\in\mathcal{B}\left(\ell^{1}_{\nu_{1}},\ell^{1}_{\nu_{2}}\right)$, where $\nu_{1},\nu_{2}>1$. Then $\mathcal{L}$ can be identified with an infinite dimensional matrix with respect to the basis $\left(e_{k}\right)_{k\in\mathbb{N}_{0}}$. More precisely, there exists unique coefficients $\left\{\mathcal{L}_{ij}\in\mathbb{C}:i,j\in\mathbb{N}_{0}\right\}$ such that

[TABLE]

Hence

[TABLE]

for any $a\in\ell^{1}_{\nu_{1}}$. In this particular setting, Proposition 2.9 can be interpreted as the statement that $\left\|\mathcal{L}\right\|_{\mathcal{B}\left(\ell^{1}_{\nu_{1}},\ell^{1}_{\nu_{2}}\right)}$ is a weighted supremum of the $\ell^{1}_{\nu_{2}}$-norms of the columns of $\mathcal{L}$. Moreover, in this case the converse of Proposition 2.9 holds as well:

Proposition 2.9**.**

Let $\nu_{1},\nu_{2}>1$ and suppose $\left\{\mathcal{L}_{ij}\in\mathbb{C}:i,j\in\mathbb{N}_{0}\right\}$ are coefficients such that the expression in (2.1) yields a well-defined linear operator $\mathcal{L}:\ell^{1}_{\nu_{1}}\rightarrow\mathbb{C}^{\mathbb{N}_{0}}$, i.e., $\left(\mathcal{L}(a)\right)_{k}$ is finite for all $a\in\ell^{1}_{\nu_{1}}$ and $k\in\mathbb{N}_{0}$. Then $\mathcal{L}\in\mathcal{B}\left(\ell^{1}_{\nu_{1}},\ell^{1}_{\nu_{2}}\right)$ if and only if $\sup_{l\in\mathbb{N}_{0}}\varepsilon_{l,\nu_{1}}\left\|\mathcal{L}_{\left(\cdot,l\right)}\right\|_{\nu_{2}}<\infty.$ Moreover, if $\mathcal{L}\in\mathcal{B}\left(\ell^{1}_{\nu_{1}},\ell^{1}_{\nu_{2}}\right)$, then

[TABLE]

Proof.

It follows directly from Proposition 2.9 that $\sup_{l\in\mathbb{N}_{0}}\varepsilon_{l,\nu_{1}}\left\|\mathcal{L}_{\left(\cdot,l\right)}\right\|_{\nu_{2}}<\infty$ whenever $\mathcal{L}\in\mathcal{B}\left(\ell^{1}_{\nu_{1}},\ell^{1}_{\nu_{2}}\right)$ and that in this case the operator norm is given by (2.2). Conversely, suppose $\sup_{l\in\mathbb{N}_{0}}\varepsilon_{l,\nu_{1}}\left\|\mathcal{L}_{\left(\cdot,l\right)}\right\|_{\nu_{2}}<\infty$. Let $a\in\ell^{1}_{\nu_{1}}$ be arbitrary, then

[TABLE]

Since the first factor in the right-hand side is finite, $\mathcal{L}:\ell^{1}_{\nu_{1}}\rightarrow\ell^{1}_{\nu_{2}}$ is thus a bounded linear operator. ∎

2.3 Multivariate sequences

In this section we introduce a space of sequences indexed by $d$-dimensional multi-indices, where $d\in\mathbb{N}$. This space will be used to analyze Taylor series of analytic functions $P:\left\{z\in\mathbb{C}:\left|z\right|\leq\nu\right\}^{d}\rightarrow\mathbb{C}^{n}$, where $\nu>0$. Such functions arise in the the analysis of local charts on (un)stable manifolds via the parameterization method developed in [MR2177465].

Formally, a sequence indexed by $d$-dimensional multi-indices is a function $p:\mathbb{N}_{0}^{d}\rightarrow\mathbb{C}$. The function $p$ is usually referred to as a $d$-dimensional array or multivariate sequence. In analogy with ordinary sequences, we shall write (as usual)

[TABLE]

Furthermore, for any multi-index $k\in\mathbb{N}_{0}^{d}$, we write $\left|k\right|=\sum_{i=1}^{d}k_{i}$, which is not to be confused with the absolute value of a (complex) number. In addition, we introduce a partial ordering $\preceq$ on $\mathbb{N}_{0}^{d}$ by

[TABLE]

We will now follow the same approach as in the previous section to set up a functional analytic framework for analyzing geometrically decaying arrays. Let $\nu>0$ and define

[TABLE]

endowed with the norm

[TABLE]

In the case that the dimension $d$ can be easily inferred from the context it will be omitted from the notation. In addition, if $n=1$ and it is clear from the context whether $p\in\ell^{1}_{\nu}$ or $p\in W^{1}_{\nu}$, we shall write $\left\|p\right\|_{W^{1}_{\nu,n,d}}=\left\|p\right\|_{\nu}$.

Next, recall that the Taylor coefficients of the product of two Taylor series is given by the one-sided discrete convolution, also referred to as the Cauchy product. More precisely, if $f,g:\left\{z\in\mathbb{C}:\left|z\right|<\nu\right\}^{d}\rightarrow\mathbb{C}$ admit power series expansions

[TABLE]

where $\tilde{f}=\left(f_{k}\right)_{k\in\mathbb{N}_{0}^{d}}$ and $\tilde{g}=\left(g_{k}\right)_{k\in\mathbb{N}_{0}^{d}}$ are $d$-dimensional arrays, then

[TABLE]

on $\left\{z\in\mathbb{C}:\left|z\right|<\nu\right\}^{d}$. In particular, the Cauchy product $\bm{\ast}$ yields a natural product structure on $W^{1}_{\nu}$. This is summarized in the following proposition.

Proposition 2.10**.**

The space $\left(W^{1}_{\nu},\bm{\ast}\right)$ is a commutative Banach algebra.

Proof.

This follows directly from the definition of $\bm{\ast}$ in (2.3) and the triangle inequality. ∎

Next, we derive an expression for the norm of a bounded linear operator on $W^{1}_{\nu}$. For this purpose we introduce a multivariate analog of the corner-points:

Definition 2.11**.**

The corner-points $\left\{\bm{\xi}_{k,d,\nu}\right\}_{k\in\mathbb{N}^{d}_{0}}\subset W^{1}_{\nu}$ of the unit ball in $W^{1}_{\nu}$ are defined by $\left(\bm{\xi}_{k,d,\nu}\right)_{l}:=\nu^{-\left|k\right|}\delta_{kl}$, where $l\in\mathbb{N}_{0}^{d}$. We shall write $\bm{\xi}_{k,d,\nu}=\bm{\xi}_{k}$ whenever there is no chance of confusion.

As before, the norm of a bounded linear operator on $W^{1}_{\nu}$ can be computed by evaluating it at the corner-points.

Proposition 2.12**.**

Let $\left(X,\left\|\cdot\right\|_{X}\right)$ be a normed vector space. If $\mathcal{L}\in\mathcal{B}\left(W^{1}_{\nu},X\right)$, then

[TABLE]

Proof.

See Proposition 2.8. ∎

3 An equivalent zero finding problem

In this section we set up a zero finding problem for establishing the existence of connecting orbits. Let us start by giving a precise description of the problem. Suppose $\tilde{p}_{0},\tilde{q}_{0}\in\mathbb{R}^{n}$ are hyperbolic equilibria of $g$. The objective is to validate an isolated connecting orbit $u$ from $\tilde{p}_{0}$ to $\tilde{q}_{0}$, which is robust with respect to “small” perturbations in $g$, by solving a boundary value problem (BVP) on a finite time domain. The method is based on the observation that a connecting orbit from $\tilde{p}_{0}$ to $\tilde{q}_{0}$ is characterized by

[TABLE]

where $L>0$ is the time of flight needed to travel from $W_{\text{loc}}^{u}\left(\tilde{p}_{0}\right)$ to $W_{\text{loc}}^{s}\left(\tilde{q}_{0}\right)$.

If $W^{u}\left(\tilde{p}_{0}\right)$ and $W^{s}\left(\tilde{q}_{0}\right)$ intersect transversally along $u$, then the connecting orbit is robust, i.e., it will persist for sufficiently “small” perturbations in $g$. In this case, the intersection $W^{u}\left(\tilde{p}_{0}\right)\cap W^{s}\left(\tilde{q}_{0}\right)\cap U$, where $U$ is a neighborhood of the connecting orbit in which it is unique, is necessarily an one dimensional manifold. Hence, by counting dimensions, a necessary condition for the existence of a transverse isolated connecting orbit is

[TABLE]

This condition is often referred to as a non-degeneracy condition for connecting orbits. We shall henceforth assume that this condition is satisfied. In particular, we do not assume a-priori that the connecting orbit is isolated and transverse. Instead, we will obtain these properties from the proof of existence (a contraction argument), see Proposition 3.19.

We start by setting up equations for local charts on the (un)stable manifolds by using the parameterization method [MR2177465] and the methodology presented in [ManifoldTaylor]. These charts will be used to supplant the boundary conditions in (3.1) with explicit equations. Next, we set up an equivalent system of equations for the differential equation by using Chebyshev series and domain decomposition as explained in [domaindecomposition]. Finally, in order for the resulting zero finding problem to be well posed, we complete the system of equations by imposing appropriate phase conditions.

3.1 Charts on the (un)stable manifolds

In this section we give a brief overview of the method developed in [ManifoldTaylor] to compute local charts on the (un)stable manifolds. The reader is referred to [ManifoldTaylor] for a more detailed exposition of the theory. We consider the computation of a local chart on the stable manifold of $\tilde{q}_{0}$. A chart on the unstable manifold of $\tilde{p}_{0}$ can be computed in the same way by reversing the sign of the vector field.

The Parameterization Method

The idea of the parameterization method [MR2177465] is to construct a diffeomorphism which conjugates the nonlinear dynamics on the stable manifold to an easier and fully understood flow $\psi$. For the sake of simplicity, let us assume that $Dg\left(\tilde{q}_{0}\right)$ is diagonalizable. This assumption is, however, not necessary, as we will explain in a moment.

Let $\lambda_{1}^{s},\ldots,\lambda^{s}_{n_{s}}\in\mathbb{C}$ be the stable eigenvalues of $Dg\left(\tilde{q}_{0}\right)$. If all eigenvalues are real and semisimple, then there exists neighborhoods $U\subset\mathbb{R}^{n}$ and $V\subset\mathbb{R}^{n_{s}}$ of $\tilde{q}_{0}$ and [math], respectively, such that the dynamics on $W^{s}\left(\tilde{q}_{0}\right)\cap U$ is conjugate to the flow

[TABLE]

If some of the eigenvalues are complex, however, special care has to be taken. Let us for the moment forget about this technicality and consider the complex dynamics generated by $u^{\prime}=g(u)$ on $\mathbb{C}^{n}$. Then the dynamics on the complex local stable manifold, which we denote by $W_{\text{loc}}^{s,c}\left(\tilde{q}_{0}\right)$, is conjugate to the flow $\psi$ restricted to the polydisk

[TABLE]

for some sufficiently small $\nu_{s}>0$.

The idea is to find an analytic map $Q:\mathbb{B}_{\nu_{s}}\rightarrow\mathbb{C}^{n}$ which conjugates the nonlinear flow $\varphi$ on $W^{s,c}_{\text{loc}}\left(\tilde{q}_{0}\right)$ to the linear flow $\psi$ on $\mathbb{B}_{\nu_{s}}$ for $t\geq 0$. In other words, we seek a map $Q$ such that

{\mathbb{B}_{\nu_{s}}}$${\mathbb{C}^{n}}$${\mathbb{B}_{\nu_{s}}}$${\mathbb{C}^{n}}$$\scriptstyle{\psi}$$\scriptstyle{Q}$$\scriptstyle{\varphi}$$\scriptstyle{Q}

commutes, i.e., $Q\left(\psi\left(t,\phi\right)\right)=\varphi\left(t,Q\left(\phi\right)\right)$ for all $(t,\phi)\in\mathbb{R}_{\geq 0}\times\mathbb{B}_{\nu_{s}}$. Differentiation of this relation at $t=0$ yields the so-called invariance equation:

[TABLE]

Note that this equation does not depend on time anymore. Moreover, it is easy to see that if the invariance equation holds, then

[TABLE]

is an orbit in $W^{s,c}\left(\tilde{q}_{0}\right)$ for any $\phi\in\mathbb{B}_{\nu_{s}}$, i.e., $Q:\mathbb{B}_{\nu_{s}}\rightarrow W^{s,c}_{\text{loc}}\left(\tilde{q}_{0}\right)$ (see [ManifoldTaylor, Lemma $2.6$]). Therefore, the problem of computing a chart is now reduced to solving (3.3).

Solving the invariance equation

Since $Q$ is assumed to be analytic on $\mathbb{B}_{\nu_{s}}$, i.e., $Q$ is analytic on a slightly larger open neighborhood of $\mathbb{B}_{\nu_{s}}$, there exist coefficients $q\in W^{1}_{\nu_{s},n,n_{s}}$ such that

[TABLE]

Observe that the zeroth order Taylor coefficient is necessarily the equilibrium, i.e., $q_{0}=\tilde{q}_{0}$. Furthermore, since $Q$ is assumed to be a diffeomorphism, it must hold that

[TABLE]

where $E_{s}$ is the stable eigenspace of $Dg\left(q_{0}\right)$. In other words, the first order Taylor coefficients $\left\{q_{k}:\left|k\right|=1,\ k\in\mathbb{N}_{0}^{n_{s}}\right\}$ are the eigenvectors of $Dg\left(q_{0}\right)$. Note that these are only determined up to a scaling.

To determine the higher order Taylor coefficients $\left\{q_{k}:\left|k\right|\geq 2,\ k\in\mathbb{N}_{0}^{n_{s}}\right\}$, we first introduce the map $C:W^{1}_{\nu_{s},n}\rightarrow W^{1}_{\nu_{s},n}$ defined by

[TABLE]

where the latter product is understood to be the one-sided discrete convolution. Furthermore, $\mathcal{A}\subset\mathbb{N}_{0}^{n}$ and $\left\{g_{\alpha}:\alpha\in\mathcal{A}\right\}\subset\mathbb{R}^{n}$ are the coefficients of $g$ in the monomial basis. In particular, observe that

[TABLE]

since the Taylor coefficients of the product of two Taylor expansions is given by the one-sided discrete convolution. Formally, we should incorporate the weight $\nu_{s}$ and dimension $n_{s}$ into the notation for $C$. However, since these parameters can usually be inferred from the context and we wish to use the same notation for the unstable manifold, we have chosen to omit them from the notation.

Substitution of the Taylor expansion for $Q$ into (3.3) yields the following system of equations:

[TABLE]

where $\left\langle\cdot,\cdot\right\rangle$ denotes the standard Hermitian inner product on $\mathbb{C}^{n_{s}}$. We shall use this system of equations to set up a zero finding problem for computing a chart on $W^{s,c}_{\text{loc}}\left(\tilde{q}_{0}\right)$. Before we proceed, observe that (3.6) is equivalent to

[TABLE]

Moreover, differentiation of (3.5) at $\phi=0$ shows that $Dg\left(q_{0}\right)q_{k}=C_{k}(q)$ for $\left|k\right|=1$. Hence (3.7) reduces to the eigenvalue/eigenvector equation for $Dg\left(q_{0}\right)$ for $\left|k\right|=1$. Similarly, repeated differentiation of (3.5) at $\phi=0$ shows that the right-hand-side of (3.7) only depends on Taylor coefficients of order strictly below $\left|k\right|$. In conclusion, the Taylor coefficients $q$ can be computed recursively up to any desired order provided

[TABLE]

The latter condition is usually referred to as a non-resonance condition and is related to the regularity of $Q$. More precisely, in the presence of a resonance, the parameterization method, as applied above with the linear “model” flow $\psi$, does not yield an analytic conjugation $Q$. It is explained in [ManifoldTaylor] how to construct an analytic conjugation in the present of a resonance. The idea is to use a nonlinear normal form for $\psi$ instead of just the linear flow in (3.2). The interested reader is referred to [ManifoldTaylor] for a detailed exposition of the resonant case. For the sake of presentation, however, we shall assume throughout this paper that there are no resonances.

We are now ready to set up a zero finding problem for computing the Taylor coefficients of $Q$ (which include the equilibrium and eigenvectors) and the eigenvalues $\lambda^{s}$:

Definition 3.1** (Taylor map for stable manifolds).**

Let $0<\tilde{\nu}_{s}<\nu_{s}$ be given weights. The Taylor map $F_{Q}:\mathbb{C}^{n_{s}}\times W^{1}_{\nu_{s},n}\rightarrow W^{1}_{\tilde{\nu}_{s},n}$ for stable manifolds is defined by

[TABLE]

Remark 3.2**.**

The latter map is well-defined since $\left(\left\langle\lambda^{s},k\right\rangle q_{k}\right)_{k\in\mathbb{N}_{0}^{n_{s}}}\in W^{1}_{\tilde{\nu}_{s},n}$ for any $q\in W^{1}_{\nu_{s},n}$ and $0<\tilde{\nu}_{s}<\nu_{s}$.

Remark 3.3**.**

If $\left(\lambda^{s},q\right)$ is a zero of $F_{Q}$, then so is $\left(\lambda^{s},\mu q\right)$, where $\mu\in\mathbb{C}^{n_{s}}$ and $\left(\mu q\right)_{k}:=\mu^{k}q_{k},$ see [ManifoldTaylor, Lemma $2.2$]. We will get rid of this extra degree of freedom by fixing the orientation and length of the eigenvectors. In particular, observe that the scaling of the eigenvectors (and in turn the “scaling” of $q$) determines the decay rate of the coefficients $q$ and hence the size of $\mathbb{B}_{\nu_{s}}$. In effect, the length of the eigenvectors determine (roughly speaking) the “size” of the patch on $W_{\text{loc}}^{s,c}\left(q_{0}\right)$ parameterized by $Q$. The interested reader is referred to [MR3437754] for a more thorough explanation where this phenomena is explored in detail.

The zero finding problem for the computation of a chart on the unstable manifold is set up in an analogous way. For the sake of completeness (and introducing notation) let us explicitly state the assumptions and the associated zero finding map. We assume that $Dg\left(\tilde{p}_{0}\right)$ is diagonalizable and that the associated eigenvalues $\lambda^{u}\in\mathbb{C}^{n_{u}}$ satisfy the non-resonance condition (3.8). The goal is to compute a parameterization $P:\mathbb{B}_{\nu_{u}}\subset\mathbb{C}^{n_{u}}\rightarrow W^{u,c}_{\text{loc}}\left(\tilde{p}_{0}\right)$ of the form $P(\theta)=\sum_{k\in\mathbb{N}_{0}^{n_{u}}}p_{k}\theta^{k}$, where $\nu_{u}>0$, by finding a zero of the following map:

Definition 3.4** (Taylor map for unstable manifolds).**

Let $0<\tilde{\nu}_{u}<\nu_{u}$ be given weights. The Taylor map $F_{P}:\mathbb{C}^{n_{u}}\times W^{1}_{\nu_{u},n}\rightarrow W^{1}_{\tilde{\nu}_{u},n}$ for unstable manifolds is defined by

[TABLE]

Symmetry

In the preceding exposition we considered the complex dynamical system $u^{\prime}=g(u)$ on $\mathbb{C}^{n}$. Our main interest, however, is the computation of invariant manifolds in the real-valued dynamical system on $\mathbb{R}^{n}$. We will now explain how we can recover charts for the (un)stable manifolds in the real system from the complex ones through the use of symmetry. We remark that one could also have set up the parameterization method in the real-valued setting from the start. However, in that case, we would have had to separate the cases between the presence of complex eigenvalues and a completely real spectrum. It is in our opinion more convenient from both a practical and theoretical point of view to develop a unified approach.

Let us consider the stable manifold again. Observe that complex eigenvalues will always appear in conjugate pairs, since $g$ is real-analytic. Suppose there are $d_{s}$ complex conjugate pairs of eigenvalues and $n_{s}-2d_{s}$ real ones. Furthermore, assume that we have ordered the eigenvalues $\lambda^{s}$ in such a way that $\lambda^{s}_{i}=\overline{\lambda^{s}_{i+1}}$ for $i\in\left\{2l+1:0\leq l\leq d_{s}-1\right\}$. Next, define the map $\Sigma:\mathbb{C}^{n_{s}}\rightarrow\mathbb{C}^{n_{s}}$ by

[TABLE]

and note that $\Sigma$ is an involution on $\mathbb{C}^{n_{s}}$. For this reason we shall frequently write $z^{\star}:=\Sigma(z)$. In particular, note that we ordered the stable eigenvalues in such a way that $\left(\lambda^{s}\right)^{\star}=\lambda^{s}$. Finally, we extend this notion of involution to $W^{1}_{\nu_{s},n}$ by defining

[TABLE]

It can be readily seen from (3.9) that $k^{\star}=\Sigma(k)$ is simply a permutation of the multi-index $k$.

The key observation for obtaining charts for the real manifolds is stated in the following lemma. The proof can be found in [ManifoldTaylor, Lemma $2.1$].

Proposition 3.5**.**

If $q\in W^{1}_{\nu_{s},n}$ is symmetric, i.e., $q^{\star}=q$, then the map $Q:\mathbb{B}_{\nu_{s}}\rightarrow\mathbb{C}^{n}$ defined by

[TABLE]

is real valued on the set $\mathbb{B}^{\text{sym}}_{\nu_{s}}:=\left\{\phi\in\mathbb{B}_{\nu_{s}}:\phi^{\star}=\phi\right\}$. In addition, if $\lambda^{s}\in\mathbb{C}^{n_{s}}$ is symmetric and $F\left(\lambda^{s},q\right)=0$, then $Q|_{\mathbb{B}^{\text{sym}}_{\nu_{s}}}$ is a parameterization of the real stable manifold $W^{s}_{\text{loc}}\left(q_{0}\right)$.

Remark 3.6**.**

Note that $\mathbb{B}^{\text{sym}}_{\nu_{s}}$ is a real manifold of dimension $n_{s}$. More precisely, we can identify $\mathbb{B}^{\text{sym}}_{\nu_{s}}$ with the (real) manifold

[TABLE]

by using the (linear) map $\iota_{s}:\mathbb{R}^{n_{s}}\rightarrow\mathbb{C}^{n_{s}}$ defined by

[TABLE]

Remark 3.7**.**

Strictly speaking, when $q^{\star}=q$ the assumption that $\lambda^{s}$ is symmetric is not necessary. To see this, let $e^{s}_{i}\in\mathbb{C}^{n_{s}}$ be the unit vector defined by $\left[e^{s}_{i}\right]_{j}=\delta_{ij}$, where $1\leq i,j\leq n_{s}$. If $F\left(\lambda^{s},q\right)=0$ for some $\lambda^{s}\in\mathbb{C}^{n_{s}}$, then

[TABLE]

In particular, if we take the complex conjugate of the lefthand-side of the above expression for $k=e^{s}_{2j-1}$, we obtain

[TABLE]

since $g$ is real-analytic and $q^{\star}=q$ by assumption. On the other hand,

[TABLE]

Hence it follows that $\left(\lambda^{s}\right)^{\star}=\lambda^{s}$.

In conclusion, in order to conclude that a point lies on the real stable manifold, it suffices to verify that $q^{\star}=q$. It is explained in Section 4 how we can verify this in practice. For now, let us mention that the verification is based on the following observation whose proof can be found in [ManifoldTaylor, Lemma $4.1$]:

Lemma 3.8**.**

The map $F_{Q}$ is compatible with $\star$, i.e., $F_{Q}\left(\left(\lambda^{s}\right)^{\star},q^{\star}\right)=F_{Q}\left(\lambda^{s},q\right)^{\star}$ for any $\left(\lambda^{s},q\right)\in\mathbb{C}^{n_{s}}\times W^{1}_{\nu_{s},n}$.

Analogous results hold for the parameterization of the unstable manifold of $\tilde{p}_{0}$. To avoid clutter in the notation we shall denote the involution associated to the unstable manifold of $\tilde{p}_{0}$ by $\star$ as well.

3.2 Chebyshev series and domain decomposition

In this section we follow the strategy in [domaindecomposition] to recast the differential equation into an equivalent zero finding problem on $\ell^{1}_{\nu,n}$. The reader is referred to [domaindecomposition] for the details. Let $\mathcal{P}_{m}:=\left\{t_{0}=0<t_{1}<\ldots<\right.$ $\left.t_{m}=1\right\}$, where $m\in\mathbb{N}$, be any partition of $[0,1]$. Then the differential equation in (3.1) is equivalent to

[TABLE]

where $2\leq i\leq m$. These boundary conditions are thus imposed on the internal nodes of the partition $\mathcal{P}_{m}$ only. The boundary conditions at the end points will be discussed in Section 3.3. If $\left(\mathbb{P}_{i}\right)_{i=1}^{m}$ admits a solution, then each $u_{i}$ is real-analytic since $g$ is. Therefore, there exists weights $\nu_{i}>1$ and real coefficients $a^{i}\in\ell^{1}_{\nu_{i},n}$ such that

[TABLE]

in $C\left(\left[t_{i-1},t_{i}\right]\right)$. Here $\left(T^{i}_{k}\right)_{k\in\mathbb{N}_{0}}$ are the shifted Chebyshev-polynomials on $\left[t_{i-1},t_{i}\right]$ defined by

[TABLE]

Next, define the map $c:\ell^{1}_{\nu_{i},n}\rightarrow\ell^{1}_{\nu_{i},n}$ by

[TABLE]

where the latter product is understood to be the symmetric discrete convolution $\ast$. Then

[TABLE]

since the Chebyshev coefficients of the product of two functions is given by the symmetric discrete convolution. Formally, we should incorporate the index $i$ into the notation for $c$ to emphasize its dependence on the weight $\nu_{i}$. However, since the domain of $c$ can usually be easily inferred from the context, we haven chosen to omit the index from the notation.

Remark 3.9**.**

Throughout this paper we will need to analyze $Dc\left(a^{i}\right)$, where $a^{i}\in\ell^{1}_{\nu_{i},n}$ and $1\leq i\leq m$, on numerous occasions. For this reason, we state here for future reference how this derivative can be computed in an efficient way. Since $\left(\ell^{1}_{\nu_{i},n},\ast\right)$ is a Banach algebra, we may use the “usual” rules of calculus to compute the derivative of $c$. In particular, direct differentiation of (3.11) with respect to $a^{i}$ shows that

[TABLE]

where $\left[c\right]_{j}$ denotes the $j$-th component of $c$ and $\hat{g}^{ijl}\in\ell^{1}_{\nu_{i}}$ are the Chebyshev coefficients of

[TABLE]

Substitution of the Chebsyshev expansions for $\left(u_{i}\right)_{i=1}^{m}$ into $\left(\mathbb{P}_{i}\right)_{i=1}^{m}$ yields an equivalent system of equations for the coefficients and gives rise to the following map:

Definition 3.10** (Chebyshev map for ODEs).**

Let $\left(\nu_{i}\right)_{i=1}^{m}$ and $\left(\tilde{\nu}_{i}\right)_{i=1}^{m}$ be collections of weights such that $1<\tilde{\nu}_{i}<\nu_{i}$ for all $1\leq i\leq m$. The Chebyshev map for ODEs is the function $F_{u}:\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}\rightarrow\ell^{1}_{\tilde{\nu}_{1},n}/\mathbb{C}^{n}\oplus\bigoplus_{i=2}^{m}\ell^{1}_{\tilde{\nu}_{i},n}$ defined by

[TABLE]

where $a=\left(a^{1},\ldots,a^{m}\right)$, $f_{1}:\ell^{1}_{\nu_{1},n}\rightarrow\ell^{1}_{\tilde{\nu}_{1},n}/\mathbb{C}^{n}$ is given by

[TABLE]

and $f_{i}:\ell^{1}_{\nu_{i-1},n}\times\ell^{1}_{\nu_{i},n}\rightarrow\ell^{1}_{\tilde{\nu}_{i},n}$ by

[TABLE]

for $2\leq i\leq m$.

Remark 3.11**.**

The map $F_{u}$ is well-defined, since $\left(ka^{i}_{k}\right)_{k\in\mathbb{N}_{0}}\in\ell^{1}_{\tilde{\nu}_{i},n}$ for any $a^{i}\in\ell^{1}_{\nu_{i},n}$ and $1<\tilde{\nu}_{i}<\nu_{i}$.

Let us stress the subtle difference between the latter map and the one constructed in [domaindecomposition]; in the current setting we allow for complex Chebyshev coefficients. The main reason for this is that the parameterization maps $P$ and $Q$ (see the previous section) are in principle complex-valued and we wish to use them to proof that $u(0)\in W^{u}_{\text{loc}}\left(\tilde{p}_{0}\right)$ and $u(L)\in W^{s}_{\text{loc}}\left(\tilde{q}_{0}\right)$. We will conclude a-posteriori that the Chebyshev coefficients are in fact real by invoking symmetry arguments. Indeed, for any $a=\left(a^{1},\ldots,a^{m}\right)$, define $\overline{a}:=\left(\overline{a^{1}},\cdots,\overline{a^{m}}\right)$ by $\left(\overline{a^{i}}\right)_{k}:=\overline{a^{i}_{k}}$. We will conclude that an element $a$ is real, i.e., $\overline{a}=a$, by using the following observation:

Lemma 3.12**.**

The map $F_{u}$ is compatible with conjugation, i.e., $F_{u}\left(\overline{a}\right)=\overline{F_{u}(a)}$.

Proof.

Let $a=\left(a^{1},\ldots,a^{m}\right)\in\prod_{j=1}^{m}\ell^{1}_{\nu_{j},n}$ be arbitrary. We will prove that $c\left(\overline{a^{j}}\right)=\overline{c\left(a^{j}\right)}$ for $1\leq j\leq m$. The desired result follows directly from this observation. Recall that a Chebyshev series is a Fourier series up to coordinate transformation. To be more precise, let $1\leq j\leq m$, set $a^{j}_{-k}:=a^{j}_{k}$ for $k\in\mathbb{N}$, and define

[TABLE]

Then

[TABLE]

by (3.10) and the definition of the Chebyshev polynomials (see Definition 2.1). Note that the latter series converges uniformly to an analytic $2\pi$-periodic function, since $a^{j}\in\ell^{1}_{\nu_{j},n}$. Furthermore, since $a^{j}_{-k}=a^{j}_{k}$ (by definition), it follows that

[TABLE]

where $c$ is defined in (3.11) and we have set $c_{-k}\left(a^{j}\right):=c_{k}\left(a^{j}\right)$ for $k\in\mathbb{N}$. In particular,

[TABLE]

On the other hand, since $g$ is real-analytic, similar reasoning shows that

[TABLE]

Therefore, since a (pointwise) convergent Fourier series is unique, we conclude that $c\left(\overline{a^{j}}\right)=\overline{c\left(a^{j}\right)}$. ∎

Finally, observe that by construction we now have the following result:

Proposition 3.13**.**

Suppose $a\in\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}$ is symmetric, i.e., $\overline{a}=a$, then $F(a)=0$ if and only if the functions $\left\{u_{i}=a^{i}_{0}+2\sum_{k=1}^{\infty}a^{i}_{k}T^{i}_{k}:1\leq i\right.$ $\left.\leq m\right\}$ constitute a solution of $\left(\mathbb{P}_{i}\right)_{i=1}^{m}$.

3.3 The connecting orbit map

In this section we set up a zero finding problem for (3.1). We have already set up appropriate zero finding mappings for the ODE and charts on the (un)stable manifolds. What remains is imposing appropriate phase conditions.

Boundary conditions

We can now replace the boundary conditions in (3.1) with explicit equations. Let $P$ and $Q$ denote the local parameterizations of the complex unstable and stable manifold as before, respectively. Then the conditions $u(0)\in W^{u}_{\text{loc}}\left(\tilde{p}_{0}\right)$ and $u(1)\in W^{s}_{\text{loc}}\left(\tilde{q}_{0}\right)$ are equivalent to the problem of finding coordinates $\theta\in\mathbb{B}^{\text{sym}}_{\nu_{u}}$ and $\phi\in\mathbb{B}^{\text{sym}}_{\nu_{s}}$ such that

[TABLE]

As mentioned before, we will verify a-posteriori that $\theta^{\star}=\theta$ and $\phi^{\star}=\phi$ so that $P\left(\theta\right)$ and $Q\left(\theta\right)$ are points on the real (un)stable manifolds.

Length of the eigenvectors

Recall that the first order Taylor coefficients of $P$ and $Q$ are determined up to a rescaling. To get rid of this extra degree of freedom we prescribe the length and orientation of the eigenvectors of $Dg\left(\tilde{p}_{0}\right)$ and $Dg\left(\tilde{q}_{0}\right)$. More precisely, we require that

[TABLE]

where $\hat{p}_{k},\hat{q}_{k}\in\mathbb{C}^{n}$ are prescribed vectors and $\epsilon_{u,k},\epsilon_{s,k}>0$. In practice, $\hat{p}_{k}$ and $\hat{q}_{k}$ are numerical approximations of the eigenvectors of $Dg\left(\tilde{p}_{0}\right)$ and $Dg\left(\tilde{q}_{0}\right)$, respectively, and $\epsilon_{u,k},\epsilon_{s,k}$ are their respective squared lengths. In order to respect the symmetry $\star$, we impose the following ordering:

[TABLE]

We recall again that the length of the eigenvectors determines the decay rate of the Taylor coefficients and hence the size of the domains of $P$ and $Q$.

Translation invariance in time

Finally, we introduce a phase condition to fix the time parameterization of the connecting orbit. In [MR635945, MR2359336] a phase condition specifically tailored for continuation is presented. The idea is to fix a reference function $\tilde{u}$ and to minimize the functional

[TABLE]

on some appropriate functions space, where $u$ is the connecting orbit. This phase condition is used in popular software packages for continuation such as AUTO and Matcont and was first suggested in [MR635945].

The intuition is that this phase condition enforces the connecting orbit to remain as close as possible (in the $L^{2}$-sense) to the reference solution with respect to small shifts in time. In practice, $\tilde{u}$ is the solution computed at the previous continuation step (or just the numerical approximation $\hat{u}$ in case we are not performing continuation). In particular, a necessary condition for the latter functional to have a minimum at $s=0$ is

[TABLE]

We shall use (3.16) to construct an appropriate phase condition in terms of Chebyshev coefficients by approximating the integral on a finite domain. First, write

[TABLE]

If the time of flight $L>0$ is sufficiently large, then

[TABLE]

Next, write

[TABLE]

For notational convenience, let us omit the superscripts from the Chebyshev coefficients and assume (for the moment) that $u_{i}$ and $\tilde{u}_{i}$ are scalar functions. Then

[TABLE]

where $\left\langle\cdot,\cdot\right\rangle_{L^{2}}$ denotes the standard complex inner product on $L^{2}\left(\left[t_{i-1},t_{i}\right]\right)$.

Now, rescale time back to $[-1,1]$, use the coordinate transformation $\theta=\arccos\left(t\right)$ and the definition of the Chebyshev polynomials to see that

[TABLE]

for any $k,l\in\mathbb{N}_{0}$. Finally, substitution of the latter expression into (3.17) yields

[TABLE]

If $u_{i}$ and $\tilde{u}_{i}$ are vector-valued, then we need to carry out the above computations component-wise and sum over the components.

In practice, we choose the Chebyshev coefficients of $\tilde{u}$ to be real, since in the end we wish to establish the existence of a real-valued connecting orbit. Altogether, this motivates the following definition:

Definition 3.14** (Phase condition for translation invariance in time).**

Let $b=\left(b^{1},\ldots,\right.$ $\left.b^{m}\right)\in\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}$ be given symmetric sequences, i.e., $\bar{b}=b$, such that $b^{i}_{k}=0$ for $k\geq N_{i}$, for some $N\in\mathbb{N}^{m}$. The phase condition for translation invariance in time is the map $\eta:\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}\rightarrow\mathbb{C}$ defined by the following truncated version of (3.18):

[TABLE]

Remark 3.15**.**

The expression for $\eta$ might seem complicated at first sight. Note, however, that $\eta$ is really just an affine linear map depending on finitely many components $a^{i}_{k}$ only.

We are now ready to set up the connecting orbit map. To this end, let

[TABLE]

be given weights such that $\tilde{\nu}<\nu$, $\nu_{u},\nu_{s}>0$, $\nu_{i}>1$ for $1\leq i\leq m$, and set

[TABLE]

Definition 3.16** (Chebyshev-Taylor map for connecting orbits).**

The Chebyshev-Taylor map $F:\mathcal{X}_{\nu}\rightarrow\mathcal{Y}_{\tilde{\nu}}$ for connecting orbits is defined by

[TABLE]

where $x=\left(\theta,\phi,\lambda^{u},\lambda^{s},a^{1},\ldots,a^{m},p,q\right)$.

Remark 3.17**.**

We shall frequently denote elements in $\mathcal{Y}_{\tilde{\nu}}$ by

[TABLE]

where

• •

$y_{t_{0}},y_{t_{m}}\in\mathbb{C}^{n}$* correspond to the equations for the boundary conditions at $t=t_{0}$ and $t=t_{m}$, respectively,*

• •

$y_{\hat{p}_{1}}\in\mathbb{C}^{n_{u}}$, $y_{\hat{q}_{1}}\in\mathbb{C}^{n_{s}}$ correspond to the equations for fixing the length and orientation of the eigenvectors,

• •

$y_{\eta}\in\mathbb{C}$* corresponds to the phase condition for fixing the time parameterization of the orbit,*

• •

$y_{a},y_{p}$* and $y_{q}$ correspond to the Chebyshev and Taylor coefficients, respectively, as before.*

The only reason for introducing the weights $\tilde{\nu}$ is to specify the codomain of $F$. These weights are irrelevant though, since we will establish the existence of a connecting orbit by analyzing a fixed point map from $\mathcal{X}_{\nu}$ into itself, see Section 4.2. For this reason, we only specify a norm on $\mathcal{X}_{\nu}$. Namely, we set

[TABLE]

where $x=\left(\theta,\phi,\lambda^{u},\lambda^{s},a^{1},\ldots,a^{m},p,q\right)$.

Symmetry revisited

Next, we examine the compatibility of $F$ with respect to the symmetries introduced in the previous sections. In particular, the involution operations on the space of Taylor and Chebyshev coefficients yield an symmetry operation $\bm{\star}$ on $\mathcal{X}_{\nu}$ defined by

[TABLE]

Similarly, we define an involution on the range $\mathcal{Y}_{\tilde{\nu}}$, also denoted by $\bm{\star}$, via

[TABLE]

Lemma 3.18**.**

The map $F$ is compatible with $\bm{\star}$, i.e., $F\left(x^{\bm{\star}}\right)=F(x)^{\bm{\star}}$ for any $x\in\mathcal{X}_{\nu}$.

Proof.

We start by considering the phase conditions associated to the unstable manifold. First, observe that

[TABLE]

by definition of $\star$ and reordering of the series associated to the unstable manifold. Furthermore, note that

[TABLE]

since $\hat{p}_{k}$ was ordered in a symmetric way, see (3.14). The computations for the stable manifold are analogous. Next, observe that $\eta\left(\overline{a}\right)=\overline{\eta\left(a\right)}$, since the Chebyshev coefficients $b$ of the reference orbit $\tilde{u}$ are real. Finally, recall that $F_{u}\left(\overline{a}\right)=\overline{F_{u}\left(a\right)}$, $F_{P}\left(\left(\lambda^{u}\right)^{\star},p^{\star}\right)=F_{P}(\lambda^{u},p)^{\star}$ and $F_{Q}\left(\left(\lambda^{s}\right)^{\star},q^{\star}\right)=F_{Q}(\lambda^{s},q)^{\star}$ by Lemmas 3.12 and 3.8, respectively. Altogether, this proves the result. ∎

We are now ready to formulate an appropriate characterization of a connecting orbit:

Proposition 3.19**.**

Suppose $F$ has a unique zero $x$ in some open neighborhood $U\subset\mathcal{X}_{\nu}$ and assume that $x^{\bm{\star}}\in U$. Then $p_{0},q_{0}\in\mathbb{R}^{n}$ are equilibria of $g$, $P\left(\theta\right)\in W^{u}_{\text{loc}}\left(p_{0}\right)$, $Q\left(\phi\right)\in W^{s}_{\text{loc}}\left(q_{0}\right)$, and the map $u$ defined by the Chebyshev coefficients is an isolated connecting orbit from $p_{0}$ to $q_{0}$.

Proof.

Suppose $F(x)=0$, then the previous lemma implies that $x^{\bm{\star}}=x$, since $x$ is the only zero in $U$ and $x^{\bm{\star}}\in U$. Consequently, the Chebyshev coefficients $a$ are real and $P(\theta)$ and $Q(\phi)$ are points on the real (un)stable manifolds by Proposition 3.5. Therefore, $u$ is a connecting orbit from $p_{0}$ to $q_{0}$ by Proposition 3.13. Moreover, the connecting orbit is isolated, since $x$ is. ∎

Remark 3.20**.**

In practice, we seek a zero of $F$ in a closed ball $B_{r}\left(\hat{x}\right)$ of radius $r>0$ centered at an approximate zero $\hat{x}$ obtained through numerical simulation. The numerical computations yield an approximate zero which is almost symmetric (up to machine precision). We enforce that $\left(\hat{x}\right)^{\star}=\hat{x}$ by going through “all” the elements of $\hat{x}$ and imposing the exact symmetry conditions. For example, for the Taylor coefficients $\hat{q}$, we determine all the multi-indices $k\in\mathcal{K}^{s}$ such that $k^{\star}\in\mathcal{K}^{s}$ and then redefine $\hat{q}_{k^{\star}}$, for each $k^{\star}\not=k$, by setting it equal to $\overline{\hat{q}_{k}}$ (if $k^{\star}=k$ we set it equal to $\text{Re}\left(\hat{q}_{k}\right)$). The symmetry implies that $\left\|x-\hat{x}\right\|_{\mathcal{X}_{\nu}}=\left\|x^{\star}-\hat{x}\right\|_{\mathcal{X}_{\nu}}$ for all $x\in\mathcal{X}_{\nu}$. Hence $B_{r}\left(\hat{x}\right)^{\star}=B_{r}\left(\hat{x}\right)$, which motivates the assumption that $x^{\star}\in U$.

Transversality

We end this section with a sufficient condition for proving that a connecting orbit is transverse. The key observation is summarized in the following lemma:

Lemma 3.21**.**

Suppose $a,\tilde{a}\in\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}$ are real. Let $u,w:[0,1]\rightarrow\mathbb{R}^{n}$ denote the maps associated to $a$ and $\tilde{a}$, respectively, i.e.,

[TABLE]

Then $\tilde{a}\in\text{ker}\left(DF_{u}\left(a\right)\right)$ if and only if $w^{\prime}(t)=LDg\left(u(t)\right)w(t)$ on $[0,1]$.

Proof.

A straightforward computations shows that $\tilde{a}\in\text{ker}\left(DF_{u}\left(a\right)\right)$ if and only if

[TABLE]

Furthermore, substitution of the expression in (3.12) for $Dc\left(a^{i}\right)\tilde{a}^{i}$ shows that the above system of equations is equivalent to

[TABLE]

which proves the statement (see Proposition 3.13). ∎

We are now ready to formulate a sufficient criterium for establishing the transversality of a connecting orbit.

Proposition 3.22**.**

Suppose $x\in\mathcal{X}_{\nu}$ is symmetric and $F(x)=0$. If $DF(x)$ is injective, then $x$ corresponds to a transverse connecting orbit.

Proof.

It is shown in Proposition 3.19 that $x$ corresponds to a connecting orbit $u$ from $p_{0}$ to $q_{0}$ with the property that $u(0)=P\left(\theta\right)\in W^{u}_{\text{loc}}\left(p_{0}\right)$ and $u(1)=Q\left(\phi\right)\in W^{s}_{\text{loc}}\left(q_{0}\right)$. To show that $u$ is transverse, first observe that the mappings

[TABLE]

are parameterizations of $W_{\text{loc}}^{u}\left(p_{0}\right)$ and $W_{\text{loc}}^{s}\left(q_{0}\right)$, respectively, by Proposition 3.5 and Remark 3.6. Hence teh amp $\tilde{\theta}\mapsto\varphi\left(t,P\circ\iota_{u}\left(\tilde{\theta}\right)\right)$, where $\varphi$ denotes the flow generated by $Lg$, is a diffeomorphism from $\mathbb{B}^{\text{sym},\text{re}}_{\nu_{u}}$ into $W^{u}\left(p_{0}\right)$ for any $t\in\mathbb{R}$. Therefore, its derivative $D_{x}\varphi\left(t,P\left(\theta\right)\right)D_{\theta}P\left(\theta\right)\iota_{u}$ (evaluated at $\tilde{\theta}=\iota_{u}^{-1}\left(\theta\right)$) is an isomorphism from $\mathbb{R}^{n_{u}}$ onto $T_{\varphi\left(t,P\left(\theta\right)\right)}W^{u}\left(p_{0}\right)$. Similarly, $D_{\phi}Q\left(\phi\right)\iota_{s}$ is an isomorphism from $\mathbb{R}^{n_{s}}$ onto $T_{Q\left(\phi\right)}W^{s}\left(q_{0}\right)$. Consequently, since $\varphi\left(1,P\left(\theta\right)\right)=Q\left(\phi\right)=u(1)$, the linear map

[TABLE]

is a surjection from $\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{s}}=\mathbb{R}^{n+1}$ onto $\left(T_{u(1)}W^{u}\left(p_{0}\right)+T_{u(1)}W^{s}\left(q_{0}\right)\right)\subset\mathbb{R}^{n}$.

Now, suppose $DF(x)$ is injective but $u$ is not transverse. Then the intersection of $W^{u}\left(p_{0}\right)$ and $W^{s}\left(q_{0}\right)$ is (in particular) not transverse at $u(1)$, since $u$ is transverse if and only if it is transverse at a point. Hence the map $\Phi_{1}:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n}$ cannot be surjective. Therefore, $\dim\left(\text{ker}\left(\Phi_{1}\right)\right)\geq 2$. Consequently, there exist two linearly independent vectors $\begin{bmatrix}\tilde{\theta}_{1}\\ \tilde{\phi}_{1}\end{bmatrix},\begin{bmatrix}\tilde{\theta}_{2}\\ \tilde{\phi}_{2}\end{bmatrix}\in\text{ker}\left(\Phi_{1}\right)\subset\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{s}}$. We will show that this leads to a contradiction by constructing a nontrivial element in the kernel of $DF(x)$.

Define $\xi_{1},\xi_{2}:[0,1]\rightarrow\mathbb{R}^{n}$ by

[TABLE]

then a straightforward computation shows that

[TABLE]

where the boundary condition at $t=1$ follows from the fact that $\begin{bmatrix}\tilde{\theta}_{1}\\ \tilde{\phi}_{1}\end{bmatrix},\begin{bmatrix}\tilde{\theta}_{2}\\ \tilde{\phi}_{2}\end{bmatrix}\in\text{ker}\left(\Phi_{1}\right)$. Further note that $\xi_{1}(t)$ and $\xi_{2}(t)$ are linearly independent for each $t\in[0,1]$, since the vectors $\begin{bmatrix}\tilde{\theta}_{1}\\ \tilde{\phi}_{1}\end{bmatrix}$ and $\begin{bmatrix}\tilde{\theta}_{2}\\ \tilde{\phi}_{2}\end{bmatrix}$ are, and the operators $D_{\theta}P\left(\theta\right)\iota_{u}$ and $D_{\phi}Q\left(\phi\right)\iota_{s}$ are injective. Consequently, since

[TABLE]

where $a^{1},\ldots,a^{m}$ are real, there exist (unique) real Chebyshev coefficients $b_{1},b_{2}\in\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}$ such that

[TABLE]

In particular, note that any linear combination of $b_{1}$ and $b_{2}$ corresponds to a solution of (3.20) and is thus an element in $\text{ker}\left(DF_{u}\left(a\right)\right)$ by Lemma 3.21.

Now, set $\tilde{\lambda}^{u}:=\mathbb{0}_{n_{u}}$, $\tilde{\lambda}^{s}:=\mathbb{0}_{n_{s}}$, $\tilde{p}:=0$, $\tilde{q}:=0$ and

[TABLE]

If $D\eta(a)b_{j}=0$ for some $j\in\{1,2\}$, where $\eta$ is the phase condition defined in (3.19), then a straightforward computation shows that $h_{j}\in\text{ker}\left(DF(x)\right)$. Otherwise, without loss of generality, we may assume that $D\eta(a)b_{1}\not=0$ and set

[TABLE]

A straightforward computation then shows that $0\neq h\in\text{ker}\left(DF(x)\right)$. Therefore, we have reached a contradiction, since $DF(x)$ is assumed to be injective. Hence $u$ must be transverse. ∎

Remark 3.23**.**

In practice, the injectivity of $DF(x)$ follows directly from our computer-assisted proof (a contraction argument) and is thus obtained for “free”, see Remark 4.10.

4 Functional analytic setup

In this section we set up a functional analytic framework for establishing the existence of an isolated zero of $F$. We start by introducing some notation and a finite dimensional reduction of $F$. We then combine numerical simulation and analysis on paper to set up a Newton-like operator $T$ whose fixed points correspond to zeros of $F$. Finally, we derive a finite number of inequalities to establish that $T$ is a contraction in a neighborhood of an approximate zero.

4.1 Projections

In this section we define projections on both the range and domain. These projections will help structure the calculations in the following sections.

Projections in $\mathcal{X}_{\nu}$

Write $x=\left(\theta,\phi,\lambda^{u},\lambda^{s},a^{1},\ldots,a^{m},p,q\right)\in\mathcal{X}_{\nu}$. Let $1\leq i\leq m$, $1\leq j\leq n$ and $\bm{k}\subset\mathbb{N}_{0}$. Define projections $\Pi_{a}^{i}:\mathcal{X}_{\nu}\rightarrow\ell^{1}_{\nu_{i},n}$,  $\Pi_{a}^{ij},\ \Pi^{ij\bm{k}}_{a}:\mathcal{X}_{\nu}\rightarrow\ell^{1}_{\nu_{i}}$ onto the Chebyshev coefficients by

[TABLE]

In particular, if $\bm{k}$ is a singleton, we identify $\Pi^{ij\bm{k}}_{a}(x)\simeq\left[a^{i}_{\bm{k}}\right]_{j}$. Similarly, let $\mathcal{I}\subset\mathbb{N}_{0}^{n_{u}}$, $\mathcal{J}\subset\mathbb{N}_{0}^{n_{s}}$ and define projections $\Pi_{P}:\mathcal{X}_{\nu}\rightarrow W^{1}_{\nu_{u,n}}$, $\Pi^{j}_{P},\ \Pi^{j\mathcal{I}}_{P}:\mathcal{X}_{\nu}\rightarrow W^{1}_{\nu_{u}}$ and $\Pi_{Q}:\mathcal{X}_{\nu}\rightarrow W^{1}_{\nu_{s,n}}$, $\Pi^{j}_{Q},\ \Pi^{j\mathcal{J}}_{Q}:\mathcal{X}_{\nu}\rightarrow W^{1}_{\nu_{s}}$ onto the Taylor coefficients of the (un)stable manifolds by

[TABLE]

As before, if $\mathcal{I}$ and $\mathcal{J}$ are singletons, we identify $\left(\Pi^{j\mathcal{I}}_{P}\left(x\right)\right)\simeq\left[p_{\mathcal{I}}\right]_{j}$, $\left(\Pi^{j\mathcal{J}}_{Q}\left(x\right)\right)\simeq\left[q_{\mathcal{J}}\right]_{j}$. Finally, we define projections $\Pi^{j}_{\theta},\Pi^{j}_{\phi},\Pi^{j}_{\lambda^{u}},\Pi^{j}_{\lambda^{s}}:\mathcal{X}_{\nu}\rightarrow\mathbb{C}$ by

[TABLE]

Remark 4.1**.**

In order to keep the notation and number of symbols to a minimum, we have used the symbol $k$ as a “dummy” index which, depending on the context, can be either an element in $\mathbb{N}_{0},\mathbb{N}_{0}^{n_{u}}$ or $\mathbb{N}_{0}^{n_{s}}$.

We shall denote the collection of projections onto the components of $\mathcal{X}_{\nu}$ by $\mathbb{P}$, i.e.,

[TABLE]

Observe that

[TABLE]

where $\left\|\cdot\right\|_{\Pi\left(\mathcal{X}_{\nu}\right)}$ denotes the norm on $\Pi\left(\mathcal{X}_{\nu}\right)\in\left\{\mathbb{C},\ell^{1}_{\nu_{i}},W^{1}_{\nu_{u}},W^{1}_{\nu_{s}}\right\}$. In the following we shall be a bit more sloppy in describing subsets of $\mathbb{P}$ by omitting the ranges for the components of the projections. For instance, whenever we write $\left\{\Pi_{\theta}^{j},\Pi_{\phi}^{j},\Pi^{j}_{\lambda^{u}},\Pi^{j}_{\lambda_{s}}\right\}$, we mean to say that this set contains all components associated to these projections, i.e., it contains $\Pi^{j}_{\theta}$ for $1\leq j\leq n_{u}$, $\Pi^{j}_{\phi}$ for $1\leq j\leq n_{s}$, etc. We will adopt the same convention for the projections into the range which are introduced at the end of this section.

Galerkin projection into $\mathcal{X}_{\nu}$

Let $N\in\mathbb{N}^{m}$ be a given truncation parameter and define operators $\Pi^{i}_{N}:\ell^{1}_{\nu_{i},n}\rightarrow\ell^{1}_{\nu_{i},n}$ by

[TABLE]

Similarly, let $K=\left(K^{u},K^{s}\right)\in\mathbb{N}^{n_{u}}\times\mathbb{N}^{n_{s}}$, set

[TABLE]

and define $\Pi_{K^{u}}:W^{1}_{\nu_{u},n}\rightarrow W^{1}_{\nu_{u},n}$, $\Pi_{K^{s}}:W^{1}_{\nu_{s},n}\rightarrow W^{1}_{\nu_{s},n}$ by

[TABLE]

respectively. Finally, define the Galerkin-projection $\Pi_{\text{dom}}^{NK}:\mathcal{X}_{\nu}\rightarrow\mathcal{X}_{\nu}$ into the domain by

[TABLE]

where $I_{2\left(n_{u}+n_{s}\right)}$ is the identity on $\mathbb{C}^{2\left(n_{u}+n_{s}\right)}$, and set $\mathcal{X}^{NK}:=\Pi^{NK}_{\text{dom}}\left(\mathcal{X}_{\nu}\right)$.

Remark 4.2**.**

In practice, we choose an ordering on the set of multi-indices and identify $\left[\Pi_{K_{u}}\left(p\right)\right]_{j}$ with a column vector in $\mathbb{C}^{\prod_{i=1}^{n_{u}}\left(K_{i}^{u}+1\right)}$ (and we do the same for the Taylor coefficients associated to the stable manifold). The Chebyshev coefficients $\Pi^{i}_{N}\left(a\right)$ are identified with a column vector in $\mathbb{C}^{nN_{i}}$ as well. Altogether, this yields an identification of $\Pi^{NK}_{\text{dom}}(x)$ with a vector in $\mathbb{C}^{\kappa}$, where

[TABLE]

since $n_{u}+n_{s}=n+1$. In particular, $\mathcal{X}^{NK}\simeq\mathbb{C}^{\kappa}$.

Projections in $\mathcal{Y}_{\tilde{\nu}}$

Recall that the range consists of elements of the form $y=\left(y_{t_{0}},y_{t_{m}},\right.$ $\left.y_{\hat{p}_{1}},y_{\hat{q}_{1}},y_{\eta},y_{a},y_{p},y_{q}\right)$, see Remark 3.17. We shall abuse notation and denote the projections onto the Chebyshev and Taylor coefficients in the range in the same way as for the domain. Furthermore, we define projections $\Pi^{j}_{t_{0}},\Pi^{j}_{t_{m}},\Pi^{j}_{\hat{p}_{1}},\Pi^{j}_{\hat{q}_{1}},\Pi_{\eta}:\mathcal{Y}_{\tilde{\nu}}\rightarrow\mathbb{C}$ by

[TABLE]

The truncation operators on the range are defined in the same way as for the domain. For this reason we shall use the same notation to denote them. There is one slight modification in the projection onto the Chebyshev coefficients associated to the first domain, however, namely we set

[TABLE]

Finally, we define the Galerkin-projection $\Pi^{NK}_{\text{ran}}$ into the range by

[TABLE]

and set $\mathcal{Y}^{NK}:=\Pi^{NK}_{\text{ran}}\left(\mathcal{Y}_{\tilde{\nu}}\right)$.

Remark 4.3**.**

Observe that $\mathcal{Y}^{NK}\simeq\mathbb{C}^{\kappa}$, since $n_{u}+n_{s}=n+1$. Hence $\mathcal{Y}^{NK}$ and $\mathcal{X}^{NK}$ have the same dimension.

4.2 An equivalent fixed-point problem

In this section we construct a Newton-like operator and set up an equivalent fixed point problem. We start by introducing a finite dimensional reduction of the zero finding problem amenable to numerical computations:

Definition 4.4** (Finite dimensional reduction).**

The finite dimensional reduction $F_{NK}:\mathcal{X}^{NK}\rightarrow\mathcal{Y}^{NK}$ of the connecting orbit map $F$ is defined by

[TABLE]

Next, we construct an approximation of $DF\left(\hat{x}\right)$ and its inverse by combining numerical computations and analysis on paper. To this end, assume that we have computed

• (A1)

an approximate zero $\hat{x}=\left(\hat{\theta},\hat{\phi},\hat{\lambda}^{u},\hat{\lambda}^{s},\hat{a},\hat{p},\hat{q}\right)\in\mathcal{X}^{NK}$ of $F_{NK}$ such that $\left(\hat{x}\right)^{\star}=\hat{x}$,

• (A2)

an approximate injective inverse $A_{NK}$ of $DF_{NK}\left(\hat{x}\right)$.

If the truncation parameters are sufficiently large and the grid is sufficiently fine, we expect the linear part of the mappings $F_{u}$, $F_{P}$ and $F_{Q}$ to be dominant in a small neighborhood of the approximate zero. This motivates the following definitions:

Definition 4.5** (Approximate derivative).**

The approximate derivative $\widehat{DF}:\mathcal{X}_{\nu}\rightarrow\mathcal{Y}_{\tilde{\nu}}$ at $\hat{x}$ is defined by

[TABLE]

where $y=DF_{NK}\left(\hat{x}\right)\Pi^{NK}_{\text{dom}}(h)$.

Definition 4.6** (Approximate inverse).**

The approximate inverse $A:\mathcal{Y}_{\nu}\rightarrow\mathcal{X}_{\nu}$ of $DF\left(\hat{x}\right)$ is defined by

[TABLE]

where $x=A_{NK}\Pi_{\text{ran}}^{NK}(h)$.

Remark 4.7**.**

Note that $A$ is injective, since $A_{NK}$ is (by assumption (A2)).

We are now ready to construct a Newton-like operator for $F$ based at the approximate zero:

Definition 4.8** (Newton-like operator).**

The Newton-like operator $T$ for $F$ based at $\hat{x}$ is defined by $T:=I-AF$.

A straightforward computation shows that $T$ maps $\mathcal{X}_{\nu}$ into itself by construction of the approximate inverse $A$. The weights $\tilde{\nu}$ are therefore irrelevant. Furthermore, observe that $T(x)=x$ if and only if $F(x)=0$, since $A$ is injective. We conclude this section with a theorem which can be used to prove that $T$ is a contraction in a neighborhood of $\hat{x}$ by checking a finite number of inequalities. The theorem is based on a parameterized Newton-Kantorovich method and is often referred to as the radii-polynomial approach (see [MR3454370] for instance).

Theorem 4.9** (Contraction mapping principle with variable radius).**

Suppose for each $\Pi\in\mathbb{P}$ there exist bounds $Y_{\Pi},Z_{\Pi}(r)>0$ such that

[TABLE]

where $\left\|\cdot\right\|_{\Pi\left(\mathcal{X}_{\nu}\right)}$ denotes the norm on $\Pi\left(\mathcal{X}_{\nu}\right)$. If there exists a radius $\hat{r}>0$ such that

[TABLE]

for all $\Pi\in\mathbb{P}$, then $T:B_{\hat{r}}\left(\hat{x}\right)\rightarrow B_{\hat{r}}\left(\hat{x}\right)$ is a contraction.

Proof.

A proof can be found in [Yamamoto]. ∎

Remark 4.10**.**

If $T$ is a contraction on $B_{\hat{r}}\left(\hat{x}\right)$, then there exists a unique zero $\tilde{x}\in B_{\hat{r}}\left(\hat{x}\right)$ of $F$. In particular, $\left(\tilde{x}\right)^{\bm{\star}}\in B_{\hat{r}}\left(\hat{x}\right)$, since $\left(\hat{x}\right)^{\bm{\star}}=\hat{x}$ and $\bm{\star}$ is norm-preserving. Therefore, since $\tilde{x}$ is unique in $B_{\hat{r}}\left(\hat{x}\right)$, we haven proven the existence of a real connecting orbit by Proposition 3.19. Moreover, $DF\left(\tilde{x}\right)$ is injective, since

[TABLE]

by (4.2) and (4.3). Hence the connecting orbit is transverse by Proposition 3.22.

Remark 4.11**.**

The radius ${\hat{r}}$ in Theorem 4.9 serves as an error bound on the solution. Indeed, since for any $\nu\geq 1$ the $\ell^{1}_{\nu}$ norm controls the $C^{0}$ norm, it follows that along the part of the heteroclinic orbit between the local invariant manifolds, which is described by the (domain decomposed) Chebyshev series, the distance in phase space between the numerical approximation and the solution is bounded by ${\hat{r}}$. An analogous bound holds for the distance between the parts of the orbit lying near the (un)stable manifolds: the distance between the orbit and the numerical approximation of the local manifold is bounded by ${\hat{r}}$. More precisely locating the position of the orbit within the local stable manifold (and analogously for the unstable one) requires solving the corresponding linear flow with initial data ${\hat{\phi}}$, see (3.2), and considering its image under the mapping $\hat{Q}:\phi\to\sum_{k\in\mathcal{K}^{s}}\hat{q}_{k}\phi^{k}$.

5 Bounds for proving contraction

In this section we compute the bounds as stated in Theorem 4.9 to prove that $T$ is a contraction in a neighborhood of $\hat{x}$. To compute these bounds, we need to project and perform analysis on the various subspaces of $\mathcal{X}_{\nu}$. Since the analysis for the unstable and stable manifold is the same, we will only write down the arguments in detail for the unstable manifold and simply state the analogous result for the stable manifold. We have aimed to compute the sharpest bounds whenever possible, but there are occasions in which we have chosen to use slightly less optimal bounds when the reduction in computational complexity outweighed the potential loss in accuracy.

5.1 $Y$-bounds

In this section we compute bounds for the residual

[TABLE]

as stated in Theorem 4.9. To this end, observe that

[TABLE]

Furthermore, for $k\in\mathbb{N}_{0}$, $1\leq i\leq m$, $1\leq j\leq n$, we have that

[TABLE]

since $\left[c_{k}\left(\hat{a}\right)\right]_{j}=0$ for all $k\geq N_{g_{j}}\left(N_{i}-1\right)+1$, see the definition of $c$ in (3.11). Similarly, for $k\in\mathbb{N}_{0}^{n_{u}}$,

[TABLE]

where

[TABLE]

since $\left[C_{k}\left(\hat{p}\right)\right]_{j}=0$ for $k\not\in\mathcal{K}^{u}$, see (3.4).

The above computations show that there are only a finite number of non-vanishing terms in $AF\left(\hat{x}\right)$. Therefore, $AF\left(\hat{x}\right)$ can be computed with the aid of a computer. It is now a straightforward task to compute the $Y$-bounds by taking the appropriate norms of the above expressions:

Proposition 5.1** ($Y$-bounds).**

Let $1\leq i\leq m$ and $1\leq j\leq n$. The bounds

[TABLE]

satisfy the estimate in (4.1).

5.2 Z-bounds

In this section we compute bounds for $DT$ as stated in Theorem 4.9. To this end, let $r>0$, $v,h\in B_{1}(0)$ be arbitrary and observe that

[TABLE]

We shall use this decomposition to compute quadratic polynomials $Z_{\Pi}(r)$ which satisfy the condition in (4.2). Furthermore, throughout this section we shall write

[TABLE]

Let us start with analyzing the easiest term which measures the quality of the approximate derivative and inverse:

Lemma 5.2**.**

Let $\Pi\in\mathbb{P}$, then

[TABLE]

Proof.

It suffices to observe that

[TABLE]

where $I_{NK}$ is the identity on $\mathcal{X}^{NK}_{\nu}$. The latter equality holds because $A$ and $\widehat{DF}$ are exact inverses of each other on the subspaces associated to the “tails” in $W^{1}_{\nu_{u},n}$, $W^{1}_{\nu_{s},n}$ and $\bigoplus_{i=1}^{m}\ell^{1}_{\nu_{i},n}$. ∎

Remark 5.3**.**

Note that the computation of the stated bound is finite for each $\Pi\in\mathbb{P}$, since $I_{NK}-A_{NK}DF_{NK}\left(\hat{x}\right)$ is a finite dimensional matrix.

5.2.1 Chebyshev series: convolution terms

In this section we develop tools for analyzing the terms

[TABLE]

which will be used extensively in Section 5.2.2 to compute the $Z$-bounds. We start with the observation that

[TABLE]

for $1\leq i\leq m$, where

[TABLE]

The goal is to construct workable matrix representations for both the finite (truncated) and tail part of (5.2.1).

To construct suitable matrix representations for (5.2.1), recall that

[TABLE]

where the coefficients $\hat{g}^{ijl}$ are defined in (3.13), see Remark 3.9. In particular, note that $\hat{g}^{ijl}_{k}=0$ for $k\geq M_{jl}\left(N_{i}-1\right)+1$, where $M_{jl}:=\text{order}\left(\dfrac{\partial g_{j}}{\partial x_{l}}\right)$, since $\hat{a}^{i}_{k}=0$ for $k\geq N_{i}$. Therefore, motivated by the above observations, we consider a sequence $a\in\ell^{1}_{\nu}$ such that $a_{k}=0$ for $k\geq\tilde{M}:=M\left(N-1\right)+1$, where $M\in\mathbb{N}$, $N\in\mathbb{N}$, $\nu>1$, and construct explicit matrix representations for the mappings $B(a),\Gamma(a):\ell^{1}_{\nu}\rightarrow\ell^{1}_{\nu}$ defined by

[TABLE]

where

[TABLE]

Remark 5.4**.**

The parameters $N$ and $\nu$ in this section are not to be confused with the vector valued ones used throughout this paper. In practice, we set $a=\hat{g}^{ijl}$, $N=N_{i}$ and $M=M_{jl}$. In particular, observe that

[TABLE]

by (5.2.1).

We begin by extending $a$ and $\tilde{a}$ to “bi-infinite” sequences, by setting $a_{-k}:=a_{k}$ and $\tilde{a}_{-k}:=\tilde{a}_{k}$ for $k\in\mathbb{N}$, and constructing a bi-infinite matrix representation for the map

[TABLE]

We then convert this bi-infinite matrix representation to an “one-sided” matrix representation by using appropriate reflections, which in turn will be used to construct the desired matrix representations for $B(a)$ and $\Gamma(a)$. To be more precise, first observe that

[TABLE]

Here we have identified elements in $\ell^{1}_{\nu}$ with bi-infinite column vectors (with respect to the ordering as depicted above). The bandwidth of this bi-infinite matrix is $\tilde{M}-1$, since $a_{k}$ vanishes for $\left|k\right|\geq\tilde{M}$. The shaded regions in grey indicate the position of the “zeroth” row and column. Set $b_{k}:=a_{k-1}-a_{k+1}$ for $-\tilde{M}\leq k\leq\tilde{M}$, then it follows from the above expression that

[TABLE]

In particular, this bi-infinite matrix has bandwidth $\tilde{M}$.

Next, we convert the latter matrix representation into a one-sided representation on $\mathbb{N}_{0}$ by “reflecting” all elements on the left hand-side of the zeroth column to the right and ignoring the rows with negative indices. This yields

[TABLE]

Altogether, the sum of the above two matrices, which we will denote by $S(a)$, constitutes an infinite dimensional matrix representation of the map

[TABLE]

Let $\mathbb{\tilde{B}}(a)\in\mathbb{C}^{(N-1)\times\left(\tilde{M}+N\right)}$ denote the finite dimensional submatrix of $S(a)$ defined by

[TABLE]

Note that we are using the convention that the indexing of the rows and columns start at zero rather than at one. In view of (5.2), set the elements in the columns of $\mathbb{\tilde{B}}(a)$ with index $0\leq j\leq N-1$ to zero and let $\mathbb{B}(a)$ denote the resulting matrix. Finally, let $\Gamma_{\infty}(a)$ denote the infinite dimensional matrix which consists of the rows of $S(a)$ with index $N$ and higher. Then

[TABLE]

by construction.

In preparation for the analysis in Section 5.2.2, we show that the operator norm of $\Gamma(a)$ can be computed by considering a sufficiently large finite dimensional submatrix of $\Gamma_{\infty}(a)$.

Lemma 5.5** (Operator norm of $\Gamma(a)$).**

Let $\bm{\Gamma}_{N}(a)\in\mathbb{C}^{\left(3\tilde{M}-N+1\right)\times\left(2\tilde{M}+1\right)}$ denote the submatrix $\left\{S(a)_{ij}:N\leq i\leq 3\tilde{M},\ 0\leq j\leq 2\tilde{M}\right\}$ of $S(a)$. Then

[TABLE]

Proof.

It follows directly from the expression in (5.21) that

[TABLE]

where $\left(\xi_{k^{\prime}}\right)_{k^{\prime}\in\mathbb{N}_{0}}$ are the corner points introduced in Definition 2.7. This implies that $\left\|\Gamma(a)\xi_{k^{\prime}}\right\|_{\nu}=\left\|\Gamma(a)\xi_{2\tilde{M}}\right\|_{\nu}$ for all $k^{\prime}\geq 2\tilde{M}$, hence

[TABLE]

by Proposition 2.9 and the definition of $\bm{\Gamma}_{N}(a)$. ∎

Remark 5.6**.**

The latter results shows that the operator norm of $\Gamma(a)$ is determined by its first $2\tilde{M}+1$ columns.

5.2.2 First order bounds

In this section we compute bounds for

[TABLE]

by projecting it onto all the relevant subspaces of $\mathcal{X}_{\nu}$. For notational convenience, we shall write $y=\left(DF\left(\hat{x}\right)-\widehat{DF}\right)h$ throughout this section. We start by computing the difference between the exact and approximate derivative.

A straightforward computation shows that

[TABLE]

Furthermore, $\Pi(y)=0$ for $\Pi\in\left\{\Pi_{\hat{p}_{1}},\Pi_{\hat{q}_{1}},\Pi_{\eta}\right\}$, since the equations associated to $\Pi F$ are linear and only depend on elements in the finite dimensional subspace $\mathcal{X}^{NK}$. Next, set

[TABLE]

where $M_{jl}=\text{order}\left(\dfrac{\partial g_{j}}{\partial x_{l}}\right)$, then it follows from (5.3), (5.4) and (5.22) that

[TABLE]

for $2\leq i\leq m$, see Section 5.2.1. Here the matrix-vector product $\mathbb{B}\left(\hat{g}^{ijl}\right)\Pi_{a}^{il\bm{k}_{ijl}}\left(h\right)$ is interpreted by using the identification

[TABLE]

The same formula holds for $i=1$ and $k\in\mathbb{N}$. In particular, if $i=1$, then there is no component to consider for $k=0$. Finally, we compute that

[TABLE]

Altogether, the above formulae give rise to the decomposition

[TABLE]

where we have set

[TABLE]

The strategy is to compute bounds for (5.23) by individually composing each term in the above decomposition with $A$ and analyzing the associated projections into the domain.

Remark 5.7**.**

Observe that $A\Pi^{\tilde{\imath}\tilde{\jmath}\left[N_{\tilde{\imath}}:\infty\right)}_{a}$, $A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{u}}_{P}$ and $A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{s}}_{Q}$, are “diagonal” and “uncoupled” in the sense that the only nonzero projections into the domain are $\Pi_{a}^{\tilde{\imath}\tilde{\jmath}}A\Pi^{\tilde{\imath}\tilde{\jmath}\left[N_{\tilde{i}}:\infty\right)}$, $\Pi^{\tilde{\jmath}}_{P}A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{u}}_{P}$ and $\Pi^{\tilde{\jmath}}_{Q}A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{s}}_{Q}$.

Boundary conditions

We start by considering the terms associated to $A\Pi^{\tilde{\jmath}}_{t_{0}}(y)$, $A\Pi^{\tilde{\jmath}}_{t_{m}}(y)$ and $A\Pi^{\tilde{\imath}\tilde{\jmath}0}(y)$, which are related to the boundary conditions.

Lemma 5.8**.**

Let $\Pi\in\mathbb{P}$, $2\leq\tilde{\imath}\leq m$ and $1\leq\tilde{\jmath}\leq n$, then

[TABLE]

where $\left\|\cdot\right\|_{\Pi\left(\mathcal{X}_{\nu}\right)}$ denotes the corresponding norm on $\Pi\left(\mathcal{X}_{\nu}\right)$.

Proof.

First observe that

[TABLE]

by Proposition 2.9. Similarly,

[TABLE]

by Proposition 2.12, where in the last line we used that $\hat{\theta}\in\mathbb{B}_{\nu_{u}}$. Finally, the above bounds and the expressions in (5.24) and (5.25) show that

[TABLE]

which proves the statement. ∎

Remark 5.9**.**

The computation of the stated bounds is finite for each $\Pi\in\mathbb{P}$, since $A_{NK}$ is a finite dimensional matrix.

Chebyshev coefficients

Next, we consider the terms associated to $A\Pi^{\tilde{\imath}\tilde{\jmath}\mathbb{N}}_{a}(y)$. We start with the observation that

[TABLE]

by (5.25). Note that $A_{NK}\Pi^{\tilde{\imath}\tilde{\jmath}\left[1:N_{\tilde{i}}-1\right]}_{a}\mathbb{B}\left(\hat{g}^{\tilde{\imath}\tilde{\jmath}l}\right)$ is a finite dimensional matrix which can be explicitly computed on a computer. In particular,

[TABLE]

corresponds to a finite dimensional matrix representation of a linear operator on $\ell^{1}_{\nu_{\tilde{\imath}}}$. Hence the computation of its operator norm is finite.

Lemma 5.10** (Scalar and Taylor projections).**

Let $1\leq\tilde{\imath}\leq m$, $1\leq\tilde{\jmath}\leq n$ and $\Pi\in\left\{\Pi^{j}_{\theta},\Pi^{j}_{\phi},\Pi_{\lambda^{u}}^{j},\Pi^{j}_{\lambda^{s}}\Pi^{j}_{P},\Pi^{j}_{Q}\right\}$, then

[TABLE]

Proof.

It suffices to observe that

[TABLE]

by construction of the approximate inverse $A$. Hence the result follows directly from (5.31). ∎

To analyze the terms $\Pi_{a}^{ij}A\Pi_{a}^{\tilde{\imath}\tilde{\jmath}\mathbb{N}}(y)$ for $1\leq i\leq m$, $1\leq j\leq n$, we first derive a more explicit expression for the tail $A\Pi_{a}^{\tilde{\imath}\tilde{\jmath}\left[N_{\tilde{\imath}}:\infty\right)}(y)$. For this purpose, define a (infinite dimensional) diagonal matrix $D^{i}_{\infty}$ by

[TABLE]

Then it follows from (5.25) and the definition of the approximate inverse that

[TABLE]

Altogether, by combining the latter result with (5.31), we conclude that if $\left(i,j\right)\not=\left(\tilde{\imath},\tilde{\jmath}\right)$, then

[TABLE]

Otherwise, if $\left(i,j\right)=\left(\tilde{\imath},\tilde{\jmath}\right)$, then

[TABLE]

We are now ready to compute the desired bounds. Before we proceed, observe that the operator norms of the infinite dimensional matrices in (5.33) can be computed by considering sufficiently large finite dimensional submatrices by the same reasoning as in Lemma 5.5. The details are given in the lemma below.

Lemma 5.11** (Chebyshev projections).**

Let $1\leq\tilde{\imath},i\leq m$, $1\leq\tilde{\jmath},j\leq n$. If $\left(\tilde{\imath},\tilde{\jmath}\right)\not=\left(i,j\right)$, then

[TABLE]

Otherwise, if $\left(\tilde{\imath},\tilde{\jmath}\right)=\left(i,j\right)$, then

[TABLE]

where

[TABLE]

and $\mathbb{\Gamma}_{N}\left(\hat{g}^{ijl}\right)$ is defined in Lemma 5.5.

Proof.

The statement for $\left(\tilde{\imath},\tilde{\jmath}\right)\not=\left(i,j\right)$ follows directly from (5.32) and the triangle inequality. To prove the result for $\left(\tilde{\imath},\tilde{\jmath}\right)=\left(i,j\right)$, consider the linear operator

[TABLE]

A similar computation as in Lemma 5.5 shows that

[TABLE]

where $\left(\xi_{k^{\prime}}\right)_{k^{\prime}\in\mathbb{N}_{0}}$ are the corner points introduced in Definition 2.7. Note that the latter quantity is decreasing for $k^{\prime}\geq 2\tilde{M}_{ijl}$. Hence, by Proposition 2.9, the operator norm of (5.34) is completely determined by its columns with index $0\leq k^{\prime}\leq 2\tilde{M}_{ijl}$. The corresponding submatrix is given by

[TABLE]

Therefore, the result now follows from (5.33) and the triangle inequality. ∎

Taylor coefficients

Finally, we consider the terms $A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{u}}_{P}(y)$ and $A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{s}}_{Q}(y)$. In particular, recall that the only nonzero projections in this case are $\Pi^{\tilde{\jmath}}_{P}A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{u}}_{P}(y)$ and $\Pi^{\tilde{\jmath}}_{Q}A\Pi^{\tilde{\jmath}\mathcal{K}_{c}^{u}}_{Q}(y)$. To study these terms, we will use the following result:

Lemma 5.12**.**

Let $M\in\mathbb{N}$ and set $\mathcal{K}^{u}_{M}:=\left\{k\in\mathbb{N}^{n_{u}}_{0}:k_{i}\leq MK^{u}_{i},\ 1\leq i\leq n_{u}\right\}$. Suppose $p\in W^{1}_{\nu_{u}}$ satisfies $p_{k}=0$ whenever $k\not\in\mathcal{K}^{u}_{M}$ and define an operator $\Lambda_{u}(p):W^{1}_{\nu_{u}}\rightarrow W^{1}_{\nu_{u}}$ by

[TABLE]

Then $\Lambda_{u}(p)$ is bounded and

[TABLE]

Proof.

It follows directly from the Banach algebra estimate that $\Lambda^{u}(p)$ is bounded. To obtain the stated bound for operator norm, first note that

[TABLE]

for any $k\in\mathbb{N}_{0}^{n_{u}}$. In particular,

[TABLE]

for any $l\in\mathbb{N}_{0}^{n_{u}}$, where $\left(\bm{\xi}_{l}\right)_{l\in\mathbb{N}_{0}^{n_{u}}}$ are the corner points introduced in Definition 2.11. Consequently,

[TABLE]

Next, observe that for any $l\in\mathbb{N}_{0}^{n_{u}}$ and $1\leq j\leq n_{u}$,

[TABLE]

where in the last line we used the fact that all the $\hat{\lambda}^{u}_{i}$ have the same sign. Therefore, the term in (5.35) is bounded by

[TABLE]

Finally, note that for any fixed $1\leq j\leq n_{u}$, the above sum is strictly decreasing for $l_{j}\geq K_{j}^{u}+1$. Hence the desired result now follows from Proposition 2.12. ∎

Remark 5.13**.**

A similar statement holds for the map associated to the stable manifold. The corresponding operator (for $q\in W^{1}_{\nu_{s}}$) is denoted by $\Lambda_{s}\left(q\right):W^{1}_{\nu_{s}}\rightarrow W^{1}_{\nu_{s}}$.

We are now ready to compute the required bounds. To this end, observe that

[TABLE]

by the reasoning in Remark 3.9, where $\left[C\right]_{j}$ denotes the $j$-th component of $C$ and $\hat{G}_{u}^{jl}$ are the Taylor coefficients of

[TABLE]

The coefficients $\hat{G}_{s}^{jl}$ associated to the stable manifold are defined similarly.

Lemma 5.14** (Projection onto the Taylor coefficients).**

Let $1\leq j\leq n$, then

[TABLE]

Proof.

It suffices to observe that

[TABLE]

by (5.26), (5.36) and the definition of the approximate inverse. ∎

The norms in the right-hand sides of (5.37) and (5.38) are controlled by using Lemma 5.12.

First order coefficients of $Z_{\Pi}(r)$

We are now ready to construct the first order terms of the (quadratic) polynomials $Z_{\Pi}(r)$ for $\Pi\in\mathbb{P}$. For this purpose, we first introduce some additional notation. We will denote the bounds in Lemma 5.2, which measure the quality of the approximate derivative and inverse, by $H_{\Pi}$. The bounds in (5.28), (5.29) and (5.30), which are related to the boundary conditions, will be denoted by $Z^{1,\tilde{\jmath}t_{0}}_{\Pi}$, $Z^{1,\tilde{\jmath}t_{m}}_{\Pi}$ and $Z^{1,\tilde{\imath}\tilde{\jmath}0}_{\Pi}$, respectively. The bounds in Lemmas 5.10 and 5.11, which are related to the differential equation, will be denoted by $Z^{1,\tilde{\imath}\tilde{\jmath}\left[1:N_{\tilde{\imath}-1}\right]}_{\Pi}$ and $Z^{1,\tilde{\imath}\tilde{\jmath}\mathbb{N}}_{\Pi^{ij}_{a}}$, respectively. Finally, the bounds in (5.37) and (5.38), which are related to the invariance equation for the charts on the (un)stable manifolds, will be denoted by $Z^{1,j}_{P}$ and $Z^{1,j}_{Q}$, respectively.

With the above notation in place, we define

[TABLE]

for $\Pi\in\left\{\Pi_{\theta}^{j},\Pi_{\phi}^{j},\Pi^{j}_{\lambda^{u}},\Pi^{j}_{\lambda_{s}}\right\}$, and

[TABLE]

5.2.3 Second order bounds

In this section we compute bounds for

[TABLE]

We will compute the desired bounds by projecting (5.39) onto the relevant subspaces of $\mathcal{X}_{\nu}$ as in the previous section. We start with the observation that

[TABLE]

by the (generalized) Mean Value Theorem. For notational convenience, we shall write $y\left(\tau\right)=D^{2}F\left(\hat{x}+\tau rv\right)\left[v,h\right]$ and $v=\left(\dot{\theta},\dot{\phi},\dot{\lambda}^{u},\dot{\lambda}^{s},\dot{a},\dot{p},\dot{q}\right)$. Furthermore, we will denote the max-norm on both $\mathbb{C}^{n_{u}}$ and $\mathbb{C}^{n_{s}}$ by $\left\|\cdot\right\|_{\infty}$.

Observe that $\Pi\left(y\left(\tau\right)\right)=0$ for $\Pi\in\left\{\Pi_{\hat{p}_{1}},\Pi_{\hat{q}_{1}},\Pi_{\eta}\right\}$, since the equations associated to these projections are linear. Furthermore, a straightforward computation shows that

[TABLE]

for $1\leq\tilde{\jmath}\leq n$. An analogous formula holds for $\Pi^{\tilde{\jmath}}_{t_{m}}\left(y\left(\tau\right)\right)$. Next, let $2\leq i\leq m$, then

[TABLE]

This formula is also valid for $i=1$ and $k\in\mathbb{N}$, but in this case there is no component to consider for $k=0$. Finally, observe that

[TABLE]

The formula for $\Pi_{Q}\left(y\left(\tau\right)\right)$ is analogous.

Altogether, the above formulae give rise to the decomposition

[TABLE]

We will now follow the same strategy as in the previous section and compute bounds for (5.39) by individually composing each term in the above decomposition with $A$ and analyzing the associated projections into the domain.

Boundary conditions

We start by considering the terms associated to

[TABLE]

which are related to the boundary conditions. To compute the desired bounds, we first analyze the two series in (5.2.3).

Lemma 5.15**.**

Suppose $\theta\in\mbox{int}\ \mathbb{B}_{\nu_{u}}\setminus\{0\}$ and let $\zeta\in\mathbb{C}^{n_{u}}$ be such that $\left\|\zeta\right\|_{\infty}\leq 1$. Define a linear map $\varphi_{u,1}:W^{1}_{\nu_{u}}\rightarrow\mathbb{C}$ by

[TABLE]

Then $\varphi_{u,1}\in\left(W^{1}_{\nu_{u}}\right)^{\ast}$ and

[TABLE]

Proof.

It is clear that $\varphi_{u,1}$ is bounded for any $\theta\in\mbox{int}\ \mathbb{B}_{\nu_{u}}$, since

[TABLE]

We will use this observation to compute a bound for the operator norm. Namely, let $j\in\mathbb{N}_{0}^{n_{u}}$ be arbitrary, then

[TABLE]

where $\left(\bm{\xi}_{l}\right)_{l\in\mathbb{N}_{0}^{n_{u}}}$ are the corner points introduced in Definition 2.11. In particular, observe that $\left\|\varphi_{u,1}\right\|_{\mathcal{B}\left(W^{1}_{\nu_{u}},\mathbb{C}\right)}=\sup_{j\in\mathbb{N}_{0}^{n_{u}}\setminus{\{0\}}}\left|\varphi_{u,1}\left(\bm{\xi}_{j}\right)\right|$ by Proposition 2.12, since $\varphi_{u,1}\left(\bm{\xi}_{0}\right)=0$. Next, note that the map $x\mapsto x\rho^{x}$, where $\rho\in(0,1)$, is strictly increasing on $\left[0,\dfrac{-1}{\log\rho}\right]$, strictly decreasing on $\left[\dfrac{-1}{\log\rho},\infty\right)$ and has a global maximum on $\left[0,\infty\right)$ at $x=\dfrac{-1}{\log{\rho}}$. Therefore, since $\left|j\right|\geq 1$ for $j\in\mathbb{N}_{0}^{n_{u}}\setminus{\{0\}}$ and $\left\|\theta\right\|_{\infty}<\nu_{u}$, it follows that $\sup_{j\in\mathbb{N}_{0}^{n_{u}}\setminus{\{0\}}}\left|\varphi_{u,1}\left(\bm{\xi}_{j}\right)\right|\leq\Phi_{u,1}\left(\theta\right),$ which proves the result. ∎

Remark 5.16**.**

The analogues of $\varphi_{u,1}$ and $\Phi_{u,1}$ in the context of the stable manifold are defined similarly and are denoted by $\varphi_{s,1}$ and $\Phi_{s,1}\left(\phi\right)$, respectively.

Remark 5.17**.**

In practice, $\theta$ is an interval enclosure of $\hat{\theta}$ (see Lemma 5.21). Recall that $\hat{\theta}$ corresponds to a (numerically obtained) coordinate on the chart of the unstable manifold at which the connecting orbit $u$ starts. This is why we may assume that $\theta\not=0$.

The latter result provides a way to bound the first term in (5.2.3). To bound the second term (5.2.3), we perform a similar analysis.

Lemma 5.18**.**

Suppose $\theta\in\mbox{int}\ \mathbb{B}_{\nu_{u}}\setminus\{0\}$ and define a linear map $\varphi_{u,2}:W^{1}_{\nu_{u}}\rightarrow\mathbb{C}$ by

[TABLE]

Then $\varphi_{u,2}\in\left(W^{1}_{\nu_{u}}\right)^{\ast}$ and $\left\|\varphi_{u,2}\right\|_{\mathcal{B}\left(W^{1}_{\nu_{u}},\mathbb{C}\right)}\leq\Phi_{u,2}\left(\theta\right),$ where

[TABLE]

and $K^{u}_{\text{min}}:=\min_{1\leq i\leq n_{u}}K^{u}_{i}$.

Proof.

The boundedness of $\varphi_{u,2}$ follows directly from the observation in (5.45) and the assumption that $\theta\in\mbox{int}\ \mathbb{B}_{\nu_{u}}$. To compute a bound for the operator norm, first observe that $\varphi_{u,2}\left(\bm{\xi}_{j}\right)=0$ for $j\in\mathcal{K}^{u}$. Hence $\left\|\varphi_{u,2}\right\|_{\mathcal{B}\left(W^{1}_{\nu_{u}},\mathbb{C}\right)}=\sup_{j\in\mathcal{K}^{u}_{c}}\left|\varphi_{u,2}\left(\bm{\xi}_{j}\right)\right|$ by Proposition 2.12. Furthermore, a straightforward computation shows that

[TABLE]

for any $j\in\mathcal{K}^{u}_{c}$, where we used that $\left\|\tilde{\theta}\right\|_{\infty},\left\|\dot{\theta}\right\|_{\infty}\leq 1$. Now, note that the map $x\mapsto x^{2}\rho^{x}$, where $\rho\in(0,1)$, is strictly increasing on $\left[0,\dfrac{-2}{\log{\rho}}\right]$, strictly decreasing on $\left[\dfrac{-2}{\log{\rho}},\infty\right)$ and has a global maximum on $\left[0,\infty\right)$ at $x=\dfrac{-2}{\log{\rho}}$. Therefore, since $\left|j\right|\geq K^{u}_{\text{min}}+1$ for $j\in\mathcal{K}^{u}_{c}$, it follows that $\sup_{j\in\mathcal{K}^{u}_{c}}\left|\varphi_{u,2}\left(\bm{\xi}_{j}\right)\right|\leq\Phi_{u,2}\left(\theta\right),$ which proves the result. ∎

Remark 5.19**.**

As before, the analogues of $\varphi_{u,2}$ and $\Phi_{u,2}$ in the context of the stable manifold are denoted by $\varphi_{s,2}$ and $\Phi_{s,2}\left(\phi\right)$, respectively.

We are now ready to compute bounds for (5.44). The computation of these bounds consists of a mixture of interval analysis on the computer and ordinary estimates derived with “pen and paper”. More precisely, as mentioned at the beginning of this paper, in order to manage the rounding errors on the computer, all the bounds in this paper are computed with interval arithmetic. Roughly speaking, this means that all the elementary operations on floating point numbers are replaced by operations on intervals with endpoints representable on a computer. In this way, one can compute rigorous bounds (and hence verify inequalities) with the aid of a computer.

Now, in order to compute bounds for (5.44), we first define interval enclosures for $\hat{\theta}$, $\hat{p}$, $\hat{\phi}$ and $\hat{q}$. Let $r^{\ast}>0$ be an upper bound for the radius $r$ and set

[TABLE]

We require that $r^{\ast}$ is sufficiently small and that $\nu_{u},\nu_{s}$ are sufficiently large so that $\bm{\hat{\theta}}\subset\mbox{int}\ \mathbb{B}_{\nu_{u}}$ and $\bm{\hat{\phi}}\subset\mbox{int}\ \mathbb{B}_{\nu_{s}}$.

Remark 5.20**.**

Strictly speaking, the endpoints of the above intervals should be floating point numbers so that we can perform rigorous computations on a computer. In practice, this amounts to computing slightly larger interval enclosures for $\hat{\theta}$, $\hat{p}$, $\hat{\phi}$ and $\hat{q}$ (compared to the ones above). To avoid clutter in the notation, however, we have chosen to ignore this rather technical (but easily solved) issue.

Lemma 5.21**.**

Let $1\leq\tilde{\jmath}\leq n$, $0<r\leq r^{\ast}$ and $\Pi\in\mathbb{P}$, then

[TABLE]

where

[TABLE]

Proof.

Let $\tau\in[0,1]$ and $1\leq\tilde{\jmath}\leq n$ be arbitrary and use Lemma 5.15 to see that the first term in (5.2.3) is bounded by $2\Phi_{u,1}\left(\hat{\theta}+\tau r\dot{\theta}\right)$. To bound the second term in (5.2.3), we first split it into two series; one over $\mathcal{K}^{u}$ and one over $\mathcal{K}^{u}_{c}$. We then use Lemma 5.18 and the fact that $\hat{p}_{k}=0$ for $k\in\mathcal{K}^{u}_{c}$ to estimate

[TABLE]

Altogether, we conclude that

[TABLE]

where we used the assumption that $r\leq r^{\ast}$. Finally, observe that $\hat{\theta}+\tau r\dot{\theta}\in\bm{\hat{\theta}}$ and $\hat{p}_{k}+\tau r\dot{p}_{k}\in\bm{\hat{p}}_{k}$, since $\|\dot{\theta}\|_{\infty},\ \left\|\dot{p}\right\|_{W^{1}_{\nu_{u},n}}\leq 1$. Hence (5.48) is contained in $\bm{\beta}^{\tilde{\jmath}}_{u}$ for all $\tau\in[0,1]$, which proves the result. ∎

Remark 5.22**.**

Observe that the computation of the bound in this lemma is finite, since $\mathcal{K}^{u}$ is a finite set of multi-indices.

Chebyshev coefficients

Next, we consider the terms associated to

[TABLE]

which are related to the equations for the Chebyshev coefficients. We start by computing a bound for (5.42). Since we will have to perform a similar analysis for the Taylor coefficients in the next paragraph, we first state a general result.

Lemma 5.23**.**

Suppose $\left(X,\ast\right)$ is a Banach algebra. Let $1\leq j\leq n$ and define $G_{j}:\bigoplus_{l=1}^{n}X\rightarrow X$ and $\tilde{g}_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ by

[TABLE]

where $\left\{g_{\alpha}:\alpha\in\mathcal{A}\right\}\subset\mathbb{R}^{n}$ are the coefficients of $g$ in the monomial basis. Then

[TABLE]

for any $x=\left(x_{l}\right)_{l=1}^{n},y=\left(y_{l}\right)_{l=1}^{n},z=\left(z_{l}\right)_{l=1}^{n}\in\bigoplus_{l=1}^{n}X$ such that $\left\|y\right\|,\left\|z\right\|\leq 1$, $\tau\in[0,1]$ and $0<r\leq r^{\ast}$. Here $\left\|\cdot\right\|$ denotes the max-norm on $\bigoplus_{l=1}^{n}X$, and $\mathbb{1}_{n}$ denotes the vector of length $n$ containing $1$-s.

Proof.

Note that we may use the “usual” rules of calculus to differentiate polynomials on $\bigoplus_{l=1}^{n}X$, since $\left(X,\ast\right)$ is a Banach algebra. In particular, a straightforward computation shows that

[TABLE]

Hence, by the Banach algebra estimate,

[TABLE]

since $\left\|y\right\|,\left\|z\right\|\leq 1$, $\tau\in[0,1]$ and $0<r\leq r^{\ast}$. ∎

Remark 5.24**.**

Note that the convolution mappings $\left[C\right]_{j}$ and $\left[c\right]_{j}$ are of the form $G_{j}$.

Next, we use the above result to compute bounds for (5.49). The key observation is stated in the next lemma.

Lemma 5.25**.**

Let $1\leq\tilde{\imath}\leq m$, $1\leq\tilde{\jmath}\leq n$ and $0<r\leq r^{\ast}$, then

[TABLE]

Proof.

Define linear operators $\sigma_{1}:\ell^{1}_{\nu_{1}}\rightarrow\ell^{1}_{\nu_{1}}/\mathbb{C}$ and $\left\{\sigma_{\tilde{\imath}}:\ell^{1}_{\nu_{\tilde{\imath}}}\rightarrow\ell^{1}_{\nu_{\tilde{\imath}}}\right\}_{\tilde{\imath}=2}^{m}$ by

[TABLE]

A direct application of Proposition 2.9 shows that these operators are bounded and that

[TABLE]

Therefore, since

[TABLE]

by (5.42), the result follows directly from (5.50) and Lemma 5.23. ∎

It is now a straightforward task to compute bounds for (5.49).

Lemma 5.26** (Scalar and Taylor projections).**

Let $1\leq\tilde{\imath}\leq m$, $1\leq\tilde{\jmath}\leq n$, $0<r\leq r^{\ast}$ and $\Pi\in\left\{\Pi^{j}_{\theta},\Pi^{j}_{\phi},\Pi_{\lambda^{u}}^{j},\Pi^{j}_{\lambda^{s}}\Pi^{j}_{P},\Pi^{j}_{Q}\right\}$, then

[TABLE]

Proof.

It suffices to observe that

[TABLE]

by construction of the approximate inverse $A$. Therefore, the result follows directly from Lemma 5.25. ∎

Lemma 5.27** (Chebyshev projections).**

Let $1\leq\tilde{\imath},i\leq m$, $1\leq\tilde{\jmath},j\leq n$ and $0<r\leq r^{\ast}$, then

[TABLE]

Proof.

It follows from the definition of the approximate inverse that

[TABLE]

Hence

[TABLE]

by Proposition 2.9. The result now follows from Lemma 5.25. ∎

Taylor coefficients

Finally, we consider the terms associated to

[TABLE]

which are related to the equations for the Taylor coefficients of the (un)stable manifolds. Observe that $\Pi^{\tilde{\jmath}}_{P}\left(y\left(\tau\right)\right)\not\in W^{1}_{\nu_{u}}$, due to the presence of the terms $\left\langle\tilde{\lambda}^{u},k\right\rangle\dot{p}_{k}+\left\langle\dot{\lambda}^{u},k\right\rangle\tilde{p}_{k}$, see (5.43). For this reason, in order to facilitate the analysis of (5.51), we introduce the term $y^{\tilde{\jmath}}_{P}\left(\tau\right):=\Pi^{\tilde{\jmath}}_{P}\left(y\left(\tau\right)\right)-y^{\tilde{\jmath}}_{P,\lambda}\left(\tau\right)$, where

[TABLE]

i.e., $y^{\tilde{\jmath}}_{P}\left(\tau\right)$ is defined by removing the linear terms $\left\langle\tilde{\lambda}^{u},k\right\rangle\left[\dot{p}_{k}\right]_{\tilde{\jmath}}+\left\langle\dot{\lambda}^{u},k\right\rangle\left[\tilde{p}_{k}\right]_{\tilde{\jmath}}$ from the tail of $\Pi^{\tilde{\jmath}}_{P}\left(y\left(\tau\right)\right)$. Hence $y^{\tilde{\jmath}}_{P}\left(\tau\right)\in W^{1}_{\nu_{u}}$. An analogous decomposition $\Pi^{\tilde{\jmath}}_{Q}\left(y\left(\tau\right)\right)=y^{\tilde{\jmath}}_{Q}\left(\tau\right)+y^{\tilde{\jmath}}_{Q,\lambda}\left(\tau\right)$ is used to analyze the terms associated to the stable manifold.

Lemma 5.28**.**

Let $1\leq\tilde{\jmath}\leq n$ and $0<r\leq r^{\ast}$, then

[TABLE]

where

[TABLE]

The bounds $\left\{\bm{\sigma}^{\tilde{\jmath}}_{s}:1\leq\tilde{\jmath}\leq n\right\}$ associated to the stable manifold are defined analogously.

Proof.

Let $\tau\in[0,1]$ and $1\leq\tilde{\jmath}\leq n$ be arbitrary. It follows directly from (5.43) that

[TABLE]

since $\left\|\dot{p}\right\|_{\nu_{u},n},\left\|\tilde{p}\right\|_{\nu_{u},n}\leq 1$. Next, we consider the first order components of $y^{\tilde{\jmath}}_{P}\left(\tau\right)$ by separately analyzing the terms in (5.43) for $\left|k\right|=1$. To this end, observe that

[TABLE]

and

[TABLE]

where we used that $r\leq r^{\ast}$, $\left\|\tilde{\lambda}^{u}\right\|_{\infty},\left\|\dot{\lambda}^{u}\right\|_{\infty},\left\|\dot{p}\right\|_{\nu_{u}.n},\left\|\tilde{p}\right\|_{\nu_{u},n}\leq 1$. It now follows from the expression in (5.43) that

[TABLE]

Finally, we consider the remainder of $y^{\tilde{\jmath}}_{P}$, i.e., the components associated to $\left|k\right|\geq 2$. First, note that

[TABLE]

since $\left\|\tilde{\lambda}^{u}\right\|_{\infty},\left\|\dot{\lambda}^{u}\right\|_{\infty},\left\|\dot{p}\right\|_{\nu_{u}.n},\left\|\tilde{p}\right\|_{\nu_{u},n}\leq 1$ and $\left|k\right|\leq\left|K^{u}\right|$ for $k\in\mathcal{K}^{u}$. Furthermore,

[TABLE]

by Lemma 5.23. Therefore,

[TABLE]

Altogether, the sum of (5.52), (5.53) and (5.54) yields an upper bound for $\left\|y^{\tilde{\jmath}}_{P}\left(\tau\right)\right\|_{\nu_{u}}$, for any $\tau\in[0,1]$ and $r\leq r^{\ast}$, and is contained in $\bm{\sigma}^{\tilde{\jmath}}_{u}$. This proves the result. ∎

It is now a straightforward task to compute bounds for (5.51).

Lemma 5.29**.**

Let $1\leq\tilde{\jmath}\leq n$ and $0<r\leq r^{\ast}$, then

[TABLE]

for $\Pi\in\left\{\Pi_{\theta}^{j},\Pi_{\phi}^{j},\Pi_{\lambda^{u}}^{j},\Pi_{\lambda^{s}}^{j},\Pi^{ij}_{a},\Pi^{j}_{Q}\right\}$ and

[TABLE]

for $\Pi\in\left\{\Pi_{\theta}^{j},\Pi_{\phi}^{j},\Pi_{\lambda^{u}}^{j},\Pi_{\lambda^{s}}^{j},\Pi^{ij}_{a},\Pi^{j}_{P}\right\}$.

Proof.

It suffices to observe that

[TABLE]

by construction of the approximate inverse $A$. Therefore, since

[TABLE]

the result follows from Lemma 5.28. ∎

Lemma 5.30**.**

Let $1\leq\tilde{\jmath},j\leq n$ and $0<r\leq r^{\ast}$, then

[TABLE]

The statement and corresponding bound for the stable manifold is analogous.

Proof.

If $\tilde{\jmath}\not=j$, then

[TABLE]

by definition of the approximate inverse and $y^{\tilde{\jmath}}_{P}\left(\tau\right)$. Therefore, in this case, the result follows immediately from Lemma 5.28. Now, assume that $\tilde{\jmath}=j$ and observe that

[TABLE]

In particular, since $y^{j}_{P}\left(\tau\right)\in W^{1}_{\nu_{u}}$, we may bound the first term in the above expression by

[TABLE]

where in the last line we used Lemma 5.28 again.

Next, we derive a more explicit expression for $\left\|\Pi^{j}_{P}A\Pi^{j}_{P}\right\|_{\mathcal{B}\left(W^{1}_{\nu_{u}},W^{1}_{\nu_{u}}\right)}$. To this end, recall the definition of $\Pi^{j}_{P}A\Pi^{j}_{P}$ (see Definition 4.6) and observe that

[TABLE]

for any $k\in\mathcal{K}^{u}_{c}$. Therefore, it follows from Proposition 2.12 that

[TABLE]

Finally, a straightforward computation shows that

[TABLE]

since $\left|\left\langle\lambda^{u},k\right\rangle\right|\leq\left\|\lambda^{u}\right\|_{\infty}\left|k\right|$ and $\left|\left\langle\hat{\lambda}^{u},k\right\rangle\right|\geq\min_{1\leq i\leq n_{u}}\left|\text{Re}\left(\hat{\lambda}^{u}_{i}\right)\right|\left|k\right|$ for any $k\in\mathbb{N}_{0}^{n_{u}}$ and $\lambda^{u}\in\mathbb{C}^{n_{u}}$. Altogether, this proves the result. ∎

Second order coefficients of $Z_{\Pi}(r)$

We are now ready to finish the construction of the quadratic polynomials $Z_{\Pi}(r)$ for $\Pi\in\mathbb{P}$. As before, we first introduce some additional notation. We will denote the bounds in (5.46), (5.47), (5.55), (5.56) and Lemmas 5.26, 5.27, 5.30 by $Z^{2,\tilde{\jmath}t_{0}}_{\Pi}$, $Z^{2,\tilde{\jmath}t_{m}}_{\Pi}$, $Z^{2,\tilde{\jmath}P}_{\Pi}$, $Z^{2,\tilde{\jmath}Q}_{\Pi}$, $Z^{2,\tilde{\imath}\tilde{\jmath}\mathbb{N}}_{\Pi}$, $Z^{2,\tilde{\imath}\tilde{\jmath}\mathbb{N}}_{\Pi^{ij}_{a}}$ and $Z^{2,\tilde{\jmath}P}_{\Pi^{j}_{P}}$, $Z^{2,\tilde{\jmath}Q}_{\Pi^{j}_{Q}}$, respectively. Finally, we set

[TABLE]

and define

[TABLE]

Then $Z_{\Pi}(r)$ satisfies (4.2) by construction.

6 Applications: traveling fronts in parabolic PDEs

In this section we use our method to prove the existence of connecting orbits in systems of ODEs which arise from the study of traveling fronts in scalar parabolic PDEs. In addition, we perform discrete continuation (discrete in the sense that we rigorously validate the solution for many parameter values, but we do not attempt to obtain a continuous parametrized branch of solutions). This also demonstrates the effectiveness of the phase condition introduced in Definition 3.14. Before we proceed to the applications, we first give a rough outline of our main procedure for validating connecting orbits. We have tried to automate as many steps as possible, but there are still certain steps which are based on experimentation.

Step 1: Compute parameterizations of the local (un)stable manifolds

• $1.1$

Compute numerical approximations $\tilde{p}_{0}$ and $\tilde{q}_{0}$ of the equilibria of interest.

• $1.2$

Compute numerical approximations $\left\{\left(\tilde{\lambda}_{k}^{u},\tilde{p}_{k}\right):\left|k\right|=1\right\}$ and $\left\{\left(\tilde{\lambda}_{k}^{s},\tilde{q}_{k}\right):\left|k\right|=1\right\}$ of the eigendata associated to $Dg\left(\tilde{p}_{0}\right)$ and $Dg\left(\tilde{q}_{0}\right)$, respectively. In this step we set the length of the approximate eigenvectors to one.

• $1.3$

Choose the number of Taylor coefficients $K^{u}\in\mathbb{N}_{0}^{n_{u}}$ and $K^{s}\in\mathbb{N}_{0}^{n_{s}}$ and compute approximate zeros $\left(\hat{\lambda}^{u},\hat{p}\right)$, $\left(\hat{\lambda}^{s},\hat{q}\right)$ of the mappings

[TABLE]

by Newton’s method.

• $1.4$

If necessary, increase the truncation parameters and rescale the eigenvectors so that validation is feasible, see Remark 6.1 below.

Step 2: Compute an accurate approximation of a connecting orbit

• $2.1$

Compute a numerical approximation of a connecting orbit. This step is based on solving the truly nonlinear part of the problem and involves experimentation. It is obviously problem dependent.

• $2.2$

Use the domain decomposition algorithm developed in [domaindecomposition] to compute a grid $\left(t_{i}\right)_{i=0}^{m}$ and an accurate approximate connecting orbit

[TABLE]

so that the decay rates of the Chebyshev coefficients $\hat{a}^{i}$ are equidistributed over the subdomains $\left[t_{i-1},t_{i}\right]$. The number of modes $N_{i}$ is chosen in such a way that $\left|\hat{a}^{i}_{N_{i}-1}\right|\approx 10^{-16}$. The number of subdomains $m$ is determined by experimentation. In general, we use as many subdomains as necessary in order to ensure high decay rates of the Chebyshev coefficients.

• $2.3$

Use $\hat{u}$ as a reference orbit to fix the time parameterization of the connecting orbit (see Definition 3.14).

Step 3: Validate the connecting orbit and (un)stable manifolds

• $3.1$

Combine the results from the previous two steps to construct a symmetric approximate zero $\hat{x}=\left(\hat{\theta},\hat{\phi},\hat{\lambda}^{u},\hat{\lambda}^{s},\hat{a},\hat{p},\hat{q}\right)$ of $F_{NK}$ (see Remark 3.20).

• $3.2$

Set $r^{\ast}=10^{-5}$, $\nu_{u}=\nu_{s}=1$ (see Remark 6.1) and compute the weights $\left(\nu_{i}\right)_{i=1}^{m}$ as explained in Remark 6.2 below.

• $3.3$

Initialize the numerical data with interval arithmetic and construct the $n(m+4)+2$ radii-polynomials

[TABLE]

• $3.4$

Determine an interval $I$ on which all the radii polynomials are negative.

If we fail to find an interval $I$ on which all the radii polynomials are negative, we try to determine which parameters (the truncation parameters, the weights $\nu_{i}$ or the “scalings” of the coefficients $\hat{p}$ and $\hat{q}$) need to be modified by “visual” inspection and try again.

Remark 6.1** (Scaling and the number of Taylor coefficients).**

Observe that for any $\nu_{u},\tilde{\nu}_{u}>0$ it holds that $p\in W^{1}_{\nu_{u}}$ if and only if $\frac{\nu_{u}}{\tilde{\nu}_{u}}p\in W^{1}_{\tilde{\nu}_{u}}$ with the scaling notation introduced in Remark 3.3. Therefore, since the parameterization mappings $F_{P}$ and $F_{Q}$ are invariant under the rescaling $p\mapsto\mu p$ (see Remark 3.3), we have chosen to set $\nu_{u}=\nu_{s}=1$ and search for appropriate scalings which ensure that $\hat{p}$ and $\hat{q}$ decay sufficiently fast to zero. To be more precise, we explain in detail how we choose the scalings of the eigenvectors and the number of Taylor coefficients for the unstable manifold (the procedure for the stable manifold is analogous).

The main idea is to choose the scalings and number of Taylor coefficients in such a way that the bound for $\Pi_{P}^{j}\left(DF\left(\hat{x}\right)-\widehat{DF}\right)$ is below some prescribed tolerance. More precisely, in light of (5.26), (5.36) and Lemma 5.12, we aim to find a truncation parameter $K^{u}\in\mathbb{N}_{0}^{n_{u}}$ and a scaling factor $\mu\in\left(0,\infty\right)^{n_{u}}$ such that

[TABLE]

where $\varepsilon_{u}>0$ (in practice we set $\varepsilon_{u}=\frac{1}{2}$). We start by determining $K^{u}\in\mathbb{N}_{0}^{n_{u}}$. To this end, observe that the scaling factor $\mu$ has no effect on $\left|\hat{G}^{jl}_{0}\right|$. For this reason, we set $K^{u}:=\left(\max_{1\leq j,l\leq n_{u}}K^{u}_{jl}\right)\mathbb{1}_{n_{u}}$, where $K^{u}_{jl}\in\mathbb{N}$ is the smallest integer such that

[TABLE]

Here $\zeta\in(0,1]$ is an additional parameter chosen through experimentation (in practice we use $\zeta=\frac{3}{4}$).

Next, we determine an appropriate scaling factor $\mu$. Let $1\leq j,l\leq n$ and approximate $\left|\hat{G}^{jl}_{k}\right|\approx\left|\hat{G}^{jl}_{0}\right|\rho_{jl}^{-\left|k\right|}$, where $\rho_{jl}=e^{-s_{jl}}$ and $s_{jl}$ is the slope of the best line through the points

[TABLE]

Recall that $M_{jl}=\text{order}\left(\dfrac{\partial g_{j}}{\partial x_{l}}\right)$. Now, if $\mu_{\imath}<\rho_{jl}$ for all $1\leq\imath\leq n_{u}$ and $\left|K^{u}\right|$ is sufficiently large, then

[TABLE]

Motivated by this observation and the inequality in (6.1), we set $\mu_{\imath}=\mu$ for $1\leq\imath\leq n_{u}$ and require that

[TABLE]

for all $1\leq j,l\leq n$. Therefore, we set

[TABLE]

We remark that one could determine a more “refined” scaling factor $\mu$, which need not be the same in each direction, by taking the decay rates of $\hat{G}^{jl}$ in each separate direction into account (as opposed to using the “uniform” rate in (6.2) which ignores the different directions of the array). In addition, one could take the different sizes of the eigenvalues into account in the definition of $\xi_{jl}$.

Remark 6.2**.**

To determine the weights $\left(\nu_{i}\right)_{i=1}^{m}$, we use a heuristic procedure slightly more refined than the one used in [domaindecomposition]. Namely, we try to ensure that the bound for the tail of $\Pi_{a}^{ij}\left(DF\left(\hat{x}\right)-\widehat{DF}\right)$ is below some prescribed tolerance (rather than requiring the residual to be below some tolerance as in [domaindecomposition]). More precisely, in light of (5.25), we require that

[TABLE]

where $\varepsilon>0$ is some prescribed tolerance (in practice we set $\varepsilon=\frac{1}{2}$).

We use the rough approximation $\left|\hat{g}_{k}^{ijl}\right|\approx\left|\hat{g}^{ijl}_{0}\right|\rho^{-k}_{ijl}$, where $\rho_{ijl}=e^{-s_{ijl}}$ and $s_{ijl}$ is the slope of the best line through the points

[TABLE]

In practice, $\rho_{ijl}$ is roughly the same for all $1\leq i\leq m$ and $1\leq j,l\leq n$ due to the choice of the grid, and we therefore write $\rho=\rho_{ijl}$. If $\nu_{i}<\rho$ and $N_{i}$ is sufficiently large, then

[TABLE]

Altogether, (6.3) and (6.4) yield the constraint

[TABLE]

Finally, we determine an interval $\left[\nu_{\min},\nu_{\max}\right]\subset\mathbb{R}_{>0}$ on which (6.5) is satisfied for all $1\leq i\leq m$ and $1\leq j,l\leq n$. Then, if $\nu_{\max}>1$, we choose a weight $\hat{\nu}\in\left[\nu_{\min},\nu_{\max}\right]$ such that $1<\hat{\nu}<\rho$ and set $\nu_{i}=\hat{\nu}$ on each subdomain (in practice we set $\hat{\nu}=\frac{1}{2}\left(\max\left\{1,\nu_{\min}\right\}+\nu_{\max}\right)$). If $\nu_{\max}\leq 1$, we increase the number of subdomains (to increase $\rho$) or use a higher number of Chebyshev coefficients $N_{i}$ and then try again.

6.1 Lotka-Volterra

We have proven the existence of connecting orbits from $\left(b,0,1-b,0\right)$ to $\left(1,0,0,0\right)$ in (1.3) for $a=5,\ b=\frac{1}{2}$, $D=3$ and different values of $\kappa$. Recall that these orbits correspond to traveling fronts of (1.2) with wave speed $\kappa$. The choices for these parameter values were somewhat arbitrary and were obtained by experimenting with the parameter values considered in [MR760975]. In particular, we chose the parameters in such a way that the stable eigenvalues associated to $(1,0,0,0)$ consisted of one complex conjugate pair of eigenvalues and one real eigenvalue.

Connecting orbit at $\kappa=-1$

We started with a numerical approximation of a connecting orbit at $\kappa=-1$ and used the steps outlined in the previous section to obtain the following computational parameters:

• •

Parameterization mappings: we used $K^{u}=\begin{bmatrix}13&13\end{bmatrix}$ and $K^{s}=\begin{bmatrix}9&9&9\end{bmatrix}$ Taylor coefficients for approximating the local (un)stable manifolds. The length of the stable and unstable eigenvectors was set to $\epsilon_{u,e_{k}}=0.0565$ and $\epsilon_{s,e_{k}}=0.0635$, respectively. The truncation parameters and the scalings of the eigenvectors were obtained via the procedure in Remark 6.1. The scalings of the eigenvectors were relatively small, since the procedure in Remark 6.1 was designed to use as little Taylor coefficients as possible to ensure that validation is feasible. If we would allow for larger truncation parameters, the scalings of the eigenvectors (and hence the “size” of the charts on the local (un)stable manifolds) could be increased substantially. However, since it is computationally cheaper to increase the integration time in comparison to increasing the truncation parameters for the (un)stable manifolds, we have chosen to keep the truncation parameters $K^{u}$ and (especially) $K^{s}$ small.

• •

Chebyshev approximations: we used $m=3$ subdomains and $N=\begin{bmatrix}50&47&50\end{bmatrix}$ Chebyshev modes. The integration time was set to $L=15$. The Chebyshev coefficients are shown in Figure 6.1a. This figure shows that the decay rates of the Chebyshev coefficients were approximately the same on each subdomain (hence the domain decomposition was successful).

• •

Validation parameters: we used $\nu_{i}=1.1967$ on each subdomain. This value was obtained from the procedure in Remark 6.2.

The dimension of the Galerkin projection was $\dim\left(\mathcal{X}^{NK}\right)=5382$. With the above choices for the computational parameters, we successfully validated a connecting orbit at $\kappa=-1$ and proved that the radii-polynomials were negative for $r\in\left[8.6070\cdot 10^{-11},r^{\ast}\right]$. We remark that it is possible to validate the connecting orbit with a smaller number of Chebyshev coefficients as well. The reason why we used more Chebyshev coefficients than strictly necessary was to get the bounds $Z^{1}_{\Pi^{ij}_{a}}$ as small as possible with an eye towards future work, namely in order to make continuation with large step-sizes feasible.

Discrete continuation

Next, we continued the connecting orbit at $\kappa=-1$ by performing pseudo-arc length continuation. At each continuation step, we tried to validate the orbit with the same computational parameters. We succeeded in validating a family of connecting orbits in this way for a finite number of wave speeds $\kappa\in\left[-1,-0.7861\right]$, see Figure 6.2a.

The reader is referred to the code for the exact parameter values $\kappa$ at which the connecting orbits were validated. Here we only give rounded values of $\kappa$ using four decimal places.

We were not able to validate the connecting orbit at the next continuation step $\kappa\approx-0.7767$ with the same computational parameters. The reason for this was that the bound for $\Pi_{Q}^{4}\left(DF\left(\hat{x}\right)-\widehat{DF}\right)$ became too large, which was related to the fact that the decay rates of $\hat{q}$ decreased as $\kappa$ increased. In addition, the real part of the stable eigenvalues decreased as well when $\kappa$ increased (see Figure 6.3), which contributed to the deterioration of the bounds for the stable manifold. On the other hand, the bounds for $\Pi_{a}^{ij}\left(DF\left(\hat{x}\right)-\widehat{DF}\right)$ did not deteriorate at all during the continuation process. The main reason for this is that the shape and time parameterization of the orbit remained roughly “the same” throughout the continuation procedure. This caused the decay rates of the Chebyshev coefficients to remain roughly the same as well, as shown in Figure 6.1.

To validate connecting orbits for $\kappa\geq-0.7767$, we recomputed the parameterization of the local stable manifold as explained in Remark 6.1. This resulted in a parameterization with $K^{s}=\begin{bmatrix}12&12&12\end{bmatrix}$ Taylor coefficients. Furthermore, the length of the stable eigenvectors was set to $\epsilon_{s,e_{k}}=0.02275$. The resulting parameterized part of the stable manifold was significantly smaller due to the new scaling of the eigenvectors. To ensure that the endpoint of the connecting orbit was contained in the smaller chart, we integrated the connecting orbit forward in time (while keeping the “initial” starting point on the unstable manifold fixed) and increased the integration time to $L=20$. We then used the procedures from the previous section again to refine the Chebyshev approximations. In particular, we used $m=3$ subdomains, $N=\begin{bmatrix}55&52&62\end{bmatrix}$ Chebyshev modes and used $\nu_{i}=1.1627$ on each subdomain. The dimension of the Galerkin projection was $\dim\left(\mathcal{X}^{NK}\right)=10258$. With these parameter values, we were able to successfully validate a finite number of connecting orbits for $\kappa\in\left[-0.7767,-0.7075\right]$, see Figure 6.2b.

Figure 6.2b shows that there is a bifurcation at $\kappa\approx-0.7071$. To understand what caused this bifurcation, we consider the (approximate) stable eigenvalues $\hat{\lambda}^{s}$ at $\kappa\approx-0.7075$:

[TABLE]

Note that $\hat{\lambda}^{s}_{1}+\hat{\lambda}^{s}_{2}\approx\hat{\lambda}^{s}_{3}$. This provides numerical evidence for the presence of an eigenvalue resonance at $\kappa\approx-0.7071$ and explains the observed bifurcation. In this relatively simple case, one can prove with pen and paper that there is in fact a resonance at $\kappa=-\frac{1}{2}\sqrt{2}$. Hence, in order to validate connecting orbits near $\kappa=-\frac{1}{2}\sqrt{2}$ (and in particular at the resonance point itself), one needs to modify the mapping $F_{Q}$ by conjugating to a nonlinear normal form (instead of just the linear one) as explained in [ManifoldTaylor]. Here we do not pursue this issue any further and leave it as a future research project.

To continue the connecting orbit past the resonance, we set $\kappa=-0.7$ and then continued further from this point. We succeeded in validating connecting orbits for $\kappa\in\left[-0.7,-0.5938\right]$ without changing the computational parameters, see Figure 6.2b. As $\kappa$ increased, the “size” of the chart on the local stable manifold kept decreasing. As a consequence, for $\kappa>-0.5938$, the endpoint of the connecting orbit was too far away from the stable equilibrium in the sense that the bounds related to the equation $u(1)=Q\left(\phi\right)$ were too large (also see Figure 6.4). Although we did not continue any further, we remark that validation for $\kappa>-0.5938$ is feasible by increasing the integration time and recomputing the Chebyshev approximations as before.

An interesting future research project would be to develop algorithms for automatically detecting when the integration time and/or manifolds need to be modified. We believe that the heuristics in Remarks 6.1 and 6.2 would be a good starting point for developing such algorithms.

6.2 Traveling fronts in a fourth order parabolic PDE

We have proven the existence of connecting orbits in (1.5) from $(-1,0,0,0)$ to $(a,0,0,0)$ for $a=-0.1$, $\kappa=-2$ and various values of $\gamma$. Recall that these orbits correspond to traveling fronts of (1.4) with wave speed $\kappa$. The values for $a$ and the wave speed $\kappa$ were chosen through experimentation. We started with a connecting orbit at $\gamma\approx 0.4557$ (see the code for exact parameter value) and then used the procedure as explained at the beginning of this section to select the following computational parameters:

• •

Parameterization mappings: we used $K^{u}=\begin{bmatrix}15&15\end{bmatrix}$ and $K^{s}=\begin{bmatrix}12&12&12\end{bmatrix}$ Taylor coefficients for approximating the local (un)stable manifolds. The length of the stable and unstable eigenvectors was set to $\epsilon_{u,e_{k}}=5.4476\cdot 10^{-2}$ and $\epsilon_{s,e_{k}}=5.3337\cdot 10^{-3}$, respectively.

• •

Chebyshev approximations: we used $m=2$ subdomains and $N=\begin{bmatrix}62&61\end{bmatrix}$ Chebyshev modes. The integration time was set to $L=30$. The decay rates of the Chebyshev coefficients were approximately the same on each subdomain (hence the domain decomposition was successful).

• •

Validation parameters: we used $\nu_{i}=1.1491$ on each subdomain.

The dimension of the Galerkin projection was $\dim\left(\mathcal{X}^{NK}\right)=10314$. With the above choices for the computational parameters, we were able to validate the connecting orbit at $\gamma\approx 0.4557$ and proved that the radii-polynomials were negative for $r\in\left[4.8332\cdot 10^{-10},2.8332\cdot 10^{-6}\right]$.

Next, we performed (non-rigorous) pseudo-arclength continuation and tried to validate the orbits at each discrete continuation step by using the same computational parameters. If the validation failed at a particular continuation step, we determined the cause (as in the previous section) and resolved the issue by modifying the computational parameters. In addition, we also checked in the case of failure whether the dimension of the Galerkin-projection could be significantly reduced by decreasing the integration time and the number of Chebyshev or Taylor coefficients (with special emphasis on reducing the number of Taylor coefficients associated to the stable manifold). The results are summarized in Table 6.1. The corresponding bifurcation curves are shown in Figure 6.5 and the associated traveling wave profiles are depicted in Figure 6.6.

We remark that validation of connecting orbits for $\gamma>10.50$ is feasible; the reason for the “obstruction” $\hat{r}>r^{\ast}$ was that the decay rates of the Chebyshev coefficients decreased as $\gamma$ increased, which eventually resulted in a bound for the residual (i.e. the bound $Y^{ij}_{a}$) that was too large. This issue can be easily resolved by using domain decomposition or increasing the number of Chebyshev coefficients (or by just increasing $r^{\ast}$). Similarly, validation of connecting orbits for $\gamma<0.4557$ (but sufficiently far away from [math]) is feasible as well; the bottleneck for small $\gamma$ is the validation of the local stable manifold. Indeed, as $\gamma$ decreases, the real parts of the stable eigenvalues decrease (see Figure 6.7) and the number of needed Taylor coefficients increases. An interesting future research project would be to determine how close one can get to the “singular” case $\gamma=0$ with the current method.