Chebyshev-Taylor parameterization of stable/unstable manifolds for periodic orbits: implementation and applications

TL;DR

This paper introduces a high-order Chebyshev-Taylor series method combined with the parameterization approach to accurately compute stable and unstable manifolds of periodic orbits in dynamical systems, with applications to celestial mechanics and chaos theory.

Contribution

It develops a recursive linear differential equation framework and spectral boundary value problem solutions for invariant manifolds, enhancing accuracy and robustness over previous methods.

Findings

Effective for Lorenz system analysis

Successful in celestial three-body problems

Able to compute cycle-to-cycle connecting orbits

Abstract

This paper develops a computational method for studying stable/unstable manifolds attached to periodic orbits of differential equations. The method uses high order Chebyshev-Taylor series approximations in conjunction with the parameterization method -- a general functional analytic framework for invariant manifolds. The parameterization method can follow folds in the embedding, recovers the dynamics on the manifold through a simple conjugacy, and admits a natural notion of a-posteriori error analysis. The key to the approach is the derivation of a recursive system of linear differential equations describing the Taylor coefficients of the invariant manifold. We find periodic solutions of these equations by solving a coupled collection of boundary value problems with Chebyshev spectral methods. We discuss the performance of the method for the Lorenz system, and for circular restricted three and four body problems. We also illustrate the use of the method as a tool for computing cycle-to-cycle connecting orbits.

Tables

Table 1. Table 1: Error for the conjugacy test for different Taylor dimension when m = 10 𝑚 10 m=10 , D = 50 𝐷 50 D=50 , K = 250 𝐾 250 K=250 and k 0 = 10 subscript 𝑘 0 10 k_{0}=10 (left). Norm of the last Taylor dimension for different choices of K 𝐾 K when m = 10 𝑚 10 m=10 , D = 50 𝐷 50 D=50 , N = 100 𝑁 100 N=100 and k 0 = 10 subscript 𝑘 0 10 k_{0}=10 (right).

$N$	$t_0={10}^{-5}$	$t_0=1$
$20$	$9.1859\times {10}^{-4}$	$4.9280\times {10}^{-1}$
$40$	$4.2205\times {10}^{-5}$	$1.2568\times {10}^{-2}$
$60$	$1.2725\times {10}^{-6}$	$2.6406\times {10}^{-4}$
$80$	$3.1854\times {10}^{-8}$	$8.0463\times {10}^{-6}$
$100$	$7.2073\times {10}^{-10}$	$3.404\times {10}^{-8}$

Equations

$$\dot{x}=g(x),$$

$$\lambda \in \mathbb{R},\text{mbox}with\lambda <0,$$

$$\frac{d}{dt}v(t)=Dg(\gamma (t))v(t)-\lambda v(t),$$

$$\frac{\partial }{\partial t}P(t,\sigma )+\lambda \sigma \frac{\partial }{\partial \sigma }P(t,\sigma )=g(P(t,\sigma )),$$

$$P(t,0)=\gamma (t),$$

$$\frac{\partial }{\partial \sigma }P(t,0)=v(t).$$

$$x(t)=P(t,e^{\lambda t}{\sigma }_0).$$

$$\frac{d}{dt}x(t)$$

$$=\frac{\partial }{\partial t}P(t,e^{\lambda t}{\sigma }_0)+\lambda {\sigma }_0{e}^{\lambda t}\frac{\partial }{\partial \sigma }P(t,e^{\lambda t}{\sigma }_0)$$

$$=g(P(t,e^{\lambda t}{\sigma }_0))$$

$$=g(x(t)),$$

$$\underset{t\to \infty }{\lim }P(t_0,e^{\lambda t}{\sigma }_0)$$

$$=P(t_0,0)$$

$$=\gamma (t_0),$$

$$P(s+t,e^{\lambda t}\sigma )=\Phi (P(s,\sigma ),t)$$

$$P(t,\sigma )=\sum _{\alpha =0}^{\infty }{A}_{\alpha }(t){\sigma }^{\alpha },$$

$$g(x,y,z)=\left(\begin{array}[]{c}\sigma(y-x)\\ x\rho-xz-y\\ xy-\beta z\end{array}\right),$$

$$P(t,s)=\sum_{\alpha=0}^{\infty}A_{\alpha}(t)s^{\alpha}=\sum_{\alpha=0}^{\infty}\left(\begin{array}[]{c}a^{(1,\alpha)}(t)\\ a^{(2,\alpha)}(t)\\ a^{(3,\alpha)}(t)\end{array}\right)s^{\alpha}$$

$$A_0(t)=\gamma (t),\text{mbox}and{A}_1(t)=v(t).$$

$$\frac{\partial }{\partial t}P(t,s)+\lambda s\frac{\partial }{\partial s}P(t,s)=\sum _{\alpha =0}^{\infty }\left(\frac{d}{dt}{A}_{\alpha }(t)+\alpha \lambda A_{\alpha }(t)\right)s^{\alpha },$$

$$g(P(t,s))$$

$$\frac{d}{dt}{A}_{\alpha }(t)+\alpha \lambda A_{\alpha }(t)$$

$$\displaystyle=\left(\begin{array}[]{c}\sigma\left(a^{(2,\alpha)}-a^{(1,\alpha)}\right)\\ \rho a^{(1,\alpha)}-a^{(2,\alpha)}-a^{(1,\alpha)}a^{(3,0)}-a^{(3,\alpha)}a^{(1,0)}\\ -\beta a^{(3,\alpha)}+a^{(1,\alpha)}a^{(2,0)}+a^{(2,\alpha)}a^{(1,0)}\end{array}\right)+\left(\begin{array}[]{c}0\\ -\sum_{\beta=1}^{\alpha-1}a^{(1,\alpha-\beta)}a^{(3,\beta)}\\ \sum_{\beta=1}^{\alpha-1}a^{(1,\alpha-\beta)}a^{(2,\beta)}\end{array}\right).$$

$$\displaystyle\left(\begin{array}[]{c}\sigma\left(a^{(2,\alpha)}-a^{(1,\alpha)}\right)\\ \rho a^{(1,\alpha)}-a^{(2,\alpha)}-a^{(1,\alpha)}a^{(3,0)}-a^{(3,\alpha)}a^{(1,0)}\\ -\beta a^{(3,\alpha)}+a^{(1,\alpha)}a^{(2,0)}+a^{(2,\alpha)}a^{(1,0)}\end{array}\right)$$

$$=Dg(\gamma (t))A_{\alpha }(t),$$

$$g_{\alpha}(A_{1}(t),\ldots,A_{\alpha-1}(t)):=\left(\begin{array}[]{c}0\\ -\sum_{\beta=1}^{\alpha-1}a^{(1,\alpha-\beta)}(t)a^{(3,\beta)}(t)\\ \sum_{\beta=1}^{\alpha-1}a^{(1,\alpha-\beta)}(t)a^{(2,\beta)}(t)\end{array}\right).$$

$$\frac{d}{dt}{A}_{\alpha }(t)-\left(Dg(\gamma (t))-\lambda \alpha \text{mbox}Id\right)A_{\alpha }(t)=g_{\alpha }(P),$$

$$P^N(t,s)=\sum _{\alpha =0}^N{A}_{\alpha }(t)s^{\alpha }.$$

$${\pi }_d(R(\mathbf{x}))=\mathbf{x},$$

$$DR(\mathbf{x})f(\mathbf{x})=F(R(\mathbf{x})),\text{mbox}forall\mathbf{x}\in U.$$

$${\pi }_{d}F(R(\mathbf{x}))=f(\mathbf{x}),$$

$$x^{\prime }$$

$$y^{\prime }$$

$$z=\frac{1}{|x|}=\frac{1}{\sqrt{x^2}},$$

Full text

Key words. periodic orbit, (un)stable manifold, parameterization method, boundary value problems, automatic differentiation, Chebyshev polynomials

1 Introduction

Periodic solutions of differential equations are the basic building blocks of recurrence in nonlinear dynamics. Moreover, hyperbolic periodic orbits and their heteroclinic/homoclinic connections are natural generators of chaotic motions. Since a heteroclinic/homoclinic orbit will approach a periodic solution along its local stable/unstable manifolds, computational methods for studying these manifolds are of great interest. A schematic description of the stable manifold of a periodic orbit, beside an actual stable manifold in the Lorenz system are illustrated in Figure 1. See any of the works Canalias and Masdemont, [2006]; Font et al., [2009]; Jorba, [1999]; Jorba and Masdemont, [1999]; Jorba and Villanueva, [1998]; Krauskopf et al., [2005]; Masdemont, [2011]; Osinga, [2003]; Simó, [1998, 1988, 1989] for more discussion, but we caution that any such list will hardly scratch the surface of the relevant literature.

The stable/unstable normal bundles of a periodic solution approximate the stable/unstable manifolds to first order, and these bundles are obtained by studying the equations of first variation – or equivalently – by solving certain periodic eigenvalue problems. Higher order jets of the invariant manifold could be studied via higher order equations of variation, however the complexity of these equations grows exponentially with the order of the jet. More efficient methods for studying the jets are obtained by reformulating the invariant manifold as the solution of an operator equation, and studying the operator equation via numerical methods.

The parameterization method is a general functional analytic framework for studying invariant manifolds Cabré et al., 2003a ; Cabré et al., 2003b ; Cabré et al., [2005]; Canadell and Haro, [2014]; Figueras and Haro, [2012, 2013]; Haro and de la Llave, 2006a ; Haro and de la Llave, 2006b ; Haro and de la Llave, [2007] whose goal is to find a chart/covering map conjugating the dynamics on the invariant manifold to a simple and well understood model (correctly choosing this model is part of the method). By viewing the conjugacy as an operator equation for the unknown chart, the problem is susceptible to numerical methods. We will see below that the operator equation of interest in the present work is actually a PDE with prescribed periodic data.

This operator equation is referred to as the invariance equation. The unknown parameterization is not required to be the graph of a function, and hence is able to follow folds in the embedding. See again Figure 1. Since the invariance equation is based on a dynamical conjugacy, the parameterization method recovers the dynamics on the manifold in addition to the embedding. By now there is a small industry devoted to the parameterization method, and a review of the literature would take us far afield. We refer the interested reader to the recent book on the subject Haro et al., [2016], where many examples and much more complete discussion of the literature is found.

The present work is not the first numerical treatment of parameterized stable/unstable manifolds attached to periodic orbits, and indeed we build on the earlier studies of Cabré et al., [2005]; Castelli et al., [2015]; Guillamon and Huguet, [2009]; Huguet and de la Llave, [2013]. As in these earlier studies, we make a formal Taylor series arguments which analytically reduces the invariance equation before making any numerical computations. Since the coefficients of the Taylor series are themselves periodic functions, the formal power matching scheme leads to so called homological equations describing the unknown coefficients. In the present case of a periodic orbit, the homological equations are linear ordinary differential equations with periodic coefficients and periodic forcing. We solve these linear homological equations recursively using numerical spectral methods. By computing the formal series solution to high order we obtain an approximate solution which describes the stable/unstable manifold far from the periodic orbit.

In the earlier works just cited the differential equations – describing the periodic orbit, describing the normal bundles, and the homological equations describing the higher order jets – are all solved using Fourier spectral methods. Fourier methods are both efficient and accurate when applied to periodic solutions of moderate length. This efficiency is due in part to the fact that differentiation is a diagonal operation in the transform domain, and in part to the fact that the FFT speeds up evaluation of nonlinearities. However, the decay rate of the Fourier coefficients gets increasingly slow as the period/harmonic complexity of the orbit grow. In practice this means that it is necessary to compute more and more Fourier coefficients, and for long enough orbits the Fourier approximation becomes impractical.

There is much recent interest in numerical methods based on Chebyshev spectral approximation of solutions to boundary value problems. We refer the interested reader to Driscoll et al., [2008]; Platte and Trefethen, [2010]; Trefethen, [2007, 2013] and the references discussed there in. The present work builds on techniques developed by a number of authors which use Chebyshev spectral methods to compute long periodic solutions of differential equations Gameiro et al., [2016]; Lessard and Reinhardt, [2014]; van den Berg and Sheombarsing, [2016]. The merit of this approach is that the Chebyshev spectral methods posess many of the advantages of Fourier series – for example differentiation is a tri-diagonal operation in the transform domain and the fast cosine transform is available for evaluating nonlinearities – but Chebyshev series apply to non-periodic boundary value problems. Treating a periodic solution as a series of coupled boundary value problems – on smaller domains – facilitates control of the decay rates of the coefficients.

Motivated by these developments, the present work applies Chebyshev methods not only the periodic orbit – but also to the computation of the normal bundles and the homological equations for the higher order jets. The result is a computational method for finding Chebyshev-Taylor expansions of the local stable/unstable manifolds attached to periodic orbits. Our method applies to more complicated orbits and their attached invariant manifolds than could be studied using only the Fourier-Taylor approach.

Remark 1.1** (Connecting orbits and extensions of local stable/unstable manifolds).**

Of course computing local stable/unstable manifolds attached to periodic orbits is only a means to an end. In applications we are often interested in either using the local manifolds to compute connecting orbits, or to grow larger portions of the invariant manifold in order to study the global dynamics. While in the present work we do consider a number of example computations for connecting orbits, we do not make any serious attempt to numerically grow larger local manifolds. This is because the literature on computational methods for growing invariant manifolds is extensive and well developed. The interested reader will want to consult the review paper Krauskopf et al., [2005] for a thorough overview of the literature, and will find other powerful methods and fuller discussion in England et al., [2005]; Osinga, [2000, 2003]. We only note that the methods developed in the present work could be combined with existing continuation methods for even better results. This is especially true for methods which exploit the curvature or other differential geometric properties of the manifold. **

Remark 1.2** (Automatic differentiation and polynomial nonlinearities).**

Multiplication of Taylor and Chebyshev series is straight forward thanks to the Cauchy product in the former case and the discrete cosine convolution operation in the later. Then formal series manipulations for polynomial nonlinearities are especially transparent in Chebyshev-Taylor bases. In the present work we are interested in applications coming from celestial mechanics which involve non-polynomial vector fields. In order to simplify matters we exploit methods of automatic differentiation and transform to the polynomial setting, albeit in a higher dimensional phase space. This is discussed in detail in Section 2.3.

The use of automatic differentiation is a convenience rather than a necessity, as the FFT could be used to evaluate general nonlinearities. In fact, even after automatic differentiation we use the fast cosine transform to evaluate higher order polynomial nonlinearities in the present work. Nevertheless, the use of automatic differentiation in the present work simplifies the implementation details of our algorithms – as all our computations are reduced to Newton’s method for large polynomial systems. Automatic differentiation also simplifies a-posteriori error analysis for the method, which when followed to its logical conclusion provides mathematically rigorous validated error bounds for the parameterizations. **

2 Review of the parameterization method

As already mentioned in the introduction, the parameterization method is much more general than what we actually use in the present work. We refer the reader again to the book Haro et al., [2016]. In the following section we review some basic notions in the very simple setting of an orientable local manifold associated with one stable/unstable Floquet exponent. Generalities such as multiple stable/unstable exponents, complex conjugate exponents, and non-orientable bundles are discussed in detail in Castelli et al., [2015]. The methods of the present work apply in these more general setting with only obvious modifications. We focus on the one dimensional case to simplify the exposition.

2.1 Parameterization of stable/unstable manifolds attached to periodic orbits

Let $\Omega\subset\mathbb{R}^{M}$ be an open set and $g\colon\Omega\to\mathbb{R}^{M}$ be a real analytic vector field. Suppose that $\gamma\colon\mathbb{R}\to\mathbb{R}^{M}$ is a $T$-periodic solution of the first order ordinary differential equation

[TABLE]

that is we assume that $\gamma^{\prime}(t)=g(\gamma(t))$ with $\gamma(t+T)=\gamma(t)$ for all $t\in\mathbb{R}$. Suppose also that $\gamma$ has one stable Floquet exponent

[TABLE]

so that (by the stable manifold theorem) there exists a two dimensional manifold of solutions which converge exponentially fast to the periodic orbit $\gamma$. Let $v\colon\mathbb{R}\to\mathbb{R}^{M}$ denote the stable normal bundle of $\gamma(t)$, associated with the exponent $\lambda$. We assume that $v$ is an orientable bundle, so that $v(t)$ is $T$ periodic as well. We note that $(v,\lambda)$ solve the eigenvalue problem

[TABLE]

subject to some normalization, perhaps $\|v(t)\|=1$ for all $t\in\mathbb{R}$ (though in numerical applications we will choose other normalizations).

The goal of the parameterization method is to find a smooth function $P\colon[0,T]\times[-1,1]\to\mathbb{R}^{M}$ solving the invariance equation

[TABLE]

subject to the first order constraints

[TABLE]

and

[TABLE]

Then geometric content of Equation (1) is illustrated in Figure 2, but one easily checks that the image of $P$ is a stable manifold.

Indeed, let $P$ be a smooth solution of Equation (1) subject to the first order constraints. Choose any $\sigma_{0}\in(-1,1)$ and define the curve $x\colon[0,\infty)\to\mathbb{R}^{M}$ by

[TABLE]

Then

[TABLE]

as $e^{\lambda t}\sigma_{0}\in(-1,1)$ for all $t\geq 0$. Then $x(t)$ is a solution curve for the differential equation having $x(0)=P(0,\sigma_{0})$. Moreover for any $t_{0}\in[0,L]$, since $P$ is continuous we have that

[TABLE]

that is, a point on the image of $P$ accumulates at the periodic orbit $\gamma$ with asymptotic phase $\gamma(t_{0})$. In particular, the image of $P$ is a local stable manifold for $\gamma$.

One can actually prove more. For example if $P$ solves Equation (1) subject to the first order constraints, then $P$ actually satisfies the flow conjugacy

[TABLE]

for all $t\geq 0$. Here $\Phi$ is the flow generated by the vector field $g$. The meaning of this flow conjugacy is illustrated in Figure 3. The proof of the flow conjugacy is given for example in Castelli et al., [2015]. In the same reference it is shown that solutions of Equation (1) are unique up to the choice of the eigenfunction $v(t)$. (Any constant multiple of $v$ is a parameterization of the stable normal bundle, but up to this choice of scaling the solution is unique). If the solution $P$ exists, it is as regular as $g$. In this case $P$ is real analytic if $g$ is Cabré et al., 2003b ; Cabré et al., [2005]. Moreover, in the case of one stable exponent, there exists a choice of scaling small enough that the solution $P$ exists, and is analytic. See Cabré et al., 2003a ; Cabré et al., [2005]; Huguet and de la Llave, [2013].

2.2 Formal series and the reduction to homological equations

Since $P$ exists and is analytic it makes sense to seek a power series solution

[TABLE]

where $\sigma\in(-1,1)$ and the functions $A_{\alpha}\colon\mathbb{R}\to\mathbb{R}^{M}$ are analytic and $T$-periodic. Plugging the power series into Equation (1), expanding the nonlinearities, and matching like powers leads to equations for the unknown Taylor coefficient functions.

Example: the Lorenz system in the specific case of the Lorenz field we have that $M=3$, and that the vector field $g$ is given by

[TABLE]

where $\beta,\rho,\sigma$ are positive constants. Suppose that $\gamma(t)=(x(t),y(t),z(t))$ is an analytic period- $T$ orbit with stable (or unstable) exponent $\lambda$ and that $v(t)=(v_{1}(t),v_{2}(t),v_{3}(t))$ is an analytic $T$-periodic parameterization of the stable normal bundle. We look for

[TABLE]

Then

[TABLE]

Since

[TABLE]

and

[TABLE]

equating and matching like powers of $s$ leads to

[TABLE]

Noting that

[TABLE]

define the functions $g_{\alpha}$ by

[TABLE]

We write $g_{\alpha}(P)$ for short. Now we seek $A_{\alpha}(t)$ the $T$-periodic solution of the equation

[TABLE]

which we refer to as the homological equation for $P$. Note that this is a linear inhomogeneous first order ordinary differential equation with periodic coefficients, and that the right hand side is independent of $A_{\alpha}$. Indeed, $g_{\alpha}$ depends recursively on lower order terms. The Floquet theory guarantees that our homological equation has a unique periodic solution for each $\alpha\geq 2$. Then we recursively solve the equations to order $N$ and have the approximate solution

[TABLE]

Remark 2.1** (A-posteriori error analysis).**

Truncation error analysis is treated carefully in Castelli et al., [2017]. Note that the analysis in that reference is independent of the basis used to represent $A_{\alpha}(t)$: only the implementation exploits that these functions are given as Fourier series. Then the methods of the work just cited apply directly to the expansions used in the present work. The key to the analysis in Castelli et al., [2017] is that the approximation $P^{N}$ have small defect. In the present work we only check the defect numerically, and postpone to an upcoming work more careful analysis of the errors for our Chebyshev-Taylor approximations. **

The calculations above generalizes to any polynomial vector field in the obvious way, and we have that the $A_{\alpha}$ satisfy homological equations of exactly the form given in Equation (6). Only the term $Dg(\gamma(t))$ and the form of the recursive functions $g_{\alpha}(P)$ depend explicitly on the form of the vector field $g$. In the other examples in the present work we simply write down the correct homological equations and leave the derivations as an exercise for the interested reader.

2.3 A dynamical perspective on automatic differentiation

The term automatic differentiation refers to a whole suite of methods for managing the complexity of problems involving high order derivatives. These methods exploit the fact that differentiation and multiplication are related through the chain rule. For a general introduction and an overview of the literature we refer to Rall and Corliss, [1996]; Bücker and Corliss, [2006], though in the present work our use of the term is much closer to that of Chapter $2$ of Haro et al., [2016] and also Jorba and Zou, [2005]. See also Chapter $4.7$ of Knuth, [1998]. The discussion of the literature in this last reference is especially illuminating, though the reference does not use the “automatic differentiation” terminology.

While automatic differentiation can be viewed as a technique for computing transcendental functions of polynomials (where the result should again be expressed as a polynomial), automatic differentiation can also bee seen as a technique transforming non-polynomial into polynomial problems by appending additional polynomial differential equations. The solutions of the appended differential equations are required to give the non-polynomial terms of the original problem. When viewed this way it is possible to use automatic differentiation to solve problems involving Fourier and Chebyshev rather than just Taylor bases. Such extensions are discussed for example in Lessard et al., [2016], where automatic differentiation is simply viewed as a change of variables which transforms the given problem to a polynomial problem.

In general a change of variables will disturb the orbit structure of a differential equation, unless the change of variables induces a dynamical (semi)-conjugacy. This issue is considered at length in Kepley and Mireles James, [2017], where a dynamical systems interpretation of automatic differentiation is introduced. We adopt this interpretation here as well. So, suppose that $U\subset\mathbb{R}^{d}$ is an open set and that $f\colon U\to\mathbb{R}^{d}$ is a real analytic vector field with non-polynomial nonlinearity. From the perspective of the present work, the goal of automatic differentiation is to find a mapping $R\colon U\to\mathbb{R}^{D}$ with $D>d$ and a polynomial vector field $F\colon\mathbb{R}^{D}\to\mathbb{R}^{D}$ having that

[TABLE]

where $\pi_{d}\colon\mathbb{R}^{D}\to\mathbb{R}^{d}$ is the projection into the lower dimensional space, and that

[TABLE]

This leads to the following observations.

• $\bullet$

Since $\pi_{d}R=\mbox{Id}_{d}$, it follows that

[TABLE]

so that the original field is recovered by projection.

• $\bullet$

Equation (7) is an infinitesimal conjugacy, and implies that $R$ maps orbits of $f$ to orbits of $F$. Moreover it says that graph of $R$ is invariant under the flow generated by $F$.

• $\bullet$

Then when we solve the differential equations given by $F$ with initial conditions on the graph of $R$, we can recover the orbits of $f$ simply by projecting.

The procedure for choosing $R$ and $F$ is best illustrated through examples, but it is worth noting that $R$ will have the same domain as $f$. So even though $F$ is polynomial and hence entire, the composition of $F\circ R$ will have the same singularities as $f$.

Example: Consider the Kepler problem

[TABLE]

We introduce the new variable

[TABLE]

and note that away from the origin one has

[TABLE]

Then taking

[TABLE]

one easily checks the conjugacy Equation (7), and also that the first two components of $F\circ R$ recover $f$. Solution curves of the vector field $F$ with initial conditions on the graph of $R$ recover solutions of the Kepler problem as long as they remain on the graph of $R$. That is, if

[TABLE]

then $x(t),y(t)$ is a solution of the Kepler problems as long as $x(t),y(t),z(t)$ is on the graph of $R$. Moreover, the trajectory leaves the graph of $R$ if and only if there is a collision in the the Kepler problem. This example is typical of automatic differentiation for problems in celestial mechanics. A thorough discussion – from the dynamical systems point of view – of automatic differentiation for the circular restricted four body problem is found in Kepley and Mireles James, [2017].

3 Chebyshev Expansion for Periodic Solutions of $\dot{u}(t)=h(u(t),t)$

The periodic orbit $\gamma$, its stable/unstable normal bundle $v(t)$, and the higher order Taylor coefficients $A_{\alpha}(t)$ for $2\leq\alpha\leq N$ are all periodic solutions of non-autonomous differential equations of the form

[TABLE]

with period $\tau$, and where $h:\mathbb{R}^{M}\times\mathbb{R}\to\mathbb{R}^{M}$. In each case we look for a solutions $v:[0,\tau]\to\mathbb{R}^{M}$ satisfying the boundary value problem

[TABLE]

We develop a Chebyshev scheme to solve this class of problems.

To begin, break the solution into sub-pieces using the following mesh. Let $0=t_{0}<t_{1}<t_{2}<\ldots<t_{D}=\tau$, and for $i=1,\ldots,D$ define $v^{i}(t)$ on $[t_{i-1},t_{i}]$ such that

[TABLE]

Thus, each $v^{i}(t)$ is a solution of $\dot{u}(t)=h(u(t),t)$. Moreover, at each of the points $t_{i}$, $1\leq i\leq D-1$, two different pieces are defined and they must agree. Then we impose the boundary conditions

[TABLE]

for $i\neq 1$, and to impose periodicity

[TABLE]

We want to expand each piece $v^{i}(t)$ using Chebyshev polynomials. To do so, we first rescale the problem to the interval $[-1,1]$. First, note that since there are no time dependence in $g$ (the original system for which we compute the periodic orbit) we can translate the time domain to some interval $[-L_{i},L_{i}]$ and then rescale time $t\mapsto\frac{t}{L_{i}}$, so that the solution satisfies

[TABLE]

for all $i=1,\ldots,D$.

For the time varying case, i.e. the case of the bundles and the homological equations, the solution will depend on the lower power $\alpha$ of the coefficient in question. Since the bundles and Taylor coefficients of the parameterization are all periodic with the same period, we choose a fixed mesh for all the problems. That is, the number of subdomains $D$ will be the same at every step. Moreover each $L_{i}$ is a fixed proportion of the global period. So $L_{i}=p_{i}L$ for $p_{i}\in(0,1]$ some given constant. In applications unless otherwise specified we use a uniform mesh, so that $p_{i}=\frac{1}{D}$ for all $i$.

We introduce a Chebyshev expansion for each sub-piece $v^{i}$ that are now defined on $[-1,1]$. Let $v^{(i,j)}:[-1,1]\to\mathbb{R}$ denote the $j-$th component of $v^{i}$ for all $i=1,\ldots,M$. For any $i=1,\ldots,D$ and $j=1,\ldots,M$, we set

[TABLE]

where $T_{k}$ is the $k-$th Chebyshev polynomial, which are defined as follows.

Definition 3.1**.**

The Chebyshev polynomials $T_{k}:[-1,1]\rightarrow\mathbb{R}$, $k=0,1,2,\dots$ are defined by $T_{0}(t)=1$ and $T_{1}(t)=t$ and the recurrence relation

[TABLE]

It is well known that these polynomials satisfy $T_{k}(\cos\theta)=\cos(k\theta)$, a fact which can be used to prove further results relating the Chebyshev series to results from Fourier analysis.

Definition 3.2**.**

We denote the set of unknown Chebyshev coefficients of the full periodic orbit $v$ by

[TABLE]

So that for every $i=1,2\ldots,M$, $A^{(i)}$ represent the set of coefficients of the $i-$th piece of $v$ defined on the interval $[t_{i-1},t_{i}]$. Moreover, for a fixed $k$, $A_{k}^{(i)}\in\mathbb{R}^{M}$ since each has an image included in $\mathbb{R}^{M}$, which are all expanded using distinct Chebyshev expansions. Finally, $a_{k}^{(i,j)}$ denotes the $k-$th coefficient of the $j-$th dimension of $v^{i}$.

To rewrite the system as an operator defined on the set $A$ we integrate (10) from $-1$ to $t$ and obtain

[TABLE]

Note that $h:[-1,1]\to\mathbb{R}^{M}$, since it depends on the lower and current term of the expansion of the parameterization which all are defined on $[-1,1]$ after rescaling time. Thus it can also be expanded using Chebyshev polynomials. We set

[TABLE]

and substitute both Chebyshev expansions in (12) to get an equation whose only time dependence is in the Chebyshev polynomials themselves/the integral. We use the recurrence formulas for the integral of the Chebyshev polynomials and rewrite the initial condition so that after simplification the Chebyshev coefficients need to satisfy a set of conditions defined in the space of Chebyshev coefficients. That is, for all $i=1,\ldots,D$, $j=1,\ldots,M$ and $k\geq 0$, we define

[TABLE]

Each $f_{k}^{(i,j)}(L,A)$ is given by

[TABLE]

Definition 3.3**.**

We write $c_{k\pm 1}$ to denote $c_{k\pm 1}:=c_{k+1}-c_{k-1}$.

We omit the derivation of the operators $f_{k}^{(i,j)}$. For further details, and to see why the term $k=0$ does not depend on the vector field but only on the boundary condition, we refer to Lessard and Reinhardt, [2014].

The boundary condition provides the initial component for the new problem, that is for the case $f_{0}^{(i,j)}$. So we use the form given in Equation (13) for every subintervals except for the case $i=1$, for which we use

[TABLE]

The operator involves coefficients of the Chebyshev expansion of $h$, which we need to write in terms of the unknowns $A$. Since the Chebyshev expansion of a sum (or difference) of two functions $f$ and $g$ will be given by the sum (or difference) of the Chebyshev coefficients of $f$ and $g$, and since we assume $h$ to be polynomial, the only other case to handle is a product. We use the following Lemma.

Lemma 3.4**.**

If $f(t)$ and $g(t)$ are expanded with Chebyshev series so that

[TABLE]

then $f(t)g(t)=(a\ast b)_{0}+2\displaystyle\sum_{k=1}^{\infty}(a\ast b)_{k}T_{k}(t)$. Here $\ast$ denotes the discrete convolution product

[TABLE]

Therefore, we can completely rewrite each $f_{k}^{(i,j)}$ without the coefficients $c_{k}$. To compute an approximation of the solution, we consider a finite dimension approximation

[TABLE]

so that $F:\mathbb{R}^{mMD}\to\mathbb{R}^{mMD}$. Note that the boundary conditions introduce a dependence between each subdomain, hence we solve the equations simultaneously. In the case $\alpha\geq 2$ one can directly use Newton’s method to find $\bar{A}$ such that $F(\bar{A})\approx 0$, where $\bar{A}$ is a finite dimensional approximation of the unknowns $A$ with same dimension as $F$.

The first order data: In the case $\alpha=0$ or $\alpha=1$, we have additional unknowns hence we also need to add phase conditions. For the periodic orbit $\gamma$, the period $\tau$ is an unknown. For the eigenvalue equation defining the stable/unstable bundles it is the eigenvalue/Floquet exponent which is unknown. In both cases we must balance the equations.

For the periodic orbit $\gamma(t)$, we replace $h(v(t),t)$ by the given vector field $g(v(t))$. Note that for any $s$, the time translated solution $\gamma_{s}(t)=\gamma(t+s)$ is still a solution of the problem since it is periodic and satisfies $\dot{\gamma_{s}}(t)=g(\gamma_{s}(t))$. In order to isolate a solution we impose a Poincaré condition

[TABLE]

for some $p_{0}\in\mathbb{R}^{M}$ with $\dot{p_{0}}:=g(p_{0})$. This condition translates into a condition on $v^{1}(-1)$ after applying the proper change of variable to Chebyshev coefficients as

[TABLE]

The normal bundle satisfies the eigenvalue problem

[TABLE]

where $Dg(\gamma(t))$ is the periodic matrix given by the derivative of $g$ evaluated at the periodic orbit. It follows that any rescaling $kv(t)$ of the bundle is again a solution associated to the same eigenvalue $\lambda$. To isolate a solution we fix $\|v(0)\|=K$, where we are free to choose $K$. In term of the Chebyshev coefficients

[TABLE]

For simplicity we truncate this condition to

[TABLE]

which still isolates an eigenfunction.

4 Examples

We introduce the following operator to simplify the expansion of the operator $F$ in each example.

Definition 4.1**.**

Let $(j_{1},\ldots,j_{n})$ be a set consisting of hyperscript corresponding to component of the solution. That is $1\leq j_{k}\leq M$ for all $k=1,\ldots,n$, we denote their Cauchy product of convolutions by

[TABLE]

Note that the case $\alpha=0$ simply returns the convolution product.

4.1 The example of the Lorenz system

All of our numerical computations use the classical parameter values $\beta=\frac{8}{3}$, $\rho=27$ and $\sigma=10$. The operator defining the unknowns $A_{\alpha}(t)$ is given by

[TABLE]

with $k\geq 1$. The formula for $k=0$ is omitted since it is already explicitly given in (13). Note that $\mathcal{C}$ involves the lower order terms $A_{\beta}$ for $\beta\leq\alpha$ but we are only solving for $A_{\alpha}$ with the lower order terms fixed.

Only the $\alpha=0$ problem is nonlinear in the unknown Chebyshev coefficients and this requires a good initial guess to obtain a periodic orbit. In the present work we use the data from Viswanath, [2003] as the input for a Newton method. For each $\alpha$ we truncate and consider the finite dimensional problem

[TABLE]

The corresponding truncated operator is such that $F:\mathbb{R}^{3mD+1}\to\mathbb{R}^{3mD+1}$ for $\alpha=0,1$ and $F:\mathbb{R}^{3mD}\to\mathbb{R}^{3mD}$ in the remaining cases.

In Figure 5 we computed the stable manifold for two periodic orbits. For these computations we use $D=10$ Chebyshev domains, $m=50$ Chebyshev coefficients per domain, and $N=50$ Taylor nodes. The time rescaling in the operator (i.e. half the period) is $L\approx 1.1530$ for the shorter orbit and $L\approx 1.9101$ for the longer one. Each color represent a different subdomain of the Chebyshev time decomposition.

In Figure 6 we extended the manifold on the right of figure 5 by integrating backward in time $100$ points evenly distributed on the boundary of the parameterization. The orbits were computed by integrating backward in time for $t=0.5$. Using the conjugacy relation (4), we have that one would need to integrate for $t=2.5212$ to go from $\sigma\approx 10^{-16}$ to the boundary $\sigma=1$. Thus, allowing the utilization of much smaller time lapse to get a good extension of the attractor for this orbit. The discretization of the continued manifold is very coarse, but the example is only included to show that nearly all the “slow” dynamics of the manifold is captured by the parameterization. Once we start integrating the local manifold orbits move away very rapidly.

A-posteriori error analysis: Since computing a local invariant manifold with the parameterization method does not necessarily involve a small parameter, its not always clear what it means to talk about “the order” of the error. This is especially true when polynomials of high degree are used to approximate the manifold in a large neighborhood of the periodic orbit. In this case it is more natural to use the notion of defect to quantify errors, and we exploit the fact that the parameterization satisfies a conjugacy equation.

That is, since the equation (4) must be satisfied, we define

[TABLE]

and $t_{0}$ is some fixed test time. Sampling points in $[0,\tau]$ with $\sigma=\pm 1$ leads to a useful and numerically accessible estimate.

Some heuristics are also helpful. For example we find that choosing $N$ and $K$ so that the norm of the last Taylor component is around machine precision leads to excellent results. A useful norm for making this assessment is

[TABLE]

as this involves only sums of the known coefficients.

We test the conjugacy for $t_{0}=10^{-5}$ and $t_{0}=1$ with $200$ different starting points evenly distributed on the parameterization of the manifold and then took the average of the resulting errors. The results are displayed in table 1. From this table, one can note that the choice of $N$ and $K$ are such that the conjugacy error remains considerably small for longer period of time while the norm of the last Taylor sequence is not too far beyond machine precision.

4.1.1 Short Connecting orbit

Following the convention established in Lessard et al., [2014] we say that there is a short connection from $\gamma_{1}$ to $\gamma_{2}$ if the local parameterization of the unstable manifold of $\gamma_{1}$ intersects the local parameterization of the stable manifold of $\gamma_{2}$. (Here we are assuming that $\gamma_{1,2}$ come equipped with some choice of local manifold parameterizations. In fact it is more correct to say that the connection is short relative to these fixed parameterizations). In this case one can establish the existence of a connection without the use of any numerical integration. The fact this is possible provides another illustration of the fact that our manifold computations are in some sense not local.

In Figure 7 we display the stable manifold of the orbit AB and the unstable manifold of ABB in Lorenz and observe several intersection of the two manifolds. The boundary of the unstable manifold crosses the stable manifold giving a connecting orbit. Let $P:[0,2L_{1}]\times[-1,1]\to\mathbb{R}^{3}$ denote the parameterization of the local stable manifold of the periodic orbit AB, whose period is $2L_{1}$. Similarly, let $Q:[0,2L_{2}]\times[-1,1]\to\mathbb{R}^{3}$ be the parameterization of the unstable manifold of the periodic orbit ABB, whose period is given by $2L_{2}$. Once the unstable manifold is restricted to one of its boundary circles only $3$ unknowns remain. We set $x=(\theta_{s},\sigma_{s},\theta_{u})$ so that the desired intersection is a zero of

[TABLE]

Therefore $S:\mathbb{R}^{3}\to\mathbb{R}^{3}$ and it is possible to apply Newton’s method to obtain an approximation of the solution. Using this approach, we found

[TABLE]

Then the connecting orbit can be computed from those coordinates in the parameter space. Since the value $\sigma_{s}$ of $x$ is already quite small we used the conjugacy relation (4) with $t=2$ and it is already sufficient. For the other half of the connecting orbit, we use (4) again to determine the “time of flight”. In this case, the unstable eigenvalue is $\lambda\approx 0.9947$, so one would need to integrate backward in time for $t\approx 34.7246$ to obtain $\sigma$ smaller than $10^{-15}$. Lets stress this point again: computing this orbit using the linear approximation of the stable/unstable manifolds would require almost $35$ units of time integration to cross from the unstable to the stable normal bundles. Using the parameterized manifolds no integration is necessary and the entire orbit is represented “locally”.

These results are displayed in figure 8. The green curve is the one obtained using the stable manifold while the trajectory in red correspond to the one using the unstable manifold and the time $t=34.7246$.

4.2 The Circular Restricted Three Body Problem

The circular restricted three body problem (CRTBP) describes the evolution of a massless particle moving in the gravitational field of two other massive bodies, called the primaries. It is assumed that the primaries move in circular orbits about their center of mass. In this work we focus on the case that the third body moves in the plane of the primaries. In co-rotating coordinates the equations of motion are

[TABLE]

with

[TABLE]

Here $\mu$ is the mass ratio of the primary bodies. In the rotating frame:

the center of mass is at the origin, 2. 2.

the motion of the primaries is fixed and they sit on the $x-$axis at $-\mu$ and $1-\mu$.

This choice of coordinates introduces the Coriolis effect, that is the coordinates are non-inertial.

The CRTBP is much studied as an example in Hamiltonian dynamics and celestial mechanics going back to the work of Poincaré. It is a useful model of the motion of a satellite or astroid influenced by a two body system such as Earth/Moon, Sun/Earth, or Sun/Jupiter. It is also one of the simplest $N$-body systems which admits chaotic motions, hence is not integrable. For further discussion of this and more general $N-$body problems, we refer the reader to the books of Belbruno, [2004, 2007]; Gómez et al., [2001]; Jorba and Masdemont, [1999]; Koon et al., [2000]; Meyer et al., [2009]. See also Alessi et al., [2009]; Belbruno et al., [2013]; Belbruno, [1981]; Canalias and Masdemont, [2006]; Font et al., [2009]; Gómez et al., [2004]; Koon et al., [2001]; Llibre et al., [1985]; Llibre and Simó, [1980]; Martínez and Simó, [2014]; Masdemont, [2005]. This list of references constitutes only the barest introduction to the relevant literature, but much more complete discussion is found in the books and papers just cited.

The system has five relative equilibrium points, also referred to as libration points, and we denote these by $\mathcal{L}_{j}$ with $j=1,2,3,4,5$. Three of these lie on the $x$-axis (the co-linear libration points). These are denoted $\mathcal{L}_{1,2,3}$. See Figure 9 for a schematic representation of the problem.

Each of the co-linear libration points are of saddle $\times$ center stability type. By an application of the Lyapunov center theorem, each of the centers gives rise to one parameter families of periodic orbits. These are known as Lyapunov families. We compute invariant manifolds attached to some Lyapunov orbits below. First we transform to a polynomial problem.

Rewriting (18) as a vector field $\dot{u}=\hat{g}(u)$ gives

[TABLE]

This system is non-polynomial and we propose a related polynomial system through the use of automatic differentiation. Following Burgos-García et al., [2017]; Kepley and Mireles James, [2017]; Lessard et al., [2016] and the discussion in Section 2.3 we derive the polynomial vector field

[TABLE]

with the additional (initial condition) constraints

[TABLE]

That is, we take $R$ defined by

[TABLE]

where $\mathbf{u}=(u^{1},u^{2},u^{3},u^{4})$. Then the conjugacy of Equation (7) is satisfied. Moreover, one easily checks that $\pi_{4}(g(R(\mathbf{u})))=\hat{g}(\mathbf{u})$. Of course $R$ is defined only on the complement of the collision set of the CRTBP.

Since periodic orbits in the CRTBP occur in one parameter families parameterized by energy, we fix $L$ and look for a periodic orbit with this half period. Moreover, we use the well known reversible symmetry of the problem to formulate a different boundary condition that will still provide periodic solutions. For example, orbits in a Lyapunov family of one of the co-linear libration points have no velocity in the $x$ direction when they cross the $x$ axis – and the time between the two crossing of the $x$ axis are separated by exactly the half orbit. So that a periodic solution $\gamma$ starting on the $x$ axis with frequency $2L$ will satisfy

[TABLE]

Remark 4.2**.**

Note that the last two conditions of (21) could also be $\gamma_{2}(L)=\gamma_{3}(L)=0$. This would allow us to compute only the half orbit, the symmetry of the problem providing the solution for the second half of the trajectory. The use of Chebyshev expansion would allow such a choice. However, since the periodic orbit is not the final product of our computations, but only an input into the higher order equations we compute full orbits in this work.

For the remaining two conditions, we simply rewrite the initial condition on $v^{5}$ and $v^{6}$ after imposing the symmetry, so that

[TABLE]

Due to the choice of the boundary condition, time translation of the solution does not satisfy this system. Thus we drop the Poincaré condition previously given and obtain a system that is still fully determined and with isolated solutions (since we don’t solve for the period we dispense with one scalar equation). Two possibilities arise from that remark, one could fix $L$ to a given value and still find an orbit as previously mentioned. The other choice is to use $L$ as a variable and obtain a fully determined system by adding a condition in which the energy level is fixed to a chosen constant. Such a choice is necessary for example if we want to compute heteroclinic connecting orbits.

We summarize the discussion in the following lemma.

Lemma 4.3**.**

Let $v^{1}(t),v^{2}(t),v^{3}(t),v^{4}(t)$ be periodic function with same period $\omega$ and such that $\dot{v}^{1}=v^{2}$, $\dot{v}^{3}=v^{4}$. Let $v^{5}(t),v^{6}(t)$ satisfy

[TABLE]

and

[TABLE]

Then $v^{5}(t)$ and $v^{6}(t)$ are periodic with period $\omega$.

Proof 4.4**.**

We first note that

[TABLE]

is a solution of (22). But by unicity of the solution it follows that $v^{5}(t)=w^{5}(t)$. Moreover, we have that $w^{5}(t)$ is periodic with period $\omega$ since $v^{1}(t)$ and $v^{3}(t)$ are periodic with period $\omega$. Thus $v^{5}(t)$ is periodic with period $\omega$, as desired. For $v^{6}(t)$, the proof is similar using the fact that

[TABLE]

is a solution of (23).

To fix the energy we have to use an integral of the CRTBP, namely

[TABLE]

known as the Jacobi integral. In our system of coordinates, this is

[TABLE]

Since $E$ is constant along any orbit of the system, we evaluate it at the endpoint of the first piece of the Chebyshev decomposition. This leads to a new phase condition that can replace the one previously exhibited at equation (16). In terms of the Chebyshev coefficients it is given by

[TABLE]

We now focus on expanding the operator we need to solve in order to obtain the coefficients of the Chebyshev expansion of the $i-$th component of the solution. In the case $k\geq 1$, each case is given by

[TABLE]

The cases $k=0$ are as given in (13) for every $i=2,\ldots,D$. In the case $i=1$, for the periodic orbit we use the condition previously given to rewrite the problem as a zero finding of an operator. Those are given by

[TABLE]

The operator is now completely determined. In order to approximate the solution we truncate the unknowns and the operator a similar way as for the Lorenz system and obtain an operator such that $F:\mathbb{R}^{6mD+1}\to\mathbb{R}^{6mD+1}$ in the case $\alpha=0,1$ and such that $F:\mathbb{R}^{6mD}\to\mathbb{R}^{6mD}$ in the higher dimensional cases. Again, only the search for the periodic orbit requires a good initial guess to obtain the approximation. Following Lessard et al., [2016] we get an initial guess on which we applied Newton’s method. To present results, we regrouped different orbits with the same level of energy.

Figure 10 illustrates two Lyapunov periodic orbits associated with the libration points $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$. Both orbit have energy $E=3.17$ and have been computed using $m=50$ Chebyshev nodes and $D=8$ Chebyshev domains. The periods are $L\approx 1.7122$ for the orbit around $\mathcal{L}_{2}$ and $L\approx 1.4242$ for the orbit around $\mathcal{L}_{1}$. See Figure 11 for an illustration of the numerically computed parameterized invariant manifolds.

We must stress that these computations were preformed for the polynomial equation $\dot{v}=g(v(t))$, where $g$ is given by (20). To show that a point laying in an invariant manifolds in the polynomial problem corresponds to a point on the manifold for the original CRTBP we use the following Theorem. Intuitively this works because of the conjugacy described in Equation (7).

Theorem 4.5**.**

Let $P(t,\sigma)$ be a parameterization of the local stable manifold of a given periodic orbit $v(t)$ satisfying $\dot{v}=g(v(t))$, where $g$ is given by (20). If $x_{0}\in\mathbb{R}^{6}$ is such that $x_{0}=P(t_{0},\sigma_{0})$, then the point $y_{0}\in\mathbb{R}^{4}$ given by the first four component of $x_{0}$ is in the stable set of the corresponding periodic orbit of $\dot{u}=f(u)$, where $f$ is given by (19).

Proof 4.6**.**

Note that using the same remark as in the proof of Lemma (4.3), we obtain from the periodic orbit $v(t)$ a periodic orbit of the CRTBP. That is a periodic solution of $\dot{u}=f(u(t))$. We denote this orbit by $u(t)$ and we have that $v^{i}(t)=u^{i}(t)$ for $i=1,2,3,4$, while

[TABLE]

Let $x(t)$ denote the trajectory obtained by flowing $x_{0}$ forward in time by $\Phi(v,t)$ the flow solution of $\dot{v}=g(v(t))$. So that, using the conjugacy relation (4), we have

[TABLE]

where $\lambda$ is the stable eigenvalue associated to the orbit. Moreover, it follows that

[TABLE]

Thus, by definition of $g$

[TABLE]

and

[TABLE]

Here $C_{1},C_{2}$ are arbitrary constants. But $x_{0}\in W^{s}(v)$, so that for any $\epsilon>0$ there exists $T$ such that for all $t\geq T$, we have that

[TABLE]

This force $C_{1}=0$. Similarly, we have that $C_{2}=0$. Now $x^{5}$ and $x^{6}$ in $g$ are rewritten with the first four component so that $g$ reduce to $f$ and $y(t)=(x^{1}(t),x^{2}(t),x^{3}(t),x^{4}(t))$ satisfies $\dot{y}(t)=f(y(t))$. Moreover, since $x_{0}\in W^{s}(v)$, we have that $y(t)$ is in the stable set of the corresponding orbit $u(t)$, as desired.

4.2.1 Connecting orbits as solutions of boundary value problems

As in Section 4.1.1 we now use the parametrized manifolds to compute connecting orbits between period orbits. In the present section we consider heteroclinic orbits for the CRTBP and we do not find any short connections. Instead, we solve a two point boundary value problem with the manifolds as boundary conditions. This strategy is standard and is for example discussed in detail in Doedel and Friedman, [1989]; Friedman and Doedel, [1993]; Doedel et al., [2008, 2009].

The references just cited obtain boundary conditions by projecting onto the linear approximation of the stable/unstable manifolds given by the stable/unstable normal bundles associated with the periodic orbit. Projecting instead onto high order parameterizations can substantially reduce the integration time and numerically stabilize the problem. The behavior of the connecting orbit on the manifold is then recovered via the flow conjugacy.

Remark 4.7**.**

In this work we compute connecting orbits by numerically integrating the system. However, one could adapt the approach developed in section 3 with $\alpha=0$ to compute any orbit solution of a given boundary value problem expanded as Chebyshev series. This has been done in Lessard and Reinhardt, [2014]; van den Berg et al., [2015]; van den Berg and Sheombarsing, [2016] and even leads to computer assisted proofs. We return to this remark in an upcoming work. We also refer to Arioli, [2002, 2004]; Capiński, [2012]; Wilczak and Zgliczyński, [2003] for further reading about computer assisted proofs in the CRTBP.

We let $P$ and $Q$ denote parameterizations of the stable and unstable manifolds, and seek $(\theta_{u},\sigma_{u})$, $(\theta_{s},\sigma_{s})$ and an integrating time $T$ – or “time of flight” – such that

[TABLE]

where $P_{0}=Q(\theta_{u},\sigma_{u})$. Equation (25) has five unknowns, namely the integrating time and the parameters on both manifold. Since orbits of the CRTBP lie in $\mathbb{R}^{4}$, Equation (25) provides only four equations and we have more unknowns than equations. As a result we cannot expect to isolate a solution. To remedy the situation we simply fix $\sigma_{u}=-1$, removing one of the variables. This corresponds to a choice of boundary components for the local unstable manifold.

Remark 4.8**.**

Recall that the energy is constant along solution curves, and it’s impossible to find a connecting orbit between two periodic solution with different energy level. We avoided this problem by introducing the energy as the phase condition when we solve for the underlying orbit. Recall the definition of $f^{(0)}$ in (24).

Using Newton’s method with the unstable manifold of the orbit on the right in figure 10 and the stable manifold of the orbit on the left we found an approximation of a solution to this problem where

[TABLE]

Both manifolds were computed with $m=50$, $D=8$, $N=50$, $K=5$ and $k_{0}=10$. To find an initial guess on which to apply Newton’s method we integrated $40$ points evenly distributed on the boundary of the unstable manifold and observed that some orbits were potentially intersecting the stable manifold. The connecting orbit and the two manifold are displayed in figure 12.

We also use the conjugacy relation to extend the connecting orbit forward and backward on the manifolds. Integrating in the parameter space until $\sigma_{s}\approx 10^{-15}$ takes

[TABLE]

and for the backward trajectory

[TABLE]

The full trajectory is displayed in figure 13. Note that out of the three pieces of the trajectory, only the one in blue was obtained by numerically integrating the system and this piece required a time step of only $2$ time units. Redoing the computation but projecting onto the linear approximations would result in a time of flight of roughly $30$ time units.

4.3 A Circular Restricted Four Body Problem

We now consider a gravitational problem consisting of three massive bodies (again called the primaries) located at the vertices of an equilateral triangle in the central configuration of Lagrange. These bodies rotate in circular orbits about their common center of mass, all with the same period, rigidly fixing the triangular formation. After changing to a co-rotating frame we are interested the motion of a massless fourth particle moving in the gravitational field of the primaries. In the present work we suppose that the massless particle moves in the plane of the primaries.

It is standard practice to normalize the masses of the primaries so that $0<m_{3}\leq m_{2}\leq m_{1}$, and

[TABLE]

The rotating coordinates are chosen so that the center of mass is at the origin, the largest primary is on the $x$-axis, the $x$-axis cuts the side of the triangle opposite the largest primary, and the smallest primary is in the first quadrant. More explicitly, the primaries are located at positions

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

Define the potential function

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

The equations of motion for the massless particle in the co-rotating coordinates are

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The problem is the subject of many studies beginning with the work of Simó, [1978]. Detailed analysis of the equilibrium solutions and their stability are found in Barros and Leandro, [2014]; Leandro, [2006]; Baltagiannis and Papadakis, 2011b ; Rusu and Santoprete, [2016]; Alvarez-Ramírez and Delgado, [2003], while periodic orbits are studied in Papadakis, 2016a ; Baltagiannis and Papadakis, 2011a ; Papadakis, 2016b ; Burgos-García and Delgado, [2013]; Burgos-García et al., [2017]. The earlier study of Blazevski and Ocampo, [2012] considers stable/unstable manifolds attached to periodic orbits (using the linear approximation given by the stable/unstable normal bundles combined with numerical integration). More complex dynamical behavior such as heteroclinic/homoclinic phenomena and transport are studied numerically in Alvarez-Ramírez and Barrabés, [2015]; Álvarez-Ramírez and Vidal, [2009]; Gidea and Burgos, [2003]; Kepley and Mireles James, [2017]. A Hill’s approximation is derived in Burgos-García and Gidea, [2015]. See also the more theoretical studies of Cheng and She, [2015, 2017]; She and Cheng, [2014].

The CRFBP has $8,9$, or $10$ equilibrium points depending on the values of the mass parameters. These will have either center $\times$ center, center $\times$ saddle, or saddle-focus stability, depending on the ratios of the masses. As in the CRTBP, the center $\times$ saddle equilibria give rise to one parameter families of hyperbolic Lyapunov orbits. We compute local stable/unstable manifolds attached to one of these below.

Using automatic differentiation we derive a related polynomial vector field

[TABLE]

The mapping $R$ used in the automatic differentiation is defined Kepley and Mireles James, [2017], but is similar to the mapping discussed above for the three body case.

The constants terms $\beta,\alpha_{1},\alpha_{2},\alpha_{3}$ are Lagrange multipliers, which are needed to isolate a periodic solution (this time we will not impose any symmetry, hence the boundary condition constraints remain and have to be balanced). The following result, whose proof is found in Burgos-García et al., [2017], explains the relation between the polynomial and non-polynomial problems.

Lemma 4.9**.**

Assume that $\beta,\alpha_{1},\alpha_{2},\alpha_{3},L\in\mathbb{R}$ are fixed constant with $L>0$ and let $\textbf{n},\textbf{p}\in\mathbb{R}^{6}$ be fixed vector. Suppose that $u:[0,2L]\to\mathbb{R}^{9}$ is a periodic solution of $\dot{v}(t)=g(v)$ with $g$ as above and

[TABLE]

and that $u^{7}(t),u^{8}(t),u^{9}(t)>0$ for all $t\in[0,2L]$. Then

$\beta=\alpha_{1}=\alpha_{2}=\alpha_{3}=0$** 2. 2.

the function $\hat{u}:[0,2L]\to\mathbb{R}^{6}$ given by

[TABLE]

is a periodic solution of the four body problem.

Yet, the variables $\alpha_{1},\alpha_{2},\alpha_{3}$ are not necessary in the Chebyshev setting since one can use the following lemma to force the initial value condition on $v^{7},v^{8},v^{9}$ without introducing additional equations. The proof is omitted since it is similar to the proof of Lemma 4.3.

Lemma 4.10**.**

*Let $v^{1}(t),v^{2}(t),v^{3}(t),v^{4}(t),v^{5}(t),v^{6}(t)$ be periodic solution with same period $\omega$ and such that $\dot{v}^{1}=v^{2}$, $\dot{v}^{3}=v^{4}$, $\dot{v}^{5}=v^{6}$. Let $v^{7}(t),v^{8}(t),v^{9}(t)$ satisfy *

[TABLE]

Then $v^{7}(t),v^{8}(t)$ and $v^{9}(t)$ are periodic with period $\omega$.

The extra condition balancing the system is the Poincaré condition which rewrites exactly as in (16). This condition rejects potential time translation of a periodic solution. The other conditions are coming from automatic differentiation and are given by

[TABLE]

The boundary condition for each Chebyshev subdomain being used for this problem define the operators $f^{(i,j)}_{\alpha,0}$ as previously given in (13) and (14). For all $1\leq i\leq D$, the case for $k\geq 1$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Again, we solve recursively the truncated operator to obtain an approximation of the manifold. The case $\alpha=0$ will have the extra variable $\beta$ from Lemma (4.9) and the case $\alpha=1$ will have the eigenvalue as an extra unknown. Note that in this case we fix the frequency $L$ to a constant, it is still possible to find a solution for the same reason as mentioned in the case of the CRTBP. In figure 15, we display an unstable manifold for a planar Lyapunov orbit about $\mathcal{L}_{8}$. The computations were done with $D=4$ and $N=60$. We used masses

[TABLE]

We remark that, for the purposes of the present demonstration, we could have taken $m_{1},m_{2},m_{3}$ to have more or less any values. Nevertheless the particular choice of mass values was suggested in Burgos-García, [2017]. The values of $m_{1},m_{2}$ correspond respectively to the mass of the primary star in the binary system Epsilon Reticuli (Henry Draper Catalogue number $27442$), and the mass of an extra solar planet discovered orbiting HD $27442$ in the year $2000$. The mass ratio of these bodies is such that the system could form a CRFBP with a third trojan object. Our value of $m_{3}$ corresponds to that of a supposed Saturn like planet forming an equilateral triangle with $m_{1}$ and $m_{2}$. We refer the reader to Schwarz et al., [2009] for more extensive discussion of exoplanets.

We consider a Lyapunov orbit near the heaviest mass (an orbit about the libration point $\mathcal{L}_{8}$), and for this choice a uniform mesh was not suitable. That is: it gets harder to obtain an accurate approximation as the orbit approach one of the heavy body. Thus, we took

[TABLE]

Note that these sum to $1$, so that the integrating time for the orbit is preserved. To understand why the accuracy is affected as the orbit approach a body, recall that the variables arising from automatic differentiation are inversely proportional to the distance between the object and the corresponding primary, thus provoking a considerable change in the amplitude in the additional variables.

4.3.1 Connecting Orbits

We now consider a homoclinic connection, and hence there are no energy considerations. For the planar Lyapunov periodic orbit displayed in figure 15 we computed both the stable and unstable manifold and apply a similar BVP approach as for the CRTBP. Again, the connecting orbit starts on the boundary of the unstable manifold and ends on the boundary of the stable manifold. The conjugacy relation is used to compute the asymptotic behavior without phase space integration. In this case both eigenvalues have the same value with opposite signs, so the integrating time forward or backward needed is the same.

We remark that to get from the boundary of the invariant manifold to $\sigma\approx 10^{-8}$ close to the periodic orbit (so that the error in the linear approximation is on the order of machine epsilon) one would need to integrate for roughly $t=342$ time units. So this problem illustrates starkly the utility of using the local parameterizations to absorb such a substantial portion of the homoclinic orbit.

Both manifolds are computed using $N=50$, $D=4$, $m=50$, $K=2$ and $k_{0}=10$. We computed one connecting orbit for each component of the local manifold boundaries. In the case of $\sigma_{u}=1$, the coordinates for the connecting orbit are

[TABLE]

For the case of $\sigma_{u}=-1$, the coordinates are

[TABLE]

The sign of the value $\sigma$ is affected by the choice of the eigenvector, i.e. the sign determines the polarity of the embedding. In this case we picked the eigenvectors so that the boundaries have the same sign when they lay on the same side of the orbit in the choice of coordinates displayed.

The reason the dynamics on the parameterized manifolds are so slow in this example is that the Floquet exponent $\lambda\approx\pm 0.0538$ is much closer to zero than in any previous example. The connecting orbits are displayed in figure 16 along with both manifolds. The extension of the first orbit using the conjugacy relation is displayed in figure 17. In both cases the coordinates displayed are $(x,y,\dot{y})$.

5 Conclusion

The methods of the present work facilitate accurate computation of local stable/unstable manifolds attached to periodic orbits of vector fields in a large neighborhood of the orbit itself. The computations exploit Chebyshev expansions, so that domain decomposition can be used to improve the accuracy of the parameterization (compared to a Fourier-Taylor expansion) without necessarily increasing the total number spectral modes used. The results approximate the local manifolds in relatively large regions of phase space. The method is based on solving an invariance equation, so that the computations are equipped with a convenient notion of defect/a-posteriori error.

The method is non-perturbative. So even though manifolds for Lyapunov orbits were computed in the CRTB and CRFB problems, the calculations do not use the fact that we were near an equilibrium. Since the parameterization method is based on finding a zero of an invariance equation, and since we use a Newton scheme to compute the numerical solution, it would be natural to develop numerical continuation methods for the manifold computations. Continuation in frequency, energy, or other system parameter would be natural. In this case the manifolds, and even the connecting orbits, do not have to be recomputed from scratch as parameters are changed. Rather, the old orbits/manifolds can be used as the initial guess for the Newton method at the new parameters. Utilizing a predictor/corrector scheme would also be natural.

An interesting topic for future work would be to compare the techniques developed in the present work with other techniques for computing high order expansions of local invariant manifolds attached periodic orbits. For example, invariant objects similar to those discussed in Section 4.2 have been computed by a number of authors using methods based on Lindstedt-Poincaré series or using high order normal forms Masdemont, [2011]; Delshams et al., [2008]; Masdemont, [2005]; Gidea and Masdemont, [2007]; Jorba and Villanueva, [1998]. One of the main differences between the methods of the works just cited, and the methods of the present work, is that both Lindstedt-Poincaré and normal form methods develop expansions valid in a full neighborhood of the periodic orbit. The stable/unstable manifolds are then obtained as suitable zero sections. That is: the number of variables used in the expansion is equal to the dimension of the phase space rather than the dimension of the underlying invariant object.

Computing an expansion of a full neighborhood of the periodic orbit is important for many applications. For example when designing a “fly-by” mission one wants to find trajectories which approach the periodic orbit along (but not on) the stable manifold, and then move away after a finite time along (but not on) the unstable manifold. On the other hand if one is primarily interested in heteroclinic and homoclinic connections then computing on a full neighborhood is much more expensive then just parameterizing the manifold as in the present work. For example the parameterizations in Section 4.2 were computed using $N=10$ Taylor order. Since the parameterizations were expanded only in the stable or unstable direction, this requires computing only ten unknown Taylor coefficient (each of which is a periodic function). If one uses instead Lindstedt-Poincaré or normal form methods, then it is necessary to expand in one angle variable and three polynomial variables, and a polynomial of order $N=10$ in three variables has $286$ unknown coefficients (again, each of these is a periodic function). This back of the envelope comparison illustrates the advantage of parameterizing only the desired manifold, and not the full neighborhood when the particular application allows.

Another interesting possibility for future improvement is to study more carefully the effects of a non-uniform subdivision strategy for choosing the Chebyshev domains. In figure 18, we computed the norm of each sequence of Chebyshev expansion in the case of the stable manifold of the orbit AB displayed in figure 4

[TABLE]

for all $\alpha=0,\ldots,100$ and $i=1,\ldots,50$. The scale of the eigenvector was chosen so that the norm of the last Taylor dimension is below machine precision in every subdomain. However, one can see from the figure that in this case several component were reaching this magnitude much earlier than some other ones. Such differences arise from the fact that the mesh was uniform in this case. One way to obtain bigger manifold without increasing the number of modes would be to use mesh adaptation, as in van den Berg and Sheombarsing, [2016]. We also remark that there is a possibility that the computations could be sped up by pre-computing the Floquet normal form as in Castelli et al., [2015]. The Floquet normal form would have to be discretized using Chebyshev rather than Fourier series, and we have not yet explored this possibility.

\nonumsection

Acknowledgments The authors wish to thank Jaime Burgos-García and J.P. Lessard for helpful conversations, and to thank an anonymous referee for carefully reading the original submission and making a number of helpful comments and suggestions. The final version of the manuscript is improved thanks to these efforts. Both Maxime Murray and J.D. Mireles James were partially supported by NSF grant DMS-1700154, and by the Alfred P. Sloan Foundation grant G-2016-7320.