Efficient algorithms for invariant tori and their manifolds

# Computation of whiskered invariant tori and their associated manifolds: new fast algorithms

Gemma Huguet Centre de Recerca Matemàtica, Apartat 50, 08193 Bellaterra (Barcelona), Spain Rafael de la Llave Department of Mathematics, The University of Texas at Austin, Austin, TX 78712  and  Yannick Sire Université Paul Cézanne, Laboratoire LATP UMR 6632, Marseille, France
###### Abstract.

In this paper we present efficient algorithms for the computation of several invariant objects for Hamiltonian dynamics. More precisely, we consider KAM tori (i.e diffeomorphic copies of the torus such that the motion on them is conjugated to a rigid rotation) both Lagrangian tori (of maximal dimension) and whiskered tori (i.e. tori with hyperbolic directions which, together with the tangents to the torus and the symplectic conjugates span the whole tangent space). In the case of whiskered tori, we also present algorithms to compute the invariant splitting and the invariant manifolds associated to the splitting. We present them both for the case of discrete time and for differential equations.

The algorithms for tori are based on a Newton method to solve an appropriately chosen functional equation that expresses invariance. Among their features we highlight:

• The algorithms are efficient: if we discretize the objects by elements, one step of the Newton method requires only storage and operations. Furthermore, if the object we consider is of dimension , we only need to compute functions of variables, independently of what is the dimension of the phase space.

• The algorithms do not require that the system is presented in action-angle variables nor that it is close to integrable.

• The algorithms are backed up by rigorous a-posteriori bounds which state that if the equations are solved with a small residual and some explicitly computable condition numbers are not too big, then, there is a true solution which is close to the computed one.

• The algorithms apply both to primary (i.e non-contractible) and secondary tori (i.e. contractible to a torus of lower dimension, such as islands). They also apply to whiskered tori.

The algorithms for invariant splittings are based on computations of proyections (rather than in graph transforms). The computations of invariant manifolds are also efficient in the sense indicated before.

The algorithms we present have already been implemented. We will report on the technicalities of the implementation and the results of running them elsewhere.

###### Key words and phrases:
Quasi-periodic solutions, whiskered KAM tori, whiskers, quasi-periodic cocycles, numerical computation
###### 2000 Mathematics Subject Classification:
Primary: 70K43, Secondary: 37J40

## 1. Introduction

The goal of this paper is to present efficient algorithms to compute very accurately several objects of interest in Hamiltonian dynamical systems (both discrete-time dynamical systems and differential equations). More precisely, we present algorithms to compute:

• Lagrangian KAM tori.

• Whiskered KAM tori.

• The invariant bundles of the whiskered tori.

• The stable and unstable manifolds of the whiskered tori.

The algorithms are very different. For example, the algorithms for tori require the use of small divisors and symplectic geometry and the algorithms for invariant bundles and invariant manifolds rely on the theory of normal hyperbolicity and dichotomies. The computation of whikered tori has to combine both.

We recall that KAM tori are manifolds diffeomorphic to a torus which are invariant for a map or flow, on which the motion of the system is conjugate to a rotation. As we will see later, this is also equivalent to quasi-periodic solutions. The tori are called Lagrangian when they are Lagrangian manifolds, which in our case amounts just to the fact that the tori have a dimension equal to the number of degrees of freedom of the system. The tori are called whiskered when the linearized equation has directions that decrease exponentially either in the future (stable) or in the past (unstable) and these directions together with the tangent to the torus and its symplectic conjugate span the whole tangent space. These invariant spaces for the linearization have non-linear analogues, namely invariant manifolds. It has been recognized since [Arn64] that whiskered tori and their invariant manifolds are very interesting landmarks that organize the long-term behavior of many systems.

The algorithms we present are based on running an efficient Newton method to solve a functional equation, which expresses the dynamical properties above. What we mean by efficient is that if we discretize the problem using Fourier coefficients, we require storage and only operations for the Newton step. Since the functions we are considering are analytic, we see that the truncation error is where is the dimension of the object. Note that, in contrast, a straightforward implementation of a Newton method would require to use storage – to store the linearization matrix and its inverse – and operations to invert.

In practical applications, using the algorithms described in this paper, computing with several million coefficients becomes quite practical in a typical desktop computer of today. Implementation details and the results of several runs will be discussed in another companion paper [HdlLS09]. Given the characteristics of today’s computers, savings in storage space are more crucial than savings in operations for these problems.

The algorithms we present here are inspired by the rigorous results of [LGJV05] – for KAM tori – and [FdlLS09b, FdlLS09a] –for whiskered tori. The algorithms to compute the stable and unstable manifolds had not been previously discussed. The rigorous results of the above papers are also based on a Newton method applied to the same functional equation that we consider here. Of course, going from a mathematical treatment to a practical algorithm requires significantly many more details. In particular, the algorithms to compute invariant splittings and invariant manifolds are different from those in the above references. This paper discusses these algorithmic issues.

The results of the papers [LGJV05, FdlLS09a] give a justification of the algorithms for tori and splittings presented here. The theorems in [LGJV05, FdlLS09a], have been formulated in an a-posteriori way, i.e. the theorems assert that if we have a function which solves approximately the invariance equation very accurately (e.g. the outcome of a successful run of the algorithms) and which also satisfies some explicit non-degeneracy conditions, then, we can conclude that there is a true solution which is close to the computed solution. Hence, by supplementing our calculations with the (very simple) computations of the non-degeneracy conditions (they play a role very similar to the condition numbers common in numerical analysis), we can be sure that the computation that we are performing is meaningful. This allows to compute with confidence even close to the limit of validity of the KAM theorem (a rather delicate boundary since the smooth KAM tori do not disappear completely but rather morph into Cantor sets).

Since the papers [LGJV05, FdlLS09a] contain estimates, in the present paper, we will only discuss the algorithmic issues. For example, we will detail how solutions of equations (whose existence was shown in the above papers) can be computed with small requirements of storage and small operation count. Note that different algorithms of a same mathematical operation can have widely different operation counts and storage requirements. (See, for example, the discussion in [Knu97] on the different algorithms to multiply matrices, polynomials, etc.) On the other hand, we will not include some implementation issues (methods of storage of arrays, ordering of loops, precision, etc.) needed to obtain actual results in a real computer. They will be given in another paper together with experimental results obtained by running the algorithms.

One remarkable feature of the algorithms presented here is that they do not require the system to be close to integrable. We only need a good initial guess for the Newton method. Typically, one uses a continuation method starting from an integrable case, where solutions can be computed analytically. However, in the case of secondary KAM tori, which do not exist in the integrable case, one can use, for instance, Lindstedt series, variational methods or approximation by periodic orbits to obtain an initial guess.

As for the algorithms to compute invariant splittings, we depart from the standard mathematical methods (most of the time based on graph transforms) and we have found more efficient to device an equation for the invariant projections. We also found some acceleration of convergence methods that give superexponential convergence. They are based on fast algorithms to solve cohomology equations which could be of independent interest (See Appendix A).

The algorithms to compute invariant manifolds are based on the parameterization method [CFL03a, CFL05]. Compared to standard methods such as the graph transform it has the advantage that to compute geometric objects of dimension , we only need to compute with functions of dimension . In contrast with [CFL03a, CFL05], which was based on contractive iterations, our method is based on a Newton iteration which we also implement without requiring large matrices and requiring only operations.

### An overview of the method

The numerical method we use is based on the parameterization methods introduced in [CFL03a, CFL03b]. In this section, we provide a sketch of the issues, postponing some important details. For brevity, we make the presentation for maps only, even if a similar sketch can be made for flows.

#### Invariant tori

We observe that if is a map and we can find an embedding in which the motion on the torus is a rotation , it should satisfy the equation

 F(K(θ))−K(θ+ω)=0. (1)

Given an approximate solution of (1), i.e.

 F(K(θ))−K(θ+ω)=E(θ),

the Newton method aims to find solving

 DF(K(θ))Δ(θ)−Δ(θ+ω)=−E(θ), (2)

so that will be a much more approximate solution.

The main idea of the Newton method is that, using the decomposition into invariant subspaces, one can decompose (2) into three components

 DF(K(θ))Δs(θ)−Δs(θ+ω)=−Es(θ)DF(K(θ))Δu(θ)−Δu(θ+ω)=−Eu(θ)DF(K(θ))Δc(θ)−Δc(θ+ω)=−Ec(θ) (3)

where the refer to the stable, unstable components and refers to the component along the tangent to the torus and its symplectic conjugate. For Lagrangian tori, only the part appears in the equations.

The algorithm requires:

• Efficient methods to evaluate the LHS of (1).

• Efficient methods to compute the splitting.

• Efficient methods to solve the equations (3).

As we will see in Section 3.6, to evaluate (1), it is efficient to use both a Fourier representation (which makes easy to evaluate ) and a real representation which makes easy to evaluate . Of course, both of them are linked through the Fast Fourier Transform (FFT from now on).

The methods to compute the splitting are discussed in Section 4.3. More precisely, we present a numerical procedure to compute the projections on the linear stable/unstables subspaces based on a Newton method. In [HdlLS09], we present an alternative procedure for the computation of the projections based on the calculation of invariant bundles for cocycles. Indeed, these algorithms require the computation of the projections over the linear subspaces of the linear cocycle.

The solution of the hyperbolic components in equation (3) is discussed in Section 4.1 and Appendix A. Indeed, equations of this form appear as well in the calculation of the invariant splitting discussed in Section 4.3. A first method is based on an acceleration of the fixed point iteration (Appendix A.1). We note that to obtain superexponential convergence for the solution of (3), we need to use both the Fourier representation and the real space representation.

In the case that the bundles are one-dimensional, there is yet another algorithm, which is even faster than the previous ones (see Appendix A.2). The algorithms are discussed for maps, and they do not have an easy analog for flows except by passing for the integration of differential equations. We think that this is one case where working with time- maps is advantageous.

The most challenging step is the solution of the center component of (3). This depends on cancelations which use the symplectic structure, involves small divisors and requires that certain obstructions vanish. Using several geometric identities that take advantage of the fact that the map is symplectic (see Section 4.2), the solution of (3) in the center direction is reduced to solving the following equation for given ,

 φ(θ)−φ(θ+ω)=η(θ). (4)

Equation (4) can be readily solved using Fourier coefficients provided that (and that is sufficiently irrational). The solution is unique up to addition of a constant.

The existence of obstructions – which are finite dimensional – is one of the main complications of the problem. It is possible to show that, when the map is exact symplectic, the obstructions for the solution are . An alternative method to deal with these obstructions is to add some new – finite dimensional – unknowns and solve, instead of (1), the equation

 F(K(θ)−K(θ+ω)+G(θ)λ=0

where is an explicit function. Even if is kept through the iteration involving approximate solutions, it can be shown that, if the map is exact symplectic, we have . This counterterm approach also helps to weaken non-degeneracy assumptions.

A minor issue that we omit in this preliminary discussion is that the solutions of (1) are not unique. If is a solution, defined by is also a solution for any . This can be easily solved by taking an appropriate normalization that fixes the origin of coordinates in the torus. In [FdlLS09a] it is shown that this is the only non-uniqueness phenomenon of the equation. Furthermore, this local uniqueness property allows to deduce results for vector fields from the results for maps.

It is important to remark that the algorithms that we will present can compute in a unified way both primary and secondary tori. We recall here that secondary tori are invariant tori which are contractible to a torus of lower dimension, whereas this is not the case for primary tori. The tori which appear in integrable systems in action-angle variables are always primary. In quasi-integrable systems, the tori which appear through Lindstedt series or other perturbative expansions starting from those of the integrable system are always primary. Secondary tori, however, are generated by resonances. In numerical explorations, secondary tori are very prominent features that have been called “islands”. In [HL00], one can find arguments showing that these solutions are very abundant in systems of coupled oscillators. As an example of the importance of secondary tori we will mention that in the recent paper [DLS06] they constituted the essential object to overcome the “large gap problem” and prove the existence of diffusion. In [DH09], one can find a detailed analysis of secondary tori.

In this paper, we will mainly discuss algorithms for systems with dynamics described by diffeomorphisms. For systems described through vector fields, we note that, taking time maps or, more efficiently, surfaces of section, we can reduce the problem with vector fields to a problem with diffeomorphisms. However, in some practical applications, it may be convenient to have a direct treatment of the system described by vector fields. For this reason, we have included the invariance equations for flows, in parallel with the invariance equations for maps and we have left for the Appendix the algorithms that are specially designed for flows.

#### Invariant manifolds attached to invariant tori

When the torus is whiskered, it has invariant manifolds attached to it. For simplicity, in this presentation we will discuss the case of one dimensional directions – even if the torus can be of higher dimension.

We use again a parameterization method. Consider an embedding in which the motion on the torus is a rotation and the motion on the stable (unstable) whisker consists of a contraction (expansion) at rate , it should satisfy the invariance equation

 F(W(θ,s))−W(θ+ω,μs)=0. (5)

Again, the key point is that taking advantage of the geometry of the problem we can devise algorithms which implement a Newton step to solve equation (5) without having to store—and much less invert—a large matrix. We first discuss the so-called order by order method, which serves as a comparison with more efficient methods based on the reducibility. Although they are based on the same idea as for the case of tori, they have not been introduced previously and constitute one of the main novelties of this paper. We present algorithms that compute at the same time the torus and the whiskers and algorithms that given a torus and the linear space compute the invariant manifold tangent to it. It is clearly possible to extend the method to compute stable and unstable manifolds in general dimensions (or even non-resonant bundles). To avoid increasing the length of this paper and since higher dimensional examples are harder numerically, we postpone this to a future paper.

### Some remarks on the literature

Invariant tori in Hamiltonian dynamics have been recognized as important landmarks in Hamiltonian dynamics. In the case of whiskered tori, their manifolds have also been crucial for the study of Arnold diffusion.

Since the mathematical literature is so vast, we cannot hope to summarize it here. We refer to the rather extensive references of [Lla01] for Lagrangian tori and those of [FdlLS09a] for whiskered tori. We will just briefly mention that [Gra74, Zeh76] the earliest references on whiskered theory, as well as most of the later references, are based on transformation theory, that is making changes of variables that reduce the perturbed Hamiltonian to a simple form which obviously presents the invariant torus. From the point of view of numerics, this has the disadvantage that transformations are very hard to implement.

The numerical literature is not as broad as the rigorous one, but it is still quite extensive. The papers closest to our problems are [HL06c, HL06b, HL07], which consider tori of systems under quasi-periodic perturbation. These papers also contain a rather wide bibliography on papers devoted to numerical computation of invariant circles. Among the papers not included in the references of the papers above because these appeared later, we mention [CdlL09], which presents other algorithms which apply to variational problems (even if they do not have a Hamiltonian interpretation). Another fast method is the “fractional iteration method” [Sim00]. Note that the problems considerd in [HL06c, HL06b, HL07] do not involve center directions (and hence, do not deal with small divisors) and that the frequency and one of the coordinates of the torus is given by the external perturbation. The methods of [HL06c, HL06b, HL07] work even if the system is not symplectic (even if they can take advantage of the symplectic structure).

The papers [JO05, JO09] present and implement calculations of reducible tori. This includes tori with normally elliptic directions. The use of reducibility indeed leads to very fast Newton steps, but it still requires the storage of a large matrix for the changes of variables. As seen in the examples in [HL07, HL06a], reducibility may fail in a codimension set (a Cantor set of codimension manifolds for elliptic tori in Hamiltonian systems). For these reasons, we will not discuss methods based on reducibility in this paper (even if it is a useful and practical tool) and just refer to the references just mentioned. Indeed, thanks to hyperbolicity, reducibility is not needed in the present paper.

The paper is organized as follows. In Section 2 we summarize the notions of mechanics and symplectic geometry we will use. In Section 3 we formulate the invariance equations for the objects of interest (invariant tori, invariant bundles and invariant manifolds) and we will present some generalities about the numerical algorithms.

Algorithms for whiskered tori are discussed in Section 4. In particular, we discuss how to compute the decomposition (3) of the linearized equation (2), and how to solve efficiently each equation in (3).

In Section 5 we discuss fast algorithms to compute rank-1 (un)stable manifolds of whiskered tori. More precisely, we present an efficient Newton method to solve equation (5).

In Appendix A one can find the fast algorithms to solve cohomology equations with non-constant coefficients that will be used in the computation of the splitting (3) as well as to solve the hyperbolic components of equations (3). In Appendices B-E, one can find the algorithms specially designed for flows, analogous to the ones for maps.

## 2. Setup and conventions

We will be working with systems defined on an Euclidean phase space endowed with a symplectic structure. The phase space under consideration will be

 M⊂R2d−ℓ×Tℓ.

We do not assume that the coordinates in the phase space are action-angle variables. Indeed, there are several systems (even quasi-integrable ones) which are very smooth in Cartesian coordinates but less smooth in action-angle variables (e.g., neighborhoods of elliptic fixed points [FGB98, GFB98], hydrogen atoms in crossed electric and magnetic fields [RC95, RC97] and several problems in celestial mechanics [CC07]).

We will assume that the Euclidean manifold is endowed with an exact symplectic structure (for some one-form ) and we have

 Ωz(u,v)=⟨u,J(z)v⟩,

where denotes the inner product on the tangent space of and is a skew-symmetric matrix.

An important particular case is when induces an almost-complex structure, i.e.

 J2=−Id. (6)

Most of our calculations do not need this assumption. One important case, where the identity (6) is not satisfied, is when is a symplectic structure on surfaces of section chosen arbitrarily in the energy surface or when is the symplectic form expressed in symplectic polar coordinates near an elliptic fixed point. When (6) holds, some calculations can be made faster.

As previously mentioned, we will be considering systems described either by diffeomorphisms or by vector-fields.

### 2.1. Systems described by diffeomorphisms

We will consider maps which are not only symplectic (i.e. ) but exact symplectic, that is

 F∗α=α+dP,

for some smooth function , called the primitive function.

We will also need Diophantine properties for the frequencies of the torus. For the case of maps, the useful notion of a Diophantine frequency is:

 D(ν,τ)={ω∈Rℓ∣∣|ω⋅k−n|−1≤ν|k|τ ∀ k∈Zℓ−{0}, n∈Z},ν>ℓ.

### 2.2. Systems described by vector fields

We will assume that the system is described by a globally Hamiltonian vector-field , that is

 X=J∇H,

where is a globally defined function on .

In the case of flows, the appropriate notion of Diophantine numbers is:

 Daff(ν,τ)={ω∈Rℓ∣∣|ω⋅k|−1≤ν|k|τ ∀ k∈Zℓ−{0}},ν≥ℓ−1
###### Remark 1.

It is well known that for non-Diophantine frequencies substantially complicated behavior can appear [Her92, FKW01]. Observing convincingly Liouvillian behaviors seems a very ambitious challenge for numerical exploration.

## 3. Equations for invariance

In this section, we discuss the functional equations for the objects of interest, that is, the invariant tori and the associated whiskers. These functional equations, which describe the invariance of the objects under consideration, are the cornerstone of the algorithms.

### 3.1. Functional equations for whiskered invariant tori for diffeomorphisms

At least at the formal level, it is natural to search quasi-periodic solutions with frequency (independent over the integers) under the form of Fourier series

 x(n)=∑k∈Zℓ^xke2πik⋅ωn , (7)

where and .

We allow some components of in (7) to be angles. In that case, it suffices to take some of the components of modulo 1.

It is then natural to describe a quasi-periodic function using the so-called “hull” function defined by

 K(θ)=∑k∈Zℓ^xke2πik⋅θ,

so that we can write

 x(n)=K(nω).

The geometric interpretation of the hull function is that it gives an embedding from into the phase space. In our applications, the embedding will actually be an immersion.

It is clear that quasi-periodic functions will be orbits for a map if and only if the hull function satisfies:

 F∘K−K∘Tω=0, (8)

where denotes a rigid rotation

 Tω(θ)=θ+ω. (9)

A modification of the invariance equations (8) which we will be important for our purpose consists in considering

 F∘K−K∘Tω−(J(K0)−1DK0)∘Tωλ=0, (10)

where the unknowns are now (as before) and . Here, denotes a given approximate (in a suitable sense which will be given below) solution of the equation (8).

It has been shown in [FdlLS09b, FdlLS09a] (the vanishing lemma that, for exact symplectic maps, if satisfy the equation (10) with close to , then at the end of the iteration of the Newton method, we have and, therefore, is a solution of the invariance equation (8). In other words, the formulations (10) and (8) are equivalent. Of course, for approximate solutions of the invariance equation (8), there is no reason why should vanish and it is numerically advantageous to solve the equation with more variables.

The advantage of equation (10) for numerical calculations is that, at the initial stages of the method, when the error in the invariance equation is large, it is not easy to ensure that certain compatibility conditions (see (20) in Section 3.7), are satisfied approximately, so that the standard Newton method has problems proceeding. On the other hand, we can always proceed by adjusting the . This is particularly important for the case of secondary tori that we will discuss in Section 3.4. We also note that this procedure makes possible to deal with tori when the twist condition degenerates.

The equations (8) and (10) will be the centerpiece of our treatment. We will discretize them using Fourier series and study numerical methods to solve the discretized equations.

It is important to remark that there are a posteriori theorems (see [FdlLS09b, FdlLS09a]) for equations (8), (10) (as well as their analogous for flows (11), (13) ). That is, theorems that ensure that given a function that satisfies (8), (10) up to a small error and that, at the same time, satisfies some non-degeneracy conditions (which are given quite explicitly), then there is a true solution close to the computed one. Hence, if we monitor the non-degeneracy conditions, we can be sure that the computed solutions correspond to some real effects and are not spurious solutions.

###### Remark 2.

Notice that for whiskered tori the dimension of the torus is smaller than half the dimension of the phase space . Hence, the algorithms presented here have the advantage that they look for a function of variables to compute invariant objects of dimension . This is important because the cost of handling functions grows exponentially fast with the number of variables. Indeed, to discretize a function of variables in a grid of side into , one needs to store real values.

###### Remark 3.

Equations (8) and (10) do not have unique solutions. Observe that if is a solution, for any , is also a solution. In [FdlLS09a], it is shown that, in many circumstances, this is the only non uniqueness phenomenon in a sufficiently small neighborhood of . Hence, it is easy to get rid of it by imposing some normalization. See Section 3.5.2 for a discussion on this issue.

### 3.2. Functional equations for whiskered invariant tori for vector-fields

In this case, one can write

 x(t)=∑k∈Zℓ^xke2πik⋅ωt

where , and then the hull function is defined by

 x(t)=K(ωt).

The invariance equation for flows is written:

 ∂ωK−X∘K=0, (11)

where denotes the derivative in direction

 ∂ω=ℓ∑k=1ωk∂θk. (12)

The modification of (11) incorporating a counterterm is:

 ∂ωK−X∘K−J(K0)−1(DX∘K0)λ=0, (13)

where is a given embedding satisfying some non-degeneracy conditions.

###### Remark 4.

Recall that, taking time maps, one can reduce the problem of vector fields to the problem of diffeomorphisms. Furthermore, since autonomous Hamiltonian systems preserve energy, we can take a surface of section and deal with the return map. This reduces by the dimension of the phase space and the parameterization of the torus requires variable less. In practice, it is much more efficient to use a numerical integrator to compute the point of intersection with the surface of section than to deal with functions of one more variable and with two more components.

### 3.3. Some global topological considerations

In our context, both the domain and the range of have topology. As a consequence, there will be some topological considerations in the way that the torus gets embedded in the phase space. More explicitly, the angle variables of can get wrapped around in different ways in the phase space.

A concise way of characterizing the topology of the embedding is to consider the lift of to the universal cover, i.e.

 ˆK:Rℓ→R2d−ℓ×Rℓ,

in such a way that is obtained from by identifying variables in the domain and in the range that differ by an integer.

It is therefore clear that

 ˆKp(θ+e)=ˆKp(θ),\@@LTX@noalign\vskip6.0pt\omitˆKq(θ+e)=ˆKq(θ)+I(e),

where denote the projections of the lift on the and coordinates of . It is easy to see that is a linear function of , namely

 I(e)i=1,…,ℓ=(ℓ∑j=1Iijej)i=1,…,ℓ (14)

with .

We note that if a function satisfies

 ˆKq(θ+e)=ˆKq(θ)+I(e) ,

the function

 ˜Kq(θ)≡ˆKq(θ)−I(θ) (15)

is periodic. The numerical methods will always be based on studying the periodic functions , but we will not emphasize this unless it can lead to confusion.

Of course, the integer valued matrix remains constant if we modify the embedding slightly. Hence, it remains constant under continuous deformation. For example, in the integrable case with , invariant tori satisfy , so that we have . Hence, all the invariant tori which can be continued from tori of the integrable system will also have .

### 3.4. Secondary tori

One can produce other -dimensional tori for which the range of is of dimension less than . These tori are known as secondary tori. It is easy to see that if we can contract to a diffeomorphic copy of . Even in the case of maximal tori , one can have contractible directions. The most famous example of this phenomenon are the “islands” generated in twist maps around resonances.

Secondary tori do not exist in the integrable system and they cannot be even continuously deformed into some of the tori presented in the integrable system. This is often described informally as saying that the secondary tori are generated by the resonances.

Perturbative proofs of existence of secondary tori are done in [LW04] and in [DLS06] and in more detail in [DH09]. In [Dua94] one can find rigorous results showing that these islands have to be rather abundant (in different precise meanings) in many classes of 2D-maps. In particular, for standard-like maps, secondary tori appear at arbitrarily large values of the parameter.

In [HL00], there are heuristic arguments and numerical simulations arguing that in systems of coupled oscillators, the tori with contractible directions are much more abundant than the invariant tori which can be continued from the integrable limit.

In view of these reasons, we will pay special attention to the computation of secondary tori.

We want to emphasize on some features of the method presented here, which are crucial for the computation of secondary tori:

• The method does not require neither the system to be close to integrable nor to be written in action-angle variables.

• The modification of the invariance equations (8) and (11) allows to adjust some global averages required to solve the Newton equations (see [FdlLS09a]).

• The periodicity of the function can be adjusted by the matrix introduced in (3.3). Hence, the rank of the matrix has to be chosen according to the number of contractible directions.

### 3.5. Equations for the invariant whiskers

Invariant tori with may have associated invariant bundles and whiskers. We are interested in computing the invariant manifolds which contain the torus and are tangent to the invariant bundles of the linearization around the torus. This includes the stable and unstable manifolds but also invariant manifolds associated to other invariant bundles of the linearization, such as the slow manifolds, associated to the less contracting directions.

Using the parameterization method, it is natural to develop algorithms for invariant manifolds tangent to invariant sub-bundles that satisfy a non-resonance condition (see [CFL03a]). This includes as particular cases, the stable/unstable manifolds, the strong stable and strong unstable ones as well as some other slow manifolds satisfying some non-resonance conditions.

To avoid lengthening the paper, we restrict in this paper just to the one-dimensional manifolds (see Section 5), where we do not need to deal with resonances as it is the case in higher dimensions. We think that, considering this particular case, we can state in a more clear and simpler way the main idea behind the algorithms. We will come back to the study of higher dimensional manifolds in future work.

#### 3.5.1. Invariant manifolds of rank 1

We once again use a parameterization to describe the whiskers. This amounts to finding a solution of the equations of motion under the form

 u(n)=W(ωn,μns)

in the discrete time case and

 u(t)=W(ωt,seμt)

in the continuous time case, where and . The function has then to satisfy the following invariance equations

 F(W(θ,s))=W(θ+ω,μs),\@@LTX@noalign\vskip6.0pt\omit∂ωW(θ,s)+μs∂∂sW(θ,s)=(X∘W)(θ,s),

for the case of maps and flows, respectively.

Note that equations (3.5.1) imply that in variables the motion on the torus consists of a rigid rotation of frequency whereas the motion on the whiskers consists of a contraction (or an expansion) by a constant ( in the case of flows). We call contractive the situation for maps (or for flows). We call expansive the case when for maps (or for flows). Note that if satisfies (3.5.1) then is a solution of the invariance equations (8) or (11).

As in the case of invariant tori, it will be convenient to consider the following modified invariance equations

 F(W(θ,s))=W(θ+ω,μs)+(J(K0)−1DK0)∘Tωλ,\@@LTX@noalign\vskip6.0pt\omit∂ωW(θ,s)+μs∂∂sW(θ,s)=(X∘W)(θ,s)+J(K0)−1(DX∘K0)λ,

where is, as before, a given approximate solution of the equations (8) and (11), respectively.

#### 3.5.2. Uniqueness of solutions of the invariance equation for whiskers

The solutions of equations (3.5.1) are not unique. Indeed, if is a solution, for any , , we have that is also a solution. This non-uniqueness of the problem can be removed by supplementing the invariance equation with a normalization condition.

Some suitable normalization conditions (in the case of maps) that make the solutions unique are

 ∫TℓW(θ,0)−θ=0,DF(W(θ,0))D2W(θ,0)=μD2W(θ,0),||DW(⋅,0)||=ρ (16)

where denotes the derivative with respect to the second argument, is any arbitratrily chosen number and stands for a suitable norm.

The fact that the solutions of (8) supplemented by (16) are locally unique is proved in [FdlLS09a]. In this paper, we will see that these normalizations uniquely determine the Taylor expansions (in ) of the function whenever the first term is fixed, and we will present algorithms to perform these computations.

The first equation in (16) amounts to choosing the origin of coordinates in the parameterization of the torus and, therefore eliminates the ambiguity corresponding to . (Check how does (16) change when we choose ).

The second equation in (16) indicates that is chosen to be a vector in the hyperbolic direction. We furthermore require that we have chosen the coordinate so that it is an eigenvector of the expanding/contracting direction.

The third equation in (16) chooses the eigenvalue. Equivalently, it fixes the scale in the variables . Observe that, setting amounts to multiplying by . Hence, setting the norm of sets the .

From the mathematical point of view, all choices of are equivalent. Nevertheless, from the numerical point of view, it is highly advantageous to choose so that the numerical coefficients of the expansion (in ) of have norms that neither grow nor decrease fast. This makes the computation more immune to round off error since it becomes more important when we add numbers of very different sizes.

### 3.6. Fourier-Taylor discretization

One of the ingredients of algorithms to solve the functional equations is to consider discretizations of functions one searches for.

In this section, we introduce the discretizations we will use. Roughly, for periodic functions, we will use both a Fourier series discretization and a real discretization on a grid. We will show that the Newton step can be decomposed into substeps which require only operations in either of the representations. Of course, one can switch between both representations using operations using FFT algorithms. For the study of invariant manifolds, we will use Taylor series in the real variables.

#### 3.6.1. Fourier series discretization

Since we are seeking functions which are periodic in the angle variable , it is natural to discretize them retaining a finite number of their Fourier coefficients

 K(θ)=∑k∈Zℓ,k∈ONcke2iπk⋅θ, (17)

where

 ON={k∈Zℓ||k|≤N}.

Since we will deal with real-valued functions, we have and one can just consider the cosine and sine Fourier series,

 K(θ)=a0+∑k∈ONakcos(2πk⋅θ)+bksin(2πk⋅θ). (18)

These Fourier discretizations have a very long history going back to classical astronomy, but have become much more widely used with computers and go under different name such as “automatic differentiation”. The manipulation of these polynomials are reviewed in [Knu97]. A recent review of their applications in dynamics – including implementation issues and examples – is [Har08].

The main shortcoming of Fourier series discretization of a function is that they are not adaptative and that for discontinuous functions, they converge very slowly and not uniformly. These shortcomings are however not very serious for our applications.

Since the tori are invariant under rigid rotations, they tend to be very homogeneous, so that adaptativity is not a great advantage. Also, it is known [FdlLS09a] that if tori are for sufficiently large , they are in fact analytic.

The fact that the Fourier series converge slowly for functions with discontinuities is a slight problem if one wants to compute tori close to the breakdown of analyticity, when the tori transform into Aubry-Mather objects. Of course, when they are far from breakdown – as it happens in many interesting problems in celestial mechanics – the Fourier coeffients converge very fast. To perform calculations close to breakdown, the a posteriori theorems in [FdlLS09a] prove invaluable help to have confidence in the computed objects.

#### 3.6.2. Fourier vs grid representation

Another representation of the function is to store the values in a regularly space grid. For functions of variables, we see that if we want to use variables, we can store either the Fourier coefficients of index up to or the values on a grid of step .

Some operations are very fast on the real space variables, for example multiplication of functions (it suffices to multiply values at the points of the grid). Also, the evaluation of is very fast if we discretize on a grid (we just need to evaluate the function for each of the points on the grid). Other operations are fast in Fourier representation. For example, it is fast to shift the functions, to take derivatives and, as we will see in (3.7), to solve cohomology equations. Hence, our iterative step will consist in the application of several operations, all of which being fast – – either in Fourier mode representation or in a grid representation. Of course, using the Fast Fourier Transform, we can pass from a grid representation to Fourier coefficients or viceversa in operations. There are extremely efficient implementations of the FFT algorithm that take into account not only operation counts but also several other characteristics (memory access, cache, etc.) of modern computers.

#### 3.6.3. Fourier-Taylor series

For the computation of whiskers of invariant tori, we will use Fourier-Taylor expansions of the form

 W(θ,s)=∞∑n=0Wn(θ)sn, (19)

where are -periodic functions in which we will approximate using Fourier series (17).

To manipulate this type of series we will use the so called automatic differentiation algorithms (see [Knu97],[Har08]). For the basic algebraic operations and the elementary transcendental functions (exp, sin, cos, log, power, etc.), they provide an expression for the Taylor coefficients of the result in terms of the coefficients of each of the terms.

### 3.7. Cohomology equations and Fourier discretization

In the Newton step to construct KAM tori, one faces solving cohomology equations, that is, given a periodic (on ) function , we want to find another periodic function solving (the first equation is a small divisor equation for flows and the second one for maps)

 ∂ωφ=η,\@@LTX@noalign\vskip6.0pt\omitφ−φ∘Tω=η.

As it is well known, equations (3.7) have a solution provided that

 ^η0≡∫Tℓη=0, (20)

and that is Diophantine in the appropriate sense. The Fourier coefficients of the solution of (3.7) are then given respectevely by

 ^φk=^ηk2πiω⋅k,\@@LTX@noalign\vskip6.0pt\omit^φk=^ηk1−e2πik⋅ω,

where are the Fourier coefficients of the function .

Notice that the solution is unique up to the addition of a constant (the average of is arbitrary).

Equations (3.7) and their solutions (3.7) are very standard in KAM theory (see the exposition in [Lla01]). Very detailed estimates can be found in [Rüs75], when is Diophantine (which is our case).

## 4. Fast Newton methods for (possibly) whiskered tori

In this section we develop an efficient Newton method to solve the invariance equations (8)-(11) and (10)-(13). We mainly focus on the case of maps (the case for vector fields being similar is described in the appendices).

We emphasize that the algorithm applies both to whiskered tori and to Lagrangian tori. Indeed, the case of Lagrangian tori is simpler. The hyperbolic part of the Lagrangian tori is just empty so that we do not need to compute the splittings. We refer to Algorithm 1 and Remark 6.

We will assume that the motion on the torus is a rigid rotation with a Diophantine frequency . As we have already shown, the invariance implies that the vectors in the range of are invariant under . The preservation of the symplectic structure implies that the vectors in the range of grow at most polynomially under iteration. We note also that tori with an irrational rotation are co-isotropic, , i.e.

 RangeDK∩Range(J∘K)−1DK={0}, (21)

and therefore . Hence, at any point of the invariant torus of dimension with motion conjugate to a rotation, we can find a -dimensional space of vectors that grow at most polynomially under iteration. As it is shown in [FdlLS09b], approximately invariant tori are approximately co-isotropic and the transversality (21) also holds.

The tori that we will consider are as hyperbolic as possible, given the previous argument. That is, we will assume that there exist directions that contract exponentially in the past or in the future, which span the complement of the tangent to the torus and its symplectic conjugate.

We will consider tori that have a hyperbolic splitting

 TK(θ)M=EcK(θ)⊕EsK(θ)⊕EuK(θ), (22)

such that there exist , satisfying , and such that for all and

 \@@LTX@noalign\vskip6.0pt\omitv∈EsK(θ)⟺|M(n,θ)v|≤Cμn1|v|∀n≥1\@@LTX@noalign\vskip6.0pt\omitv∈EuK(θ)⟺|M(n,θ)v|≤Cμn2|v|∀n≤1\@@LTX@noalign\vskip6.0pt\omitv∈EcK(θ)⟺|M(n,θ)v|≤Cμn3|v|∀n∈Z

where is the cocycle with generator and frequency , i.e. is given by

 M(n,θ)=⎧⎨⎩Z(θ+(n−1)ω)⋯Z(θ)n≥1,Idn=0,Z−1(θ+(n+1)ω)⋯Z−1(θ)n≤1. (23)

We will also assume that

 dimEcK(θ)=2ℓ,dimEsK(θ)=dimEuK(θ)=d−ℓ. (24)

The assumption (24) implies that the only non-hyperbolic directions are those spanned by the tangent to the torus and its symplectic conjugate, that is, there are no elliptic directions except those that are forced by the symplectic structure and the fact that the motion on the torus is a rotation.

We associate to the splitting (22) the projections , and over the invariant spaces , and .

It is important to note that since is symplectic (i.e. ), for all and

 Ω(u,v)=Ω(DFnu,DFnv),

so that, if have rates of decrease, by taking limits in the appropriate direction we obtain that is zero. That is, we get

 Ω(Es,Es)=0, Ω(Eu,Eu)=0, Ω(Ec,Es)=0, Ω(Ec,Eu)=0.

Therefore, we have

 (J(K(θ)))−1EcK(θ) = EcK(θ), (J(K(θ)))−1EsK(θ) = EuK(θ), (J(K(θ)))−1EuK(θ) = EsK(θ).

In [HdlLS09], we provide a method to compute the rank-1 bundles by iterating the cocycle. Of course, once we have computed the vector spanning the rank-1 (un)stable bundle it is very easy to obtain the projections. In Section 4.3 we discuss an alternative to compute the projections by means of a Newton method. In this case we do not need to assume that the bundle is 1-dimensional.

### 4.1. General strategy of the Newton method to solve the invariance equation

In this section we will design a Newton method to solve the invariance equation (8) and the modified one (10), and discuss several algorithms to deal with the linearized equations.

We first define the following concept of approximate solution.

###### Definition 1.

We say that (resp. ) is an approximate solution of equation (8) (resp. (10)) if

 F∘K−K∘Tω=E,\@@LTX@noalign\vskip6.0pt\omit(resp.F∘K−K∘Tω−((J∘K0)−1DK0)∘Tωλ=E),

where is small.

The Newton method consists in computing in such a way that setting and expanding the LHS of (1) in up to order , it cancels the error term .

###### Remark 5.

Throughout the paper, we are going to denote some norms in functional spaces without specifying what they are exactly. We refer the reader to [LGJV05, FdlLS09a], where the whole theory is developed and the convergence of the algorithms is proved. Recall that one of the key ideas of KAM theory is that the norms are modified at each step.

Performing a straightforward calculation, we obtain that the Newton procedure to solve equation (8) and (11), given an approximate solution , consists in finding satisfying

 (DF∘K)Δ−Δ∘Tω=−E. (25)

For the modified invariance equation (10), given an approximate solution , the Newton method consists in looking for in such a way that and eliminate the error in first order. The linearized equation in this case is

 (DF∘K)Δ−Δ∘Tω−((J∘K0)−1DK0)∘Tωδ=−E, (26)

where one can take .

As it is well known, the Newton method converges quadratically in and the error at step is such that

 ∥˜E∥≤C∥E∥2

where is the error at the previous step.

In order to solve the linearized equations (25) and (26), we will first project them on the invariant subspaces and , and then solve an equation for each subspace.

Thus, let us denote

 \@@LTX@noalign\vskip6.0pt\omitΔs,c,u(θ)=Πs,c,uK(θ)Δ(θ),\@@LTX@noalign\vskip6.0pt\omitEs,c,u(θ)=Πs,c,uK(θ+ω)E(θ),

such that . Then, by the invariant properties of the splitting, the linearized equations for the Newton method (25) and (26) split into:

 DF(K(θ))Δc(θ)−Δc∘Tω(θ)=−Ec(θ),\@@LTX@noalign\vskip6.0pt\omitDF(K(θ))Δs(θ)−Δs∘Tω(θ)=−Es(θ),\@@LTX@noalign\vskip6.0pt\omitDF(K(θ))Δu(θ)−Δu∘Tω(θ)=−Eu(θ),

and

 DF(K(θ))Δc(θ)−Δc∘Tω(θ)+ΠcK(θ+ω)(J∘K0(θ+ω))−1DK0(θ+ω)δ=−Ec(θ),\@@LTX@noalign\vskip6.0pt\omitDF(K(θ))Δs(θ)−Δs∘Tω(θ)+ΠsK(θ+ω)(J∘K0(θ+ω))−1DK0(θ+ω)δ=−Es(θ),\@@LTX@noalign\vskip6.0pt\omitDF(K(θ))Δu(θ)−Δu∘Tω(θ)+ΠuK(θ+ω)(J∘K0(θ+ω))−1DK0(θ+ω)δ=−Eu(θ).

Notice that once is obtained, the equations (4.1) on the hyperbolic spaces reduce to equations of the form (4.1). More precisely,

 DF(K(θ))Δs,u(θ)−Δs,u∘Tω(θ)=−˜Es,u(θ) (27)

where

 ˜Es,u=Es,u(θ)+Πs,uK(θ+ω)(J∘K0(θ+ω))−1DK0(θ+ω)δ.

Equations (4.1) and (4.1) for the stable and unstable spaces can be solved iteratively using the contraction properties of the cocycles on the hyperbolic spaces given in (4). Indeed, a solution for equations (27) is given explicitly by

 Δs(θ)=˜Es∘T−ω(θ)+∞∑k=1(DF∘K∘T−ω(θ)×⋯×DF∘K∘T−kω(θ))(˜Es∘T−(k+1)ω(θ)) (28)

for the stable equation, and

 Δu(θ)=−∞∑k=0(DF−1∘K(θ)×⋯×DF−1∘K∘Tkω(θ))(˜Eu∘Tkω(θ)) (29)

for the unstable direction. Of course, the contraction of the cocycles guarantees the uniform convergence of these series.

The algorithms presented in Appendix A allow us to compute the solutions of equations (27) efficiently.

In Section 4.2 we discuss how to solve equations (4.1) and (4.1) in the center direction.

Hence, the Newton step of the algorithm for whiskered tori that we summarize here will be obtained by combining several algorithms.

###### Algorithm 1.

Consider given , , and an approximate solution (resp. ), perform the following operations:

1. Compute the invariant splittings and the projections associated to the cocycle and using the algorithms described in Section 4.3 (or in [HdlLS09]).

2. Project the linearized equation to the hyperbolic space and use the algorithms described in Appendix A to obtain .

3. Project the linearized equation on the center subspace and use the Algorithm 3 in Section 4.2 to obtain and .

4. Set and

Of course, since this is a Newton step, it will have to be iterated repeatedly until one reaches solutions up to a small tolerance error.

We will start by some remarks on the different steps of Algorithm 1 and, later, we will provide many more details on them.

###### Remark 6.

It is important to remark that the above Algorithm 1 also applies to the case of Lagrangian tori. It suffices to remark that in that, case, the center space is the whole manifold, so that there is no need to compute the splitting. Hence, for Lagrangian tori, the steps A) and B) of Algorithm 1 are trivial and do not need any work.

###### Remark 7.

The main issue of the Newton method is that it needs a good initial guess to start the iteration. Any reasonable algorithm can be used as an input to the Newton method. Indeed, our problems have enough structure so that one can use Lindstedt series, variational methods, approximation by periodic orbits, frequency methods, besides the customary continuation methods.

###### Remark 8.

As we have mentioned in Remark 3, the solutions of (8) and (11) are not unique. Therefore, in order to avoid dealing with non-invertible matrices in the Newton procedure, we will impose the normalization condition

 ∫Tℓ(K(θ)−I(θ))⋅νi=0

where is a basis for ( being the dimension) and is a linear function introduced in (14).

### 4.2. Fast Newton method for (whiskered) tori: the center directions

We present here the Newton method to solve the equations on the center subspace in the case of maps.

For Lagrangian tori, the hyperbolic directions are empty and the study of the center direction is the only component which is needed. Hence, the algorithms discussed in this section allow to solve, in particular, equations (25) and (26) in the case of Lagrangian tori. For a discussion of the center equations for Hamiltonian flows, we refer the reader to Appendix B.

The key observation is that the linearized Newton equations (25) and (26) are closely related to the dynamics and therefore, we can use geometric identities to find a linear change of variables that reduces the Newton equations to upper diagonal difference equations with constant coefficients. This phenomenon is often called “automatic reducibility”.

The idea is stated in the following proposition:

###### Proposition 2 (Automatic reducibility, see [FdlLS09b, FdlLS09a]).

Given an approximation of the invariance equation as in (1), denote

 α(θ)=DK(θ)\@@LTX@noalign\vskip6.0pt\omitN(θ)=([α(θ)]⊥α(θ))−1\@@LTX@noalign\vskip6.0pt\omitβ(θ)=α(θ)N(θ)\@@LTX@noalign\vskip6.0pt\omitγ(θ)=(J∘K(θ))−1β(θ)

and form the following matrix

 M(θ)=[α(θ)∣γ(θ)], (30)

where by we denote the matrix obtained by juxtaposing the two matrices that are in the arguments.

Then, we have

 (31)

where

 A(θ)=β(θ+ω)⊥[(DF∘K(θ))γ(θ)−γ(θ+ω)], (32)

and in the case of (25) or in the case of (26).

###### Remark 9.

If the symplectic structure is almost-complex (i.e. ), we have that

 β(θ+ω)⊥γ(θ+ω)=0,

since the torus is isotropic. Then has a simpler expression given by

 A(θ)=β(θ+ω)⊥(DF∘K)(θ)γ(θ).

Once again, we omit the definition of the norms used in the bounds for . For these precisions, we refer to the paper [FdlLS09a], where the convergence of the algorithm is established.

It is interesting to pay attention to the geometric interpretation of the identity (31). Note that, taking derivatives with respect to in (1), we obtain that

 (DF∘K)DK−DK∘Tω=DE,

which means that the vectors are invariant under (up to a certain error). Moreover, are the symplectic conjugate vectors of , so that the preservation of the symplectic form clearly implies (31). The geometric interpretation of the matrix is a shear flow near the approximately invariant torus. See Figure 1.

To be able to use the change of unknowns via the matrix previously introduced on the center subspace, one has to ensure that one can identify the center space with the range of . This is proved in [FdlLS09a] to which we refer.

For our purposes it is important to compute not just the invariant spaces, but also the projections over invariant subspaces. Knowing one invariant subspace is not enough to compute the projection, since it also depends on the complementary space chosen.

In the following, we will see that the result stated in Proposition 2 allows us to design a very efficient algorithm for the Newton step.

Notice first that if we change the unknowns in (25) and (26) and we use (31) we obtain

 M(θ+ω)(IdA(θ)0Id)W(θ)−M(θ+ω)W(θ+ω)\@@LTX@noalign\vskip6.0pt\omit−(J(K0(θ+ω)))−1DK0(θ+ω)δ=−E(θ)

Of course, the term involving has to be omitted when considering (25).

Multiplying (4.2) by and using the invertibility of the matrix (see [FdlLS09b, FdlLS09a]), we are left with the system of equations

 W1(θ)+A(θ)W2(θ)−B1(θ)δ−W1(θ+ω)=−˜E1(θ)\@@LTX@noalign\vskip6.0pt\omitW2(θ)−W2(θ+ω)−B2(θ)δ=−˜E2(θ)

where

 ˜E(θ)=(M(θ+ω)⊥J(K(θ+ω))M(θ+ω))−1M(θ+ω)⊥J(K(θ+ω))E(θ)\@@LTX@noalign\vskip6.0pt\omitB(θ)={(M⊥J(K)M)−1M⊥J(K)(J(K0))−1DK0}∘Tω(θ)

and the subindices indicate symplectic coordinates.

When is close to , we expect that is close to the -dimensional identity matrix and is small.

The next step is to solve equations (4.2) for (and ). Equations (4.2) are equations of the form considered in (3.7) and they can be solved very efficiently in Fourier space.

More precisely, the second equation of (4.2) is uncoupled from the first one and allows us to determine (up to a constant) and . The role of the parameter is now clear. It allows us to adjust some global averages that we need to be able to solve equations (4.2). Indeed, we choose so that the term has zero average (which is a necessary condition to solve small divisor equations as described in Section 3.7). This allows us to solve equation (3.7) for . We then denote

 W2(θ)=˜W2(θ)+¯¯¯¯¯¯W2

where has average zero and .

Once we have , we can substitute in the first equation. We get imposing that the average of

 B1(θ)δ−A(θ)˜W2(θ)−A(θ)¯¯¯¯¯¯W2(θ)−˜E1(θ)

is zero and then we can find up to a constant according to (3.7).

We therefore have the following algorithm to solve (3) in the center direction,

###### Algorithm 3 ( Newton step in the center direction).

Consider given , , and an approximate solution (resp. ). Perform the following calculations

• Compute

• Compute

• Compute the invariant projections, .

• Set (resp. set