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

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 are based on a Newton method to solve an appropriately chosen functional equation that expresses invariance. The algorithms are efficient: if we discretize the objects by $N$ elements, one step of the Newton method requires only O(N) storage and $O(N \ln(N))$ operations. Furthermore, if the object we consider is of dimension $\ell$, we only need to compute functions of $\ell$ 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 \emph{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 have already been implemented. We will report on the technicalities of the implementation and the results of running them elsewhere.


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.
• 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 quasiperiodic 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 N Fourier coefficients, we require O(N ) storage and only O(N ln(N )) operations for the Newton step. Since the functions we are considering are analytic, we see that the truncation error is O(exp(−CN 1/d )) where d is the dimension of the object. Note that, in contrast, a straightforward implementation of a Newton method would require to use O(N 2 ) storageto store the linearization matrix and its inverse -and O(N 3 ) 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 N log(N ) 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 F is a map and we can find an embedding K 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.
The main idea of the Newton method is that, using the decomposition into invariant subspaces, one can decompose (2) into three components where the s, u refer to the stable, unstable components and c refers to the component along the tangent to the torus and its symplectic conjugate. For Lagrangian tori, only the E c 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 K(θ + ω)) and a real representation which makes easy to evaluate F (K(θ)). 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-1 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 η = 0 (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 F is exact symplectic, the obstructions for the solution are O(||E|| 2 ). An alternative method to deal with these obstructions is to add some new -finite dimensional -unknowns λ and solve, instead of (1), the equation where G(θ) 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 λ = 0. 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 K is a solution,K defined byK(θ) = K(θ + σ) is also a solution for any σ ∈ R ℓ . 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−1 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 W 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 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 1 set (a Cantor set of codimension 1 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.

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 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 actionangle 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 M is endowed with an exact symplectic structure Ω = dα (for some one-form α) and we have where ·, · denotes the inner product on the tangent space of M and J(z) is a skew-symmetric matrix.
An important particular case is when J induces an almost-complex structure, i.e.
Most of our calculations do not need this assumption. One important case, where the identity (6) is not satisfied, is when J is a symplectic structure on surfaces of section chosen arbitrarily in the energy surface or when J 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 F : U ⊂ M → M which are not only symplectic (i.e. F * Ω = Ω) but exact symplectic, that is F * α = α + dP, for some smooth function P , 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: 2.2. Systems described by vector fields. We will assume that the system is described by a globally Hamiltonian vector-field X, that is where H is a globally defined function on T * M.
In the case of flows, the appropriate notion of Diophantine numbers is: 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.

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 quasiperiodic solutions with frequency ω (independent over the integers) under the form of Fourier series where ω ∈ R ℓ and n ∈ Z.
We allow some components of x in (7) to be angles. In that case, it suffices to take some of the components of x modulo 1.
It is then natural to describe a quasi-periodic function using the so-called "hull" function K : T ℓ → M defined by so that we can write x (n) = K(nω). The geometric interpretation of the hull function is that it gives an embedding from T ℓ 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 F if and only if the hull function K satisfies: where T ω denotes a rigid rotation A modification of the invariance equations (8) which we will be important for our purpose consists in considering where the unknowns are now K : T ℓ → M (as before) and λ ∈ R ℓ . Here, K 0 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 (K, λ) satisfy the equation (10) with K 0 close to K, then at the end of the iteration of the Newton method, we have λ = 0 and, therefore, K 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 (24) 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 nondegeneracy 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 2d. Hence, the algorithms presented here have the advantage that they look for a function K 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 h into R 2d , one needs to store (1/h) ℓ · 2d real values.
Remark 3. Equations (8) and (10) do not have unique solutions. Observe that if K is a solution, for any σ ∈ R ℓ , K • T σ 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 K. Hence, it is easy to get rid of it by imposing some normalization. See Section 3.5.2 for a discussion on this issue.

Functional equations for whiskered invariant tori for vectorfields.
In this case, one can write where ω ∈ R ℓ , t ∈ R and then the hull function K is defined by The invariance equation for flows is written: where ∂ ω denotes the derivative in direction ω The modification of (11) incorporating a counterterm is: where K 0 is a given embedding satisfying some non-degeneracy conditions.
Remark 4. Recall that, taking time−1 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 2 the dimension of the phase space and the parameterization of the torus requires 1 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 T ℓ and the range of K have topology. As a consequence, there will be some topological considerations in the way that the torus T ℓ gets embedded in the phase space. More explicitly, the angle variables of T ℓ 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 K to the universal cover, i.e.
in such a way that K is obtained from K by identifying variables in the domain and in the range that differ by an integer.
It is therefore clear that ∀ e ∈ Z ℓ where K p , K q denote the projections of the lift on the p and q coordinates of R 2d−ℓ × R ℓ . It is easy to see that I(e) is a linear function of e, namely with I ij ∈ Z. We note that if a function K q satisfies The numerical methods will always be based on studying the periodic functions K q , but we will not emphasize this unless it can lead to confusion.
Of course, the integer valued matrix I = {I ij } ij remains constant if we modify the embedding slightly. Hence, it remains constant under continuous deformation. For example, in the integrable case with ℓ = d, invariant tori satisfy K q (θ) = θ, so that we have I = Id. Hence, all the invariant tori which can be continued from tori of the integrable system will also have I = Id.
3.4. Secondary tori. One can produce other ℓ-dimensional tori for which the range of I is of dimension less than ℓ. These tori are known as secondary tori. It is easy to see that if rank(I) < ℓ we can contract K(T ℓ ) to a diffeomorphic copy of T rank(I) . Even in the case of maximal tori ℓ = d, 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 standardlike 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 K can be adjusted by the matrix I introduced in (14). Hence, the rank of the matrix I has to be chosen according to the number of contractible directions.
3.5. Equations for the invariant whiskers. Invariant tori with ℓ < d 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 nonresonance 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 onedimensional 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 u of the equations of motion under the form in the discrete time case and The function W has then to satisfy the following invariance equations for the case of maps and flows, respectively. Note that equations (17) imply that in variables (θ, s) 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 µ (e µ in the case of flows). We call contractive the situation |µ| < 1 for maps (or µ < 0 for flows). We call expansive the case when |µ| > 1 for maps (or µ > 0 for flows). Note that if W (θ, s) satisfies (17) then W (θ, 0) 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 where K 0 is, as before, a given approximate solution of the equations (8) and (11) 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 where D 2 W denotes the derivative with respect to the second argument, ρ > 0 is any arbitratrily chosen number and . stands for a suitable norm.
The fact that the solutions of (8) supplemented by (19) are locally unique is proved in [FdlLS09a]. In this paper, we will see that these normalizations uniquely determine the Taylor expansions (in s) of the function W whenever the first term W 1 (θ) ≡ D 2 W (θ, 0) is fixed, and we will present algorithms to perform these computations.
The first equation in (19) amounts to choosing the origin of coordinates in the parameterization of the torus and, therefore eliminates the ambiguity corresponding to σ. (Check how does (19) change when we choose σ).
The second equation in (19) indicates that W 1 (θ) 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 (19) chooses the eigenvalue. Equivalently, it fixes the scale in the variables s. Observe that, setting b amounts to multiplying W 1 by b. Hence, setting the norm of DW sets the b.
From the mathematical point of view, all choices of ρ are equivalent. Nevertheless, from the numerical point of view, it is highly advantageous to choose ||W 1 || so that the numerical coefficients of the expansion (in s) of W 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 O(N ) operations in either of the representations. Of course, one can switch between both representations using O(N ln(N )) 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 K which are periodic in the angle variable θ, it is natural to discretize them retaining a finite number of their Fourier coefficients where Since we will deal with real-valued functions, we have c k =c −k and one can just consider the cosine and sine Fourier series, 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 examplesis [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 C r for sufficiently large r, 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 K is to store the values in a regularly space grid. For functions of ℓ variables, we see that if we want to use N variables, we can store either the Fourier coefficients of index up to O(N 1/ℓ ) or the values on a grid of step O(N −1/ℓ ).
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 F • K is very fast if we discretize on a grid (we just need to evaluate the function F 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 -O(N ) -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 O(N ln N ) 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 where W n are 1-periodic functions in θ which we will approximate using Fourier series (20).
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 T ℓ ) 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) As it is well known, equations (23) have a solution provided that and that ω is Diophantine in the appropriate sense. The Fourier coefficientŝ ϕ k of the solution ϕ of (25) are then given respectevely bŷ whereη k are the Fourier coefficients of the function η.
Notice that the solution ϕ is unique up to the addition of a constant (the averageφ 0 of ϕ is arbitrary).
Equations (23) and their solutions (25) 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).

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 ω ∈ R ℓ . As we have already shown, the invariance implies that the vectors in the range of DK are invariant under DF . The preservation of the symplectic structure implies that the vectors in the range of (J • K) −1 DK grow at most polynomially under iteration. We note also that tori with an irrational rotation are co-isotropic, and therefore dim(Range DK ⊕ Range (J • K) −1 DK) = 2ℓ. Hence, at any point of the invariant torus of dimension ℓ with motion conjugate to a rotation, we can find a 2ℓ-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 (26) 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 such that there exist 0 < µ 1 , µ 2 < 1, µ 3 > 1 satisfying µ 1 µ 3 < 1, µ 2 µ 3 < 1 and C > 0 such that for all n ≥ 1 and where M(n, θ) is the cocycle with generator Z(θ) = DF (K(θ)) and frequency ω, i.e. M : Z × T ℓ → GL(2d, R) is given by We will also assume that The assumption (30) 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 (27) the projections Π c K(θ) , Π s K(θ) and Π u K(θ) over the invariant spaces E c K(θ) , E s K(θ) and E u K(θ) . It is important to note that since F is symplectic (i.e. F * Ω = Ω), for all n ≥ 1 and n ≤ −1 Ω(u, v) = Ω(DF n u, DF n v), so that, if u, v have rates of decrease, by taking limits in the appropriate direction we obtain that Ω is zero. That is, we get Therefore, we have . 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 K (resp. (K, λ)) is an approximate solution of equation (8) where E is small.
The Newton method consists in computing ∆ in such a way that setting K ← K + ∆ and expanding the LHS of (31) in ∆ up to order ∆ 2 , it cancels the error term E.
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 K, consists in finding ∆ satisfying For the modified invariance equation (10), given an approximate solution (K, λ), the Newton method consists in looking for (∆, δ) in such a way that K + ∆ and λ + δ eliminate the error in first order. The linearized equation in this case is where one can take K 0 = K.
As it is well known, the Newton method converges quadratically in E and the error E at step K + ∆ is such that where E is the error at the previous step.
In order to solve the linearized equations (32) and (33), we will first project them on the invariant subspaces E c , E u and E s , and then solve an equation for each subspace.
The algorithms presented in Appendix A allow us to compute the solutions ∆ s,u of equations (37) efficiently.
In Section 4.2 we discuss how to solve equations (35) and (36) 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 F , ω, K 0 and an approximate solution K (resp. K, λ), perform the following operations: A) Compute the invariant splittings and the projections associated to the cocycle Z(θ) = DF • K(θ) and ω using the algorithms described in Section 4.3 (or in [HdlLS09]). B) Project the linearized equation to the hyperbolic space and use the algorithms described in Appendix A to obtain ∆ s,u . C) Project the linearized equation on the center subspace and use the Algorithm 3 in Section 4.2 to obtain ∆ c 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 is a basis for Range(I) (L being the dimension) and I is a linear function introduced in (15).

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 (32) and (33) 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 (32) and (33) 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 K of the invariance equation as in (31), denote and form the following matrix where by [· | ·] we denote the 2d × 2ℓ matrix obtained by juxtaposing the two 2d × ℓ matrices that are in the arguments.
Then, we have where and E ≤ DE in the case of (32) or E ≤ DE + |λ| in the case of (33).
Remark 9. If the symplectic structure is almost-complex (i.e. J 2 = − Id), we have that β(θ + ω) ⊥ γ(θ + ω) = 0, since the torus is isotropic. Then A(θ) has a simpler expression given by Once again, we omit the definition of the norms used in the bounds for E. 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 (42). Note that, taking derivatives with respect to θ in (31), we obtain that  To be able to use the change of unknowns via the matrix M previously introduced on the center subspace, one has to ensure that one can identify the center space E c K(θ) with the range of M . 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.
When K is close to K 0 , we expect that B 2 is close to the ℓ-dimensional identity matrix and B 1 is small. The next step is to solve equations (45) for W (and δ). Equations (45) are equations of the form considered in (23) and they can be solved very efficiently in Fourier space.
More precisely, the second equation of (45) is uncoupled from the first one and allows us to determine W 2 (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 (45). Indeed, we choose δ so that the term B 2 (θ)δ − E 2 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 (23) for W 2 . We then denote W 2 (θ) = W 2 (θ) + W 2 where W 2 (θ) has average zero and W 2 ∈ R.
Once we have W 2 , we can substitute W 2 in the first equation. We get W 2 imposing that the average of is zero and then we can find W 1 up to a constant according to (25).
We therefore have the following algorithm to solve (3) in the center direction, Algorithm 3 ( Newton step in the center direction). Consider given F , ω, K 0 and an approximate solution K (resp. K, λ). Perform the following calculations We denote E 1 , E 2 the components of E along DK and along J −1 DK. (3.9) Compute as indicated in (43) 4. (4.1) Solve for W 2 satisfying Set W 2 such that the average is 0.) 5. (5.1) Compute A(θ)W 2 (θ) (5.2) Solve for W 2 satisfying Notice that steps (1.2), (3.1), (3.6), (4.1) (resp. (4.2 ′ )), (5.3) (resp. (5.3 ′ )) in Algorithm 3 are diagonal in Fourier series, whereas the other steps are diagonal in the real space representation. Note also that the algorithm only stores vectors whose size is of order N .
Remark 10. Using the symplectic properties of the matrix M , step (3.8) can be sped up.
When the torus is exactly invariant we have that the invariant torus is co-isotropic. Hence DK ⊥ J • KDK = 0. Hence, when the torus is invariant, we have M (θ + ω) ⊥ J(K(θ + ω))M (θ + ω) = 0 N −N 0 so that the inverse is easy to calculate.
For the purposes of a Newton Method, we can use the same expression for the inverse in step (3.8) and still obtain a quadratically convergent algorithm.

4.3.
A Newton method to compute the projections over invariant subspaces. In this section we will discuss a Newton method to compute the projections Π c K(θ) , Π s K(θ) and Π u K(θ) associated to the linear spaces E c K(θ) , E s K(θ) and E u K(θ) where K is an (approximate) invariant torus. More precisely, we will design a Newton method to compute Π s K(θ) and Π cu K(θ) = Π c K(θ) + Π u K(θ) . Similar arguments allow us to design a Newton method to compute Π u K(θ) and Π cs K(θ) = Π c K(θ) + Π s K(θ) . Then, of course, Π c K(θ) is given by Let us discuss first a Newton method to compute Π s K(θ) and Π cu K(θ) . To simplify notation, from now on, we will omit the dependence in K(θ).
We are going to design a Newton method to solve equations (46)-(47) and use equations (48)-(49) as constraints. In this context, by approximate solution of equations (46)-(47), we mean a solution (Π s , Π cu ) such that where E i denotes the error in a certain component. Notice that the error in equation (53) has components only on the center and unstable "approximated" subspaces and we denote it by E cu . The same happens with the equation (54) but on the "approximated" stable subspace. We assume that E cu and E s are both small. As standard in the Newton method, we will look for increments ∆ s and ∆ cu in such a way that setting Π s ← Π s + ∆ s and Π cu ← Π cu + ∆ cu , the new projections solve equations (46) and (47) up to order E 2 where E = E s + E cu for some norm . .
The explicit expressions for ∆ s cu and ∆ s s are and Remark 11. Notice that by the way N cu (θ) is defined, it is a matrix which does not have full rank. Therefore, we understand N −1 cu (θ) as the "pseudoinverse" matrix.
Finally, let us check that ∆ s = ∆ s cu + ∆ s s also satisfies the constraints. In order to check that constraint (59), which is equivalent to (63), is satisfied we will use the expressions (69) and (70). Notice first that and Moreover, from (53) and using (56) one can see that Then, from expressions (69) and (70) and the above expression (71), (72) and (73), it is clear that Now, using equation (58) we get ∆ cu (θ) = −(∆ s s (θ) + ∆ s cu (θ)) and the improved projections are The new error for equations (46) and (47) is now E ≤ C E 2 where E = E cu + E s . Of course equation (48) is clearly satisfied but (49) is satisfied up to an error which is quadratic in E . However it is easy to get an exact solution for (49) and the correction is quadratic in ∆ s (and therefore in ∆ cu ). To do so, we just take the SVD decomposition of Π s and we set the values in the singular value decomposition to be either 1 or 0.
In this way we obtain new projections Π s new and Π cu new = Id −Π s new satisfying new − Π cu < ∆ cu 2 , so that the error for equations (46) and (47) is still quadratic in E . Moreover, they satisfy equations (49) and, of course, (48) exactly.
Hence, setting Π s ← Π s new and Π cu ← Π cu new we can repeat the procedure described in this section and perform another Newton step.
Notice that the matrix multiplication is diagonal in real space representation, whereas the phase shift is diagonal in Fourier space. A discussion on how to perform step 4 efficiently is given in Appendix A.

Computation of rank-1 whiskers of an invariant torus
In this section, we present algorithms to compute the whiskers associated to an invariant torus, that is the invariant manifolds that contain the torus and are tangent to the invariant bundles.
For the sake of simplicity and in order to state in a clear way the main idea behind the methods we will only discuss the case when the invariant whiskers are one-dimensional (i.e. d−ℓ = 1). The same idea can be extended to compute invariant manifolds of any rank. However, there are several new phenomena (resonances) that can appear and need to be discussed. We plan to come back to this issue in the future.
As we already mentioned in Section 3.5.1, we will look for the whiskers by finding a parameterization for them, so we will search for a function W : T ℓ × (V ⊂ R d−ℓ ) → M and a scalar µ satisfying equation (17).
We will consider three different methods to solve equation (17). We will first discuss the order by order method. The other two methods are based on the philosophy of quasi-Newton methods. Using the phenomenon of "automatic reducibility", we are able to design an efficient Newton method. The first method allows to compute simultaneously the invariant tori and the whiskers, whereas the second one assumes that the invariant tori and the tangent bundles are already known.
We detail only the case of maps because the same ideas work for the case of vector fields and refer the reader to the appendices for the case of flows.
Similar algorithms were developed and implemented in [HL06b,HL07] for the slightly simpler case of quasi-periodic maps.
5.1. The order by order method. In this section we adapt the parameterization method introduced in [CFL05]. The convergence of the Fourier-Taylor series in this paper can be easily adapted to the present case. We focus on the case of maps and refer the reader to Appendix C for the case of flows.
We will find a solution (W, µ) of the invariance equation (17) discretizing it in Fourier-Taylor series. Hence, we will look for W as a power series and match similar coefficients in s n on both sides of equation (17). For n = 0, we obtain which is equation (8) for the invariant torus. Therefore, we have W 0 (θ) = K(θ), where K is a parametrization of the invariant torus. For n = 1, we obtain so that W 1 (θ) is an eigenfunction with eigenvalue µ of the operator M(1, θ) defined in equation (29). Equation (76) states that the bundle spanned by W 1 is invariant for the linearization of F , and the dynamics on it is reduced to a contraction/expansion by a constant µ. Therefore, one can compute W 1 and µ using the algorithms given in Section 4.3.
Remark 12. Notice that if W 1 (θ) is a solution of equation (76), then bW 1 (θ), for any b ∈ R, is also a solution. See Section 3.5.2 for a discussion on how to choose b.
For n ≥ 2, we obtain where R n is an explicit polynomial in W 0 , . . . , W n−1 whose coefficients are derivatives of F evaluated at W 0 . Equation (77) can be solved provided that µ is such that µ n is not in the spectrum of the operator M(1, θ). This condition is clearly satisfied in the case of (un)stable bundles which are one-dimensional but it can also be satisfied by other bundles.
Equation (77) can be solved using the large matrix method. It consists on considering a discretization of equation (77) using Fourier series and reduce the problem to solving a linear system in Fourier space, where the unknowns are the Fourier coefficients of the matrix W n .
There are also efficient algorithms which are variants of the methods devoted in the previous sections. The equation (77) is equivalent to which, for large enough n is a contraction, so that we can apply the fast methods of Section A.1. In particular Algorithm A.8. In the case that the stable and unstable directions are one dimensional -which is the one we discuss in this paper -this is enough (remember that we always have n ≥ 2. When the bundles are higher dimensional, we may need to find a splitting corresponding to the cocycle generated by Z(θ) = (DF • K(θ)) −1 µ n .
Remark 13. Notice that once W 0 (θ) and W 1 (θ) are fixed, the solution W n (θ) for n ≥ 2 of equation (77) is uniquely determined. It is then clear that any analytic solution is unique. The existence of analytic solutions is discussed in [CFL05].
Remark 14. Notice that the equations to compute the new term W n do not involve small divisors.
Remark 15. In this case, we have not considered the modified equation (18) with the counterterm, because for this method there are no obstructions to deal with as it is in the case of the Newton method, see Remark 16. Indeed, the vanishing Lemma in [FdlLS09a] guarantees that for exact symplectic maps λ = 0 in (18) once the Newton method has converged.

5.2.
A Newton method to compute simultaneously the invariant torus and the whiskers. In this Section we present an algorithm to solve equation (17) using a Newton method, instead of solving it step by step as we discussed in the previous section. As before, we only deal with the case of maps and refer the reader to Appendix D for the case of flows. We do not prove the convergence of the algorithm here for sake of length (and purpose of the paper) but it can be done using techniques similar to those developed in [FdlLS09b,FdlLS09a].
We start with an initial approximation (W, µ) (resp. (W, µ, λ)) satisfying and we look for an improved solution Remark 16. As in the Newton method for invariant tori, the role of the parameter λ is to adjust some averages to solve the equations for the case n = 1. More precisely, δ will be chosen in equation (93) so that equation (91) can be solved without any obstruction.
We will try to solve equation (79) by discretizing it in Fourier-Taylor series. Notice that equation (79) is not diagonal when discretized in Fourier-Taylor series because of the term DF • W . However, there is a way to make it diagonal using the geometric identities which are a direct generalization of those used for the automatic reducibility.
We first give the idea of the automatic reducibility when W is such that Taking derivatives with respect to θ and s, we see that From the above equations, we read that the quantity D θ W (θ, s) remains invariant under DF • W (θ, s), whereas the vector ∂ s W (θ, s) is multiplied by a factor µ.
where A(θ, s) and B(θ, s) are some matrices, which will be computed as before. See Proposition 5. Therefore, one can see that in the basis Following the proofs in [FdlLS09a] for instance, one can prove the following proposition, which generalizes Proposition 2: and form the following matrix where we denote by [· | · | · | ·] the 2d × 2d matrix obtained by juxtaposing the two 2d × ℓ matrices and the two 2d × (d − ℓ) matrices that are in the arguments. Then Similarly, equation (92) do not involve small divisors nor obstructions. However, equation (91) does have obstructions and small divisors. In order to overcome this problem, we denote by F and G the solutions of where D 1 4 (θ) = − E 1 4 (θ) + δ µ C 1 4 (θ). This gives V 1 4 (θ) = δF (θ) + G(θ). Taking the average of equation (91), we have that so we can solve for δ provided that H 0 3 − B 0 F = 0. The other orders do not have any problem. For s n , with n ≥ 2, we have and equations (94) can be solved uniquely for V n 1 , V n 2 , V n 3 and V n 4 , for n = 2, . . . , L, where L is the degree for the Taylor expansion. Hence, we have obtained δ, δ µ ∈ R and Remark 18. For L = 1, the algorithm allows us to compute simultaneously the invariant torus and the associated linear subspaces.
The algorithm for the simultaneous computation of the whiskers and the invariant torus is Algorithm 6 (Computation of tori and whiskers simultaneously). Consider given F , ω, K 0 and a fixed order L. Given an approximate solution (W, µ, µ), perform the following calculations (11.6) Solve for V n 3 satisfying µV n 3 − µ n V n 3 • T ω = − E n 3 + δH n−1 A Newton method to compute the whiskers. Assuming that we have computed an invariant torus K(θ) with the associated stable direction V s (θ) (resp. unstable direction V u (θ)) and the rate of contraction (resp. expansion) µ, we can use a Newton method to compute the whiskers. We concentrate on the case of maps, referring to Appendix E for the flows. We consider the invariance equation (17), and we assume that we have an initial approximation W for the whiskers, expressed as a power series W n (θ)s n and such that W 0 (θ) = K(θ) and W 1 (θ) = V s (θ) (the unstable case is analogous).
Remark 19. Again, as discussed in Remark 15, we do not need to consider the modified invariance equation (18) with the counterterm. The fact that K(θ) is a solution of equation (8) for exact symplectic maps, implies that λ = 0 in (18). This is guaranteed by the vanishing Lemma in [FdlLS09a].
Then, it is clear that the error E for the initial approximation W is such that E(θ, s) = n≥2 E n (θ)s n , since this is exact for the terms of order 0 and 1. For a given function G : where Using this notation, the linearized equation for the Newton method is Similarly as we did in the previous section, we can perform the change of coordinates given by the matrix M (θ, s) in (82) and reduce the problem to solving for V (θ, s) the following equation, which is diagonal in Fourier-Taylor series, with R(θ, s) given in (83), E(θ, s) = M (θ + ω, µs) −1 E(θ, s) and ∆ = M V . Notice that in this case, we do not have to solve the system of equations for order 0 and 1 and we can go straight to order n ≥ 2. We use subindices to denote coordinates and superindices to denote the order. Hence, for order n ≥ 2 we need to solve the system of equations The solution of (96) for n = 2, 3 provides an exact solution of the invariance equation up to order 4. That is, if we set where V 2 and V 3 are obtained by solving equations (96), then the improved solutionW is given bȳ where M (θ, s) was introduced in (82). The function W approximates the solution of the invariance equations with an errorĒ such that This process can be iterated and at each step we solve the invariance equation exactly up to an order which is the double of the one we had for the initial approximation. Thus, if we assume that the initial guess W is such that the error in (78) satisfies The algorithm in this case is: Algorithm 7 (Computation of whiskers). Given F , ω as well as K, V s , µ and an approximate solution W such that with L ≥ 2 and W (θ, 0) = K(θ) and ∂ s W (θ, 0) = V s (θ). Perform the following calculations: Fast algorithms to solve difference equations with non constant coefficients In this section we present fast algorithms to solve for ∆(θ) the cohomology equation with non constant coefficients for given A(θ), B(θ) and η(θ) satisfying either A < 1, B −1 < 1 or A −1 < 1, B < 1.
Equations of this form appear in the Netwon step for whiskered tori (See the informal description in Section 1). Equations of this form also appear in the calcultion of the invariant splitting (see (67)-(68)).
We will present two algorithms. The first one is an iterative method with an accelerated convergence and the second one very fast (see Section A.1). The second one is only for the case of one-dimensional bundles and it is faster (computations are O(N ))(see Section A.2).
A.1. Fast iterative algorithms for the cohomology equation. In this section we will present a fast algorithm to solve equation (97) using iterative methods. We refer the reader to [HdlLS09] where a similar idea is used to compute iteration of cocycles.
The convergence of the algorithm is guaranteed by the contraction of A −1 and B. The cost of computing 2 n terms in the sum is proportional to n since it involves only n steps of the algorithm. The iterative algorithm is the following: Algorithm A.8 (Solution of difference equations with non constant coefficient). Given A(θ), B(θ) such that A −1 (θ) · B(θ) ≤ κ < 1, and η(θ) perform the following operations: The case when A < 1 and B −1 < 1 can be done similarly. In this case, we multiply (97) by B −1 (θ) on the LHS so that we obtain Before applying the iterative scheme we shift by −ω. In this way, we have The iterative algorithm in this case is Algorithm A.9. Given A(θ), B(θ) A(θ) B −1 (θ) ≤ κ < 1 and η(θ), perform the following operations: This algorithm gives ∆(θ + ω). Shifting it by −ω we get ∆(θ).
A.2. Fast algorithm for the 1-D cohomology equation with nonconstant coefficients. In this section we present an efficient algorithm for the one-dimensional version of equation (97). It is an adaptation of methods used in [Her83].
If we change the unknowns in (101) by W (θ) = C(θ)∆(θ), we are left with the equation Of course, if |ν| = 1, equation (102) can be solved in Fourier space. That is, we can obtain the Fourier coefficients of W as: and the solution is unique. Notice that whenever |ν| = 1, equation (102) involves small divisors, which is not the case for the iterative methods that will be discussed in the following section.
Taking log ν to be the average of log A(θ) − log B(θ), the problem reduces to solve for log C(θ) an equation of the form (23). Then C(θ) and ν can be obtained by exponentiation. Notice that log C(θ) is determined up to a constant. We will fix it by assuming that it has average 0.
Hence, we have the following algorithm: Algorithm A.10 (Solution of difference equations with non constant coefficient (1D)). Given A(θ), B(θ) and η(θ). Perform the following instructions: In this section, we provide the numerical algorithm to solve the invariance equation (11) and the modified one (13) using a Newton method analogous to the one described in Section 4.2.
The automatic reducibility can also be proved in this context (see [FdlLS09a]) and we provide here the algorithm only.

Appendix C. The order by order method for whiskers in Hamiltonian flows
In this section we present the result analogous to the one described in Section 5.1 to solve the invariance equation (17) for the whiskers in the case of Hamiltonian flows.
As in Section 5.1 we look for W as a power series W (θ, s) = ∞ n=0 W n (θ)s n , and match similar coefficients in s n on both sides of equation (17). For n = 0, one obtains which admits the solution W 0 (θ) = K(θ), where K is a parametrization of the invariant torus. For n = 1, we obtain from where we read that W 1 (θ) is an eigenfunction with eigenvalue −µ of the operator L ω L ω := ∂ ω − DX • K(θ). Again, we note that, multiplying a solution of (104) by a scalar b ∈ R, we also obtain a solution. See Remark 12.
Notice that, in this case, equation (105) can be solved provided that nµ is not in the spectrum of the operator L ω (this is a non-resonance condition which is clearly satisfied since the stable spaces are 1-dimensional). As in the case of maps, the previous equation can be solved using the large matrix method.
Appendix D. A Newton method to compute simultaneously the invariant torus and the whiskers for flows As in Section 5.2, we start with an initial approximation (W, µ) (resp. (W, µ, µ) for the invariance equation (17) (resp. (17)), that is X(W (θ, s)) − ∂ ω + µs ∂ ∂s W (θ, s) = E(θ, s), and we will look for an improved solution W → W + ∆ λ → λ + δ µ → µ + δ µ by solving the linearized equation Once again, we will use a reducibility argument similar to the automatic reducibility of Lagrangian tori. This will lead to a diagonal equation. Applying the operators D θ and ∂ s to equations (106), we obtain We expand the terms in (109) as a power series up to some order L and match coefficients of the same order on both sides of the equation. We use subindices to denote coordinates and superindices to denote the order.
Then, it is clear that the error E for the initial approximation W is such that E(θ, s) = n≥2 E n (θ)s n because the approximation is exact for the terms of order 0 and 1.