Quasi-periodic motions in a special class of dynamical equations with dissipative effects: a pair of detection methods

We consider a particular class of equations of motion, generalizing to n degrees of freedom the"dissipative spin--orbit problem", commonly studied in Celestial Mechanics. Those equations are formulated in a pseudo-Hamiltonian framework with action-angle coordinates; they contain a quasi-integrable conservative part and friction terms, assumed to be linear and isotropic with respect to the action variables. In such a context, we transfer two methods determining quasi-periodic solutions, which were originally designed to analyze purely Hamiltonian quasi-integrable problems. First, we show how the frequency map analysis can be adapted to this kind of dissipative models. Our approach is based on a key remark: the method can work as usual, by studying the behavior of the angular velocities of the motions as a function of the so called"external frequencies", instead of the actions. Moreover, we explicitly implement the Kolmogorov's normalization algorithm for the dissipative systems considered here. In a previous article, we proved a theoretical result: such a constructing procedure is convergent under the hypotheses usually assumed in KAM theory. In the present work, we show that it can be translated to a code making algebraic manipulations on a computer, so to calculate effectively quasi-periodic solutions on invariant tori. Both the methods are carefully tested, by checking that their predictions are in agreement, in the case of the so called"dissipative forced pendulum". Furthermore, the results obtained by applying our adaptation of the frequency analysis method to the dissipative standard map are compared with some existing ones in the literature.


Introduction
Why the Moon shows us always the same side? This is one of the most ancient scientific questions raised by the observation of the sky. The data made available by modern spatial missions clearly showed that the spin-orbit periodic motion is a rather common phenomenon in our solar system. Here, a p:q spin-orbit resonance means that the satellite turns on its spin axis p times while doing q revolutions around its star/planet. Actually, more than 20 planet-satellite pairs have been observed to stay into the 1:1 spin-orbit resonant state, while just one planet (Mercury) shows a different periodic behavior, because it rotates three times on itself during two complete revolutions around the Sun. A convincing explanation of the capture in resonance for the case of Mercury is provided in [27], where its present state is explained as a consequence of the fact that in the past the Mercury's orbit was much more eccentric. This scenario is discussed within the framework of a spin-orbit model including a dissipative force depending linearly on the relative angular velocity (for its introduction see also [35], [36], [49] and [55]). In this model, different periodic orbits can coexist and the measure of their basins of attraction can be evaluated both in a numerical and in an analytic way (see [18] and [6]). Within the different context of a viscoelastic model of the satellite, it has been recently shown that the capture into the 1:1 spin-orbit resonance is the generic final fate of such a dissipative system (see [3] and [4]).
In the last few years, the data about the rotational motion of some planets and satellites (e.g., Mercury, Titan and Europa) has been related to the study of their internal structure; this renewed the interest for the rotational dynamics of a non-rigid celestial body. In this context, an important role is played also by small oscillations around periodic orbits, which are also due to the perturbations exerted by other planets (see, e.g., [29]). Therefore, more and more sophisticated numerical tools are required to analyze this kind of weakly-dissipative systems.
In the present work, we adapt a numerical method and a semi-analytic one usually devoted to the study of Hamiltonian systems, in such a way to improve the description of the invariant attractors in the dissipative framework. Namely, they are the frequency map analysis and the constructive algorithm of the Kolmogorov's normal form, respectively.
The frequency map analysis has been originally designed by J. Laskar to study conservative systems (see, e.g., [40] and [41] for an introduction, while [53] is devoted to an interesting alternative approach). It is a powerful tool that has been applied to investigate the chaotic regions and those filled by invariant tori in several Hamiltonian systems (see, e.g., [25], [30], [43], and [54]) as well as in symplectic mappings (see [42]). In particular, the study of the variation of the fundamental frequencies allows to make a detailed cartography of the regular and chaotic regions in the Solar System (see [57]).
In the present work, we mainly focus on the model of the so called dissipative forced pendulum, whose the Newton equation can be written in the following form: where x ∈ T is an angle, ε is a (small) parameter and the potential U depends periodically both on x and the time t . Let us highlight the peculiar structure of the friction term η(ẋ−Ω) appearing in (1): it is linearly depending on the momentumẋ and it contains the so called external frequency parameter Ω . The Newton equations of both the dissipative spin-orbit model and the dissipative forced pendulum are of type (1), but the numerical explorations of the latter system require less computational resources than those needed by the former one. By the way, let us recall that the KAM-like theorems described in [19] and [60] apply to models described by the equation (1); moreover, dynamical systems including dissipative terms (which, of course, are not Hamiltonian) have been extensively studied in the last decades (see, e.g., [7] and [8]). Our numerical approach is based on the study of the regularity of the map Ω → ω 1 (Ω) , where the frequency ω 1 is related to the eventually existing quasi-periodic solution t → x(ω 1 t , t) of equation (1). We will show that our investigation method is very similar to that focusing on the action-frequency map, which is commonly studied, when conservative systems are considered; moreover, our approach applies also to dissipative mappings. In that context, a different frequency analysis has already been used in [22], for the study of the map relating the frequency ω 1 to the dissipation coefficient η , for various fixed values of the perturbing parameter ε .
One of the issues of our numerical method concerns the determination of the breakdown threshold (with respect to the small parameter ruling the size of the perturbing terms) of the invariant tori. This will allow us to compare our results with those given by other techniques, which have been widely tested in the literature. Among the known methods, Greene's technique (for an introduction, see [37] and [50]) provides the smallest uncertainity on the value of the breakdown threshold for symplectic mappings, but it is not so effective when dissipative terms are taken into account. This is due to the fact that the Greene's method is based on the calculation of a quantity (usually called residue), that is related to the eigenvalues of the monodromy matrix associated to a full cycle of a periodic orbit. Unfortunately, when the dissipation is introduced, each periodic orbit of fixed frequency ω 1 exists if and only if the external frequency parameter Ω ∈ [Ω ω 1 ; − , Ω ω 1 ; + ] ; moreover, when the order of resonance related to ω 1 is increased, the interval [Ω ω 1 ; − , Ω ω 1 ; + ] gets smaller and smaller. Considering the mid value of the interval [Ω ω 1 ; − , Ω ω 1 ; + ] is a good way to adapt the Greene's method to the dissipative standard map (as explained in [11] and [20]), but this interval is more and more difficult to locate for high order resonances; this limits the strength of the method. Another (recently established) technique evaluates the breakdown threshold, by studying the Sobolev norms of the function parametrizing the solution (see [14]). This approach apparently does not suffer any particular drawback, when dissipative terms are taken into account; therefore, it is able to determine the breakdown threshold with many significant digits (see [11]). Thus, the comparison with those results on dissipative mappings will represent a challenging test for our numerical method.
The semi-analytic method developed in the present paper strictly concerns with KAM theory adapted to dissipative systems. It is well known that the original versions of the KAM theorem ensure the existence of invariant tori filled by quasi-periodic orbits, in the context of both Hamiltonian systems and symplectic mappings, which are slightly perturbed with respect to some integrable approximations (see [39], [52] and [2]). In his first and last article on KAM theory, Kolmogorov pointed Celestial Mechanics as a field where such result could be naturally applied; his vision was definitely fruitful (see, e.g., [45], [28], [9], [10] and [17]). Actually, his proof scheme is based on the construction of sequences of canonical transformations and of the corresponding Hamiltonians, which are proved to converge (under suitable hypotheses) to the so called Kolmogorov's normal form (see [5] and [26]). In [32], it is proved that such a constructive algorithm can be rewritten according to a classical scheme (being O(ε r ) the size of the generating function of the r-th canonical transformation), in such a way to avoid the original quadratic convergence analogous to the Newton method (where the generating functions are O(ε 2 r ) at r-th normalization step). In [23], such a reformulation of the procedure constructing the Kolmogorov's normal form is shown to be highly effective in practical applications; moreover, it is well suited to locate invariant tori in Celestial Mechanics realistic problems (see [46], [31], [47], [48] and [58]). Let us also stress that the Kolmogorov's normal form can be used so to ensure the effective stability in a neighborhood of an invariant KAM torus, because the drift motion of the eventual diffusion can be estimated to be extremely slow (see [51] and [33]).
In [60], we have shown that equations of type (1) can be treated in the more general context of pseudo-Hamiltonian action-angle structures with n 1 + n 2 ≥ 2 degrees of freedom, where there are n 2 ≥ 0 external frequencies and the friction terms are linear and homogeneous with respect to the n 1 ≥ 1 actions; therefore, by using a technique of quadratic type, we proved the convergence of the algorithm constructing the Kolmogorov's normal form adapted to this pseudo-Hamiltonian framework, if the perturbation is small enough. In the present work, we reformulate the constructive procedure according to a classical formal scheme; moreover, we explicitly calculate the expansions of the Hamiltonians defined up to some finite normalization step, by algebraic manipulations on a computer. It is now rather common to say that such a method is semi-analytic, where we mean that we are going to use a constructive formal algorithm whose the convergence (at least for small perturbations) might be ensured by an analytic rigorous proof, but we limit us to show it, by directly checking the expansions produced on a computer.
For the sake of completeness, let us recall that recently the existence of quasi-periodic solutions for dissipative systems has been proved also in the more general context of conformally symplectic systems (see [12]). Such a result is based on a technique designed also to produce powerful applications to realistic models. Furthermore, that approach can be extended so to describe also the locally attracting dynamics in the neighborhood of the quasi-periodic solutions, although the proof scheme does not ensure the existence of any normal form (see [13]).
The present paper is organized as follows. In section 2, we define the models, which will be studied by our numerical investigations. In section 3, we adapt the frequency map analysis method to dissipative systems and, as a first stressing test, we compare our results about the breakdown threshold of invariant tori for the dissipative standard map, with those obtained by computing the Sobolev norms. Section 4 is devoted to the exploration of the dissipative forced pendulum model, by applying our adaptation of the frequency analysis. In section 5, the algorithm constructing the Kolmogorov's normal form is adapted to the general pseudo-Hamiltonian framework and it is applied to the dissipative forced pendulum model, so to check the agreement with the numerical results previously obtained in the present paper.
2 Introducing the models: dissipative standard map and forced pendulum In the present work we consider two simple systems: the dissipative standard map and the forced pendulum with dissipation. The standard map S ε : R × T → R × T is certainly the most famous symplectic map; here, we add a dissipation which is linear in the action variable. Thus, we consider the model defined by the equations where ε ≥ 0 is the perturbing parameter controlling the size of the perturbation, η ≥ 0 is the friction coefficient ruling the dissipation rate and Ω is an external forcing frequency. Let us remark that in the conservative case (i.e., when η = 0) the formula above covers also the definition of the usual standard map S ε . Moreover, in the unperturbed case (i.e., when ε = 0), the set {y = Ω , x ∈ T} is an invariant global attractor of the dynamics and Ω is also the frequency value of the angular motion on that torus. The dissipative standard map has been widely studied, like for example in [22], where it is defined as follows: In that case, the obvious correspondence between variables and parameters appearing in (2) and in (3) is given by the equations x = 2πX , y = 2πY , η = 1 − b and ηΩ = 2πc .
The second system considered here is the dissipative pseudo-Hamiltonian model of the forced pendulum. In order to properly define it, let us introduce the autonomous Hamiltonian describing the forced pendulum (with the variable q 2 playing the role of time) where (p 1 , p 2 ) ∈ R 2 , (q 1 , q 2 ) ∈ T 2 and ε is a small positive parameter. Let us simplify the notation, by introducing the "Hamiltonian vector field operator" V H , which acts on a dynamical function g : R n × T n → R (being n a generic number of degrees of freedom) so that Therefore, our pseudo-Hamiltonian model of the dissipative forced pendulum is described by the following equation: where the meaning of the symbols η and Ω is the same as in (2). The dissipative forced pendulum introduced above is substantially defined by a system of three differential equations depending on the variables q 1 , q 2 , p 1 ; the evolution of the action p 2 is actually irrelevant, because it does not have any influence on the behavior of the other variables. Moreover, once the law of motion t → q 1 (t), q 2 (t), p 1 (t) is known, the function t → p 2 (t) can be determined by computing an integral. When one is interested in investigating numerically the behavior induced by the differential equation (6), it is natural to consider the corresponding Poincaré map. This allows us to decrease the number of degrees of freedom (and to reduce the numbers of variables from 3 to 2), by sampling the state of the system at times which are multiple integers of the period of the variable q 2 , that is 2π. In other words, we are going to study the Poincaré map M ε , η , Ω : R × T → R × T , that is defined so that M ε , η , Ω (p 1 , q 1 ) = Φ 2π ε , η , Ω (p 1 , 0, q 1 ) , where Φ δ ε , η , Ω : R × T 2 → R × T 2 is the δ-time flow induced by the equation (6) and we do not take into account its effect on p 2 .
One can easily check that in the conservative case M ε , 0 , Ω is a symplectic map. Let us recall that, apart a further rescaling of the parameters, the dissipative standard map is nothing but a very rough approximation of M ε , η , Ω that is produced by a single step of the so called semi-implicit Euler method, covering a time interval equal to 2π .
From a practical point of view, in all the numerical experiments described in the present paper, the Poincaré map M ε , η , Ω is approximated by a numerical integration of the equations of motion (6), using the Taylor 1 method (see [38]). In our tests, such a software package is able to numerically integrate the flow Φ 2π ε , η , Ω , performing less than 30 steps, each of them is affected by an uncertainity not greater than the round-off error on double type variables of the C programming language. The precision of the integration scheme could be further improved, by using long double type variables or multiple precision arithmetic, that can be very well performed also by using the TIDES software package (see [1] for an introduction). We consider that our numerical results should be very slightly modified by such a further improvement and, so, it has not been implemented.
An important feature of the dissipative systems is that they need a relaxation time before converging to the invariant attractor. From the computational point of view, this means that a certain number W of preliminary iterations is necessary, in addition to those required by the frequency map analysis. In order to provide a criterion for the choice of the value of W , let us consider the unperturbed case of the dissipative standard map (2) and assume the initial value of the ordinate is y 0 ; then, one can easily check that the sequence of the iterated points is such that y n = (1 − η) n (y 0 − Ω) + Ω , ∀ n ∈ N . In a numerical experiment, the value of y 0 is determined so that the initial point (x 0 , y 0 ) is rather close to the wanted invariant attractor. Therefore, we define W so that (1 − η) n is at most of the order of the machine precision ∀ n ≥ W , i.e., being ⌈α⌉ the smallest integer greater than or equal to α ∈ R .
3 Adapting the frequency map analysis to dissipative systems

Frequency map analysis for Hamiltonian systems: a short overview
Since our investigation approach for dissipative systems is strongly reminiscent of the method designed by Laskar to study conservative systems (see, e.g., [40] and [41]), we think that it is convenient to recall some of its features in the present subsection. This will allow us to introduce our adaptation for dissipative systems in a more natural way. Let us consider an n-d.o.f. quasi-integrable system, described by an analytic Hamiltonian where (I, θ) ∈ G × T n (being G ⊂ R n an open set) are action-angle variables. According to KAM theory (see, e.g., [56]), if the following conditions are satisfied: (A) the integrable part h(I) is non-degenerate (i.e., the determinant of the hessian of h is different from zero ∀ I ∈ G), (B) the parameter ε is small enough; then, there exists a diffeomorphism Ψ : B × T n → G × T n having the following properties: (I) Ψ(ω, ϕ) is invertible and it is C ∞ with respect to ω ∈ B and analytic in ϕ ∈ T n , (II) there is a Cantor set B ε ⊂ B such that for each (diophantine) frequency ω ∈ B ε the law of motion t → (I(t), θ(t)) = Ψ(ω, ωt) is a solution of Hamilton's equations on an invariant (KAM) torus, (III) when G is bounded, the Lebesgue measure of B \ B ε tends to zero for ε → 0 . Let us recall that here the non-degeneracy condition on the integrable part h can be replaced by the so called isoenergetical non-degeneracy (see, e.g., [16] for a definition). Moreover, while practically doing numerical explorations, the (very restrictive) smallness condition on the parameter ε can be ignored, because it is known that in a neighborhood of a generic invariant KAM torus, there is a canonical transformation leading the Hamiltomian to a form such that the above conditions (A) and (B) are satisfied (see [51]).
Let t → z(t) a signal (depending on time) in the complex plane, where z = z(I, θ) is a function defined on the phase space. Let us suppose that the law of motion t → (I(t), θ(t)) is quasi-periodic and is characterized by the frequency vector ω , being its corresponding orbit on an invariant KAM torus; then the property (II) above allows us to assume that the signal t → z(I(t), θ(t)) is as follows: The basic software package implementing the numerical analysis of the fundamental frequencies (see, e.g., [42]) must allow us to calculate a suitable truncation of the expansion above. Actually, the values of the frequencies ζ l are numerically found by looking for the local maxima of the following function: where w is a weight function, i.e., an analytic, non-negative, even map w : is meant to be a time interval where the signal t → z(t) has been preliminarly computed (usually, by a numerical integration approximating the law of motion t → (I(t), θ(t)) ). In practical applications, it is natural to sample such a signal with N uniform subintervals of [t 0 − T, t 0 + T ] , being their width equal to ∆ = 2T /N . Therefore, the integral appearing in formula (11) can be approximated by using the trapezoidal rule or a similar quadrature formula. Of course, one expects that the numerical calculation of the Fourier decomposition (10) becomes better and better when the value of T increases. We suggest to the reader (once again) the reviews [40] and [41] for the careful discussion about the accuracy of the numerical results and their dependence on the parameter T , ∆ and the weight function w . In all our numerical experiments (described below), we used the Hanning window filter w(u) = 1 + cos(πu) ; here, we just recall that with this kind of weight function, the difference between the computed value of the frequency vector ω and the true one is A first test to check if an orbit lies on an invariant KAM torus can be made by controlling that the Fourier decomposition (10) of a corresponding signal holds true, within the limitations due to the unavoidable numerical errors. Let us remark that the frequency vector ω is not given a priori, but its detection is often not so difficult, in practical applications. For instance, let us consider a quasi-integrable Hamiltonian of the type (9) and satisfying the conditions (A) and (B); ∀ j = 1, . . . , n , let us study the signal z(t) = I j (t) exp iθ j (t) ; then the pointσ T corresponding to the absolute maximum of the function (11) gives an approximation of ω j , that gets more and more accurate for T → ∞ .
Another natural numerical investigation concerns the local regularity and invertibility of the action-frequency map I → ω(I) such that ω(I) = Ψ −1 (I, θ) for any fixed value of θ ∈ T n . The frequency map analysis mainly aims to obtain directly, in a numerical manner, the map I → ω(I) . The procedure, can be summarized as follows: we first arbitrarily fix the initial values of the angles θ 0 ; we pick up the initial actions I 0 from a regular grid J of values. For each initial condition, we consider the corresponding motion law t → (I(t), θ(t)) and we analyze the n signals z . . , n ; for each signal, we find the value of ω j corresponding to the absolute maximum of the function (11). Following this procedure, we can then calculate the frequency ω = ω(I 0 ) for all the initial values of the actions I 0 ∈ J .
The analysis of the frequency map can distinguish among three different dynamical situations: (a) when the values of the initial conditions (I 0 , θ 0 ) are such that the corresponding motions are chaotic ∀ I 0 ∈ J , then the map I 0 → ω(I 0 ) looks highly irregular; (b) when the initial conditions are such that the orbits are on the regular manifolds inside a resonant region (these are the so called "librational" maximal tori in the neighborhood of a stable equilibrium point or an elliptic lower dimensional torus), some components ω j of the frequency vector are constant while I 0 is changed; (c) when the initial conditions are in a region (nearly) filled by KAM tori, a thin enough enlargement of the frequency map highlights a quasi-linear and invertible behavior. This is in agreement with the property (I) (that is described above and proved in [56]) joined with the approach described in, e.g., [51]: the signature of the existence of a KAM torus is the local regularity and invertibility of the action-frequency map I → ω(I) such that ω(I) = Ψ −1 (I, θ) for any fixed value of θ ∈ T n .
These three different regimes can be sharply highlighted with some numerical experiments on symplectic maps. For instance, let us consider the standard map S ε , as it is defined by formula (2) when η = 0 . The results plotted in Figure 1 are obtained by analyzing the signal n → z(n) with z(n) = y n exp(ix n ) , where the pair (y n , x n ) is obtained by n iterations of S ε , starting from the initial condition (y 0 , x 0 ) . Let us remark that we are assuming that the signal is sampled in a trivial way, so that the "elapsed time" between an iteration of the standard map and the next one is ∆ = 1 . By looking at the definitions in (10) and (11), one can easily realize that changing the definition of ∆ would imply a harmless rescaling of the found value of the frequency ω by a factor 1/∆ . We can appreciate that the archetypical behaviors described at the points (a)-(c) are clearly detected by the numerical experiments, whose results are plotted in Figures 1a-c, respectively.

Frequency map analysis for dissipative systems
In some dissipative systems, there is just one global attractor for the dynamics. For instance, if we consider the unperturbed case of the dissipative forced pendulum, that is described by the equation (6) setting ε = 0 in (4), the solution for the motion of the action p 1 can be written as Figure 1: Frequency map analysis of the standard map S ε , with ε = 0.97 . In the plots above, three different sets of orbits are considered; each orbit is computed starting from an initial condition of the type (x, y) = (0, y 0 ) and its value y 0 is reported in abscissa.
The corresponding value of the main frequency ω 1 has been calculated by analyzing the signal n → y n exp(ix n ) with 0 ≤ n ≤ N ; see the text for further definitions and details.
To obtain the values plotted in Figure 1a we fixed N = 2 18 , while for Figures 1b-c N = 2 16 . Figure 1a refers to some orbits in a neighborhood of the chaotic zone related to the resonance 610/987 ; in the central part of Figure 1b, some chains of regular islands contouring the stable periodic orbit of frequency 2π 377/610 are considered; Figure 1c focuses on a neighborhood of the "golden" invariant torus of frequency ω 1 = 2π( By the way, let us remark that the Hamiltonian (4) when ε = 0 describes nothing but a rotator plus a clock. Looking at the equation above, it is obvious that p 1 (t) → Ω for t → ∞ . The motion law on the global attractor is given by the following equations: In the unperturbed case with ε = 0 , one can provide examples of weakly dissipative systems, where there are more than one single attractor (see, e.g., [18]). Since each basin of attraction usually contains open sets of the phase space, one immediately realizes that the study of the map ω(I) = Ψ −1 (I, θ) loses sense for dissipative systems. In fact, there are many initial values of the actions I corresponding to the same final frequency ω , therefore, the action-frequency map is obviously not invertible even when the final attractor is an invariant torus. However, the trivial example of the unperturbed case can help us to explain the simple idea underlying our new approach. The solution (13) highlights that the frequency of the invariant attractor is ω 1 =q 1 = Ω , thus, when ε = 0 the map Ω → ω 1 (Ω) is obviously regular and invertible, because it is nothing but the identity. It is natural to expect that the properties of being both regular and invertible are preserved by such a map, also in the perturbed case for small values of ε. Actually, this is guaranteed by the main result of [19], at least for what concerns systems of type of the dissipative forced pendulum, that is defined by the equations (4)- (6). In fact, when the perturbing terms are small enough, Theorem 1 of [19] claims also the following relation between the frequency ω 1 of the quasi-periodic motion (on an invariant torus) and the external forcing one Ω : when ω ∈ D γ,τ , being D γ,τ the set of diophantine numbers such that for some fixed values 0 < γ < 1 and τ ≥ 1 . Among the diophantine numbers, a special class are the "noble" numbers, having their continued fraction expansion ending with only 1; in particular we will consider here the "golden number" φ = [1; 1 ∞ ] = ( √ 5+1)/2 ≃ 1.618 . Moreover, the function ω 1 → Ω(ω 1 ) is Whitney 2 C ∞ on the Cantor set D γ,τ . Therefore, the equation (14) leads us to conclude that the map Ω → ω 1 (Ω) is regular and locally invertible in the neighborhood of an invariant torus, if ε is small enough. As it has been claimed by Celletti and Chierchia in remark (v) of section (1.1) of [19], it is expected that such kind of results can be extended to systems with more degrees of freedom. Thus, we conjecture that when a dissipative system is governed by the following equations of motion and its Hamiltonian part H satisfies the hypotheses (A)-(B) (described in the subsection 3.1), then there exists a diffeomorphism Ξ(ω, ϕ) which satisfies the same properties(I)-(III) (holding for the diffeomorfism Ψ of the conservative case), with the action vector I replaced by the external forcing frequency vector Ω . The previous discussion leads us to conclude that the frequency map analysis can be adapted to the dissipative systems of the type (16), by simply using the external forcing frequency vector Ω instead of the initial value of the action I 0 . Actually, here we can consider a set of motions, each of them corresponds to a different value of Ω , got from a regular grid A . For the sake of simplicity, we postpone to the next sections further details, about the description of the procedure calculating the corresponding frequencies of the motion on the attractors. By analogy with the points (a)-(c) of subsection 3.1, we guess that also here the analysis of the frequency map can distinguish among three different dynamical situations: (a ′ ) when the values of the external forcing frequency vector are such that (for some set of initial conditions) the corresponding orbits converge to a strange attractor ∀ Ω ∈ A , then the map Ω → ω(Ω), should look highly irregular; actually, both a strange attractor and a chaotic orbit have fractal dimension larger than the degrees of freedom (see Figure 5 in [21]), therefore, it is natural to argue that the behavior will be the same as for the chaotic motions in Hamiltonian systems (see point (a) of the previous subsection 3.1); (b ′ ) when the external forcing frequency vectors are such that some orbits converge to attractors that are regular manifolds inside a resonant region, some components ω j of the frequency vector are constant as Ω varies; for instance, this statement is well supported by the numerical experiments on the dissipative standard map; in that case, let us recall that each periodic orbit of fixed frequency ω 1 exists if and only if the external frequency parameter Ω ∈ [Ω ω 1 ; − , Ω ω 1 ; + ] ; moreover, when the order of resonance related to ω 1 increases, the interval [Ω ω 1 ; − , Ω ω 1 ; + ] gets smaller and smaller (see, e.g., [20]); (c ′ ) when the values of the external forcing frequency vector are in a region of the regular grid A such that the corresponding attractors are invariant KAM tori, then a thin enough enlargement of the frequency map should highlight a quasi-linear and invertible behavior. This is in agreement with our conjecture about the systems governed by the equation of motion (16), when its Hamiltonian part H satisfies the hypotheses (A)-(B) (described in the previous subsection 3.1).
Finally, let us remark that it is natural to expect that the behaviors described at the previous points (a ′ )-(c ′ ) should hold, also when the non-degeneracy condition (A) of subsection 3.1 is replaced by the weaker one we considered in [60] (i.e., condition (b) of theorem 3.1). Moreover, it is natural to guess that the analysis of the frequency map should highlight the same situations also for dissipative maps, that are obtained as a Poincaré map of the continuous flow induced by equations of motion of the type (16).

Numerical experiments on the dissipative standard map
The interpretation of the frequency map described in the previous subsection allows us to investigate the breakdown threshold of invariant tori for dissipative maps. This allows us to submit our method to some challenging test, because we can compare our results with some existing ones in literature.
Let us focus on the dissipative standard map (2). Figure 2 shows the frequency maps Ω → ω 1 (Ω) for ε = 0.9719 and ε = 0.9721; in both cases the friction parameter η has been set equal to 0.1 . Those maps have been drawn by analyzing the signals of the type n → z(n) with z(n) = y n exp(ix n ) , where the pair (y n , x n ) is obtained by n iterations of the dissipative standard map (2), starting from the initial condition (y 0 , x 0 ) = (Ω, 0) . Each point plotted in Figure 2 is actually related to a single analysis of an orbit corresponding to the value of Ω reported in abscissa; such an analysis considers all the values of index n ranging in [W, N + W ], with W given by formula (8) and N = 2 19 = 524288 . This means that we perform a "waiting" number W of preliminary iterations (that are needed to let the orbit approach very closely an invariant attractor), before starting the calculation of the value ω 1 (Ω) corresponding to the absolute maximum of the map (11). For definiteness, let us recall that the integral appearing in formula (11) is approximated numerically by using the trapezoidal rule with N subintervals (all with the same width) in [t 0 −T, t 0 +T ] , being t 0 = W + N/2 and T = N/2 ; moreover, the "weight" function w is the Hanning window filter w(u) = 1 + cos(πu) .
Before discussing the results, we need to make some preliminary remarks about the definition of the frequencies. In the case of the dissipative standard map, as for any discrete time map, frequencies that differ by an integer multiple of 2π are equivalent, in the sense that their dynamics are undistinguishable. In the special case of the dissipative standard map, it is also obvious that the dynamics is also 2π-periodic in the action. These two properties mean that tori whose frequencies differ by any multiple of 2π are equivalent,  Figure 2: Frequency analysis for the dissipative standard map (2), with the friction parameter η = 0.1. The range of abscissas (related to the external frequency values of Ω) has been determined so to focus on a neighborood of the golden value ω 1 /2π = 2 − φ = (3 − √ 5)/2 . The left plot refers to the case ε = 0.9719 and the right one to ε = 0.9721 . and in particular all the "golden tori" with frequencies 2πφ , 2π(φ − 1) , 2π(φ − 2) , . . . have the same shape and break in the same way. Moreover, for the dissipative standard map, the dynamics is also invariant when changing x → −x, y → −y and Ω → −Ω; this also implies that tori having opposite frequencies ω 1 and −ω 1 also behave in the same way. We thus decided to perform our numerical experiments on the torus with frequency ω 1 /2π = 2 − φ = (3 − √ 5)/2 ≃ 0.381966 , which is the only golden torus with positive frequency ω 1 in [0, π] . In Figure 2 (both on the left and on the right), the dashed horizontal lines are located in correspondence to the value of the ordinate equal to the "golden torus" frequency ω 1 = 2π(2 − φ) = 2π[(3 − √ 5)/2] . According to our discussion in the previous subsection 3.2, we are led to conclude that the attractor related to the golden mean frequency exists, if the map Ω → ω 1 (Ω) looks regular (i.e., quasi-linear) in a small neighborhood of the intersection with the dashed line; otherwise, when the map shows sudden jumps where it is crossing the dashed line, then that invariant torus does not exist. The left panel of Figure 2 clearly shows that "golden torus" still persists for ε = 0.9719 , while the right panel makes evident that it is destroyed when the parameter ruling the perturbation is ε = 0.9721 . This allows us to conclude that the breakdown threshold ε c should be in the interval (0.9719 , 0.9721) when ω 1 = 2π(2 − φ) and η = 0.1 . By the way, let us remark that in Figure 2 the large "plateaus" appearing in the left plot and in the right one correspond to the resonant values of ω 1 /(2π) equal to 2584/6765 , 1597/4181 (left), 987/2584 , and 1597/4181 (right).
In Table 1 we collect some results obtained by applying our method to compute the critical values ε c of the breakdown threshold for a pair of invariant tori and a few different values of the friction parameter η . Let us precise that we repeatedly use the same procedure described both for the case of the golden torus with ω 1 /(2π) = 2 − φ and for that concerning ω 1 /(2π) = [0; 2, 5, 3, 1 ∞ ] ≃ 0.4567 . The cases studied here can be directly   Tables I-III in [11]). Since the results listed in Table 1 are in agreement with both those based on the computation of the Sobolev norms and those obtained by applying the Greene's method, we consider that this comparison strongly support the validity of our approach, that was heuristically motivated in the previous subsection 3.2.
A more detailed comparison of the results with those provided in [11] highlights that the most performing method (to determine the breakdown threshold) is that based on the computation of the Sobolev norms; in fact, it provides the largest number of significant digits (about five). We emphasize that our approach can be nicely visualized (as in Figure 2), but, as an evident drawback, it is not easy to make the whole procedure very automatic. For instance, the determination of a suitable range of abscissas often requires many trials and errors; moreover, the numbers of trials significantly increases when a high precision is required. This is because we limited ourselves to compute the breakdown threshold up to the third significant digit in all the cases listed in Table 1, except that illustrated in Figure 2.

Numerical results about the dissipative forced pendulum
Let us now focus on the equation of motion (6) for the dissipative forced pendulum, where the Hamiltonian part H ε is defined in (4). In this section, our study aims also to determine the values of some parameters (for instance, the breakdown threshold ε c and the external forcing frequency Ω), which must be known in advance, before starting any explicit calculation of Kolmogorov's normalization algorithm, that will be discussed in the next section.

Breakdown of invariant tori in the dissipative forced pendulum
It is natural to adopt exactly the same approach used in subsection 3.3 to study the dissipative standard map, in order to investigate the behavior of the Poincaré map (7), related to the dissipative forced pendulum. In practice, after having somehow fixed the values of the parameters ε , η and Ω , we can produce frequency maps by analyzing the signals of the type n → z(n) with z(n) = y n exp(ix n ) , where the pair (y n , x n ) is obtained by n iterations of the dissipative map M ε , η , Ω , that is defined in (7). Each signal n → z(n) is analyzed so to determine the value ω 1 of the absolute maximum point of the map (11). Such an integral is numerically approximated in the same way as we did in subsection 3.3.
In particular, the endpoints of the interval [t 0 −T, t 0 +T ] are fixed so that t 0 = 2π(W +N/2) and T = 2π(N/2) , where W is given by formula (8) and N = 2 16 = 65536 . As a first numerical investigation about the dynamics of the dissipative forced pendulum, we study the breakdown of the invariant "golden torus". Let us remark that since we decided to sample the continuous dynamics with a timestep ∆ = 2π , the relation (14) still holds, but with frequencies ω 1 now in the interval [−0.5, 0.5]. The frequency maps Ω → ω 1 (Ω) for ε = 0.0371 and ε = 0.0373 are showed in Figure 3, in both cases the friction parameter η is equal to 0.1 . Both in the left plot and in the right one, the dashed horizontal lines are located in correspondence to the value of the ordinate equal to the frequency ω 1 = 2 − φ = (3 − √ 5)/2 ≃ 0.381966 . According to the discussions in subsections 3.2-3.3, we can provide a clear interpretation of the results illustrated in Figure 3: the breakdown threshold of the wanted invariant torus is ε c = 0.0372 ± 10 −4 .
In the case of the dissipative forced pendulum, we think that it is interesting to study the dependence of the breakdown threshold ε c (ω 1 , η) on the friction coefficient η . More precisely, we want check if the function η → ε c (ω 1 , η) behaves according to the KAM-like analytical estimates. The discussion about some functional properties of the theoretical threshold ε ⋆ (which depends on many parameters characterizing the system) is included in sect. 3 of [60] and it is reported below, for the sake of completeness.
(A) There is a range of "small" values of the friction parameter, with 0 ≤ η ≤ η ⋆ 1 , for which ε ⋆ is a constant.
(C) The value of η ⋆ 2 depends on both the non-resonance assumptions about the frequency vector and those guaranteeing the non-degeneracy; actually, η ⋆ 2 ≤ 10 5 η ⋆ 1 and the limit case η ⋆ 2 = 10 5 η ⋆ 1 should hold true just when the latter conditions are much weaker than the former ones.
Let us recall that the previous points (A)-(D) cover also the so called anti-dissipative case with η < 0 , because all the analytical estimates depend just on the absolute value of the friction parameter.
Here, we limit ourselves to investigate the function η → ε c (ω 1 , η) just in the case of the "golden torus" with frequency ω 1 = 2 − φ = (3 − √ 5)/2 (adopted for consistency with the experiment on the dissipative standard map), so to compare the predictions for the theoretical breakdown threshold ε ⋆ with the behavior of the numerical one ε c . Figure 4 includes (on the left) a table with some values of the correspondence η → ε c , numerically determined by applying our frequency analysis approach to the study of the dissipative forced pendulum. When the friction coefficient η gets smaller and smaller then the corresponding value of ε c seems to converge toε c ≃ 0.0275856 , that is the breakdown threshold of the golden torus in the conservative case (see, e.g., [24]). Since the map η → ε c looks regular and we can guess that it is an even function (let us recall that the value of ε ⋆ is preserved by the simmetry η → −η ), the fact that the scaling law for η → 0 , is clearly superlinear suggests that η → ε c (2 − φ, η) has a quadratic minimum in the origin. This is in agreement with the behavior of the theoretical breakdown threshold ε ⋆ described at points (A)-(B). Moreover, the plot in logarithmic scale (on both axes) of Figure 4 highlights the fast growth of the numerical breakdown threshold ε c when larger values of the friction coefficient are considered. Its behavior when η ≫ 1 definitely need additional investigations, so to check if it fits well also with the predictions (C)-(D).

Numerical determination of the forcing frequency
In order to perform explicitly the algorithm constructing the Kolmogorov's normal form (as described in the next section) for dissipative systems, we need to preliminarly determine the external forcing frequency. To fix the ideas, we limit ourselves to consider again the dynamics of the dissipative forced pendulum. The aim of this subsection is to determine, for a fixed invariant torus, the corresponding value of the parameter Ω (appearing in the equation (6), where the Hamiltonian part H ε is given in (4)). As discussed above, the frequency map provides the frequency ω 1 of the quasi-periodic motion on an invariant torus as a function of the parameter Ω . Now, the problem is the following: we fix the frequency ω 1 = ω * 1 related to an invariant torus and we need to approximate numerically the corresponding forcing frequency Ω * . Thus, denoting again the frequency map by Ω → ω 1 (Ω) , we want to find numerically the solution Ω * of the equation This requires to invert the frequency map, or, more simply, to find the real zero of the function f (Ω) = ω 1 (Ω) − ω * 1 . For this purpose, we implement explicitly a Newton's method which, as it is well known, is an iterative method to find numerically the solutions for this kind of problems. Let us stress that we expect to find a locally unique solution Ω * of the equation ω 1 (Ω) − ω * 1 = 0 , because we obviously apply Newton's method for values of the parameter ε smaller than the breakdown threshold ε c of the invariant torus related to the frequency ω * 1 . Therefore, in a neighborhood of the unknown value Ω * , the function ω 1 (Ω) has a quasi-linear behavior, which looks strictly monotone, except in the resonant zones (see the left plots in Figure 2-3 and the discussions about them). Thus, if the initial approximation belongs to the region about the solution where the map Ω → ω 1 (Ω) looks mostly quasi-linear and monotone, Newton's method is expected to be very efficient.
In order to be more definite, in the following we describe our procedure in detail. We denote withΩ (j) , the j-th approximation of the solution Ω * ; then, the single step of Newton's algorithm, applied to f (Ω) = ω 1 (Ω) − ω * 1 prescribes that the next approxima-tionΩ (j+1) is given bỹ where the derivative f ′ Ω (j) is replaced by the finite difference ω 1 (1 + α)Ω (j) − ω 1 Ω (j) / αΩ (j) and α is a small parameter to be conveniently fixed so to ensure the numerical stability of this procedure. Of course, in formula (17) the values of ω 1 Ω (j) and ω 1 (1 + α)Ω (j) are numerically calculated, by using the frequency analysis. We stop the iterations when the relative correction on the value of Ω is below a fixed precision β , i.e., when where β is a small parameter that can be conveniently chosen so to be not much greater than the machine precision (let us recall that this is about 2.2 × 10 −16 for the standard double precision type numbers). Let us remark that it is very unlikely that at some step the finite difference ω 1 (1 + α)Ω (j) − ω 1 Ω (j) / αΩ (j) is close to 0 , because it occurs that both (1 + α)Ω (j) and Ω (j) are in the same resonant region. In fact, the sizes of the resonant "plateaus" are smaller and smaller, when approaching the invariant torus; this fact can be seen in the left plots of Figures 2-3 and it has been clearly shown in the conservative framework (see, e.g., the numerical investigations in [44]). Thus, we limited ourselves to insert a test in our code, to stop the running if the finite difference above is too small. This event is so rare that it never happened in our calculations; from a practical point of view, in such a case of failure, one has to look for a better initial approximationΩ (0) before restarting the procedure.
For instance, let us discuss an explicit case. We want to determine the value of the external forcing frequency such that the equations of motion (6) has the golden torus as an invariant attractor, when the values of the parameters are fixed so that Let us recall that the value of the small parameter ε is chosen smaller than the breakdown threshold related to ω * 1 = 2 − φ = (3 − √ 5)/2 and η = 0.1 (see the corresponding value of ε c in the table appearing in Figure 4). We fix α = 10 −6 and β = 10 −15 and we start the Newton's algorithm takingΩ (0) = ω * 1 as initial approximation. In this case, the algorithm ends successfully after just 5 steps with As an internal test of our result, we performed the decomposition of the Fourier spectrum as in formula (10). Actually, we considered the motion on the invariant attractor for the dissipative forced pendulum defined by equation (6), with ε = 0.03 , η = 0.1 and Ω = Ω * ≃Ω (5) , with the value ofΩ (5) given in (19). Since the golden torus is expected to be the invariant attractor, we tried to explain every frequency ζ l as a linear combination of the components of the vector ω = (3 − √ 5)/2 , 1 . The relevant quantities involved in the decomposition of the Fourier spectrum are reported in Table 6.2 of [59], where the numerical numerical results definitely show that the invariant attractor is the golden torus related to the frequency ω * In order to avoid the proliferation of too many symbols, starting from equation (26) we do not introduce another set of variables for the new coordinates introduced after each canonical transformation; this is done by abuse of notation. Moreover, in the functional equation (28) the Lie series operator exp L χ (r) 1 appears, where L χ g = {g, χ} , being {·, ·} the classical Poisson bracket, g a generic function defined on the phase space and χ any generating function.
Since we point to a Hamiltonian part of type (23), first, we determine the generating function χ (r) 1 so to remove both the main perturbing terms of degree 0 and those that are linear with respect to the actions but do not depend on the angles. Thus, we solve with respect to X (r) (q) and ξ (r) the equations where the n 1 ×n 1 matrix C (r−1) is such that 1 2 (p 1 , . . . , p n 1 )·C (r−1) (p 1 , . . . , p n 1 ) T = f we can easily write the solution of the first homological equation appearing in (29), i.e., Let us emphasize that the solution above is well defined when the friction coefficient η = 0 ; in the conservative case (i.e., η = 0), it is enough to use the non-resonance condition (41), that will be explicitly adopted to solve the second homological equation.
We must now provide the expressions of the functionsf whereĤ (r) is defined by the functional equation (28). To this aim, we will redefine many times the same quantity without changing the symbol. In our opinion, such a repeated Let us recall that the canonical transformation K (r) inducing the Kolmogorov's normalization up to the step r is explicitly given by This conclude the r-th step of the algorithm that can be further iterated. Let us stress that the next step can be completely carried out, if both the non-degeracy condition (31) and the non-resonant inequality (41) still hold true, when the index r is replaced with r + 1 . Actually, for what concerns the former assumption, this is usually ensured by requiring that both the quadratic part of the initial Hamiltonian is non-degereate (i.e. det(C (0) ) = 0) and the parameter ε is small enough. Moreover, the latter non-resonance condition is usually provided for all indexes r , by choosing the frequency vector ω so that it is Diophantine. Let us emphasize that terms having different orders of magnitude with respect to the small parameter ε are not separatedly kept in our expansions. In particular, we prescribed to perform operations like the "reordering of the terms", which explicitly requires to sum contributions corresponding to the same polynomial degree and Fourier harmonic, but with different orders in ε . The main advantage of this formulation is to save most of the memory occupation, when the algorithm is translated in any programming language (see the discussion at the end of section 4.1 of [31]). As a consequence of such a gain in the memory handling, more normalization steps can be performed (and in a faster way); let us stress that this can definitely improve the final accuracy of the results. From a practical point of view, when H ε is expressed in the initial form (21), the parameter ε must be replaced by its numerical value. Therefore, in all the expansions written in the present subsection (whenever they are converging in some suitable domains), the sup-norm of the functions of type f (r,s) l andf (r,s) l is geometrically decreasing with respect to both the polynomial degree l and the index s , that is related to the their trigonometric degree sK .

Semi-analytic results
Perhaps, it could be astonishing that both the equations of motion (36) and the Hamiltonian part (43) defined at the r-th normalization step have exactly the same structure as those introduced by the previous step, which are written in (24) and (25), respectively. Indeed, performing the algorithm described in the previous subsection is advantageous, because the unwanted Hamiltonian terms of degree 0 and 1 in the actions get smaller and smaller, when r is increased (under the usual KAM hypotheses). This implies that the algorithm is successful if and only if also the generating functions decrease with r (recall the equations (29) and (39)); this remark can be easily translated in a numerical test concerning the construction of the normal form.
The behavior of the sequence of the external forcing frequency vectors Ω (r) r≥0 deserves a particular discussion. Let us recall that the main theorem in [60] actually proves the existence of a pair of objects: an initial frequency vector Ω (0) and a canonical transformation ψ (∞) such that the equations of motion (20) are conjugated to those in (22), where the Hamiltonian part is in the Kolmogorov's normal form (23). Thus, the proof scheme determines Ω (0) a posteriori, i.e., as a result of the normalization procedure (actually, here we followed the approach originally designed in [19]). Of course, this is unpractical when we focus on comparisons with numerical results, because any integrator of the equations of motion (20) requires that Ω (0) must be known in advance. This is the main reason why, in subsection 4.2 we developed a Newton method based on the frequency analysis, so to provide a good approximation of the initial Ω (0) , corresponding to the fixed angular velocity vector ω of the quasi-periodic motion on the wanted invariant torus. If such a vector Ω (0) would be perfectly determined, then the algorithm constructing the conjugacy to (22) (where any external frequency is not appearing) requires that Ω (r) → 0 for r → ∞ . This condition is going to be numerically tested.
From a practical point of view, let us proceed to a further study of the equations of motion (6), where the Hamiltonian part H ε is given in (4) and the values of the parameters are fixed so that For such a system, let us recall that the frequency analysis results reported in subsection 4.2 clearly show the existence of an attracting invariant torus, characterized by a quasi-periodic motion related to the golden frequency ω * 1 = 2 − φ = (3 − √ 5)/2 . In our opinion, the explicit construction of the normal form related to that torus is rather challenging, because the system is so perturbed that the small parameter ε is larger than the breakdown threshold value for the conservative case (i.e.,ε c ≃ 0.0275856) and not so far from the same threshold corresponding to the chosen friction coefficient η = 0.1 (i.e., ε c ≃ 0.0372 , see the table appearing in Figure 4).
Of course, it is convenient to reformulate the equations of motion in the form ṗ,q = V H H (0) − η p − Ω (0) , 0 where all the terms appearing in the expansion of the Hamiltonian H (0) are determined as in the discussion following formula (21). In particular, we have , , 0 = (0.0047704825942882482 , 0) .
Starting from these settings, we explicitly performed 60 steps of the normalization procedure described in subsection 5.1, by using the software package X̺óνoς, that is designed for making computer algebra, with a special care to its possible applications to Celestial Mechanics problems (see [34] for an introduction to its main concepts). Since none of the canonical transformations prescribed by the algorithm increases the polynomial degree, all the expansions of the Hamiltonian parts are rather compact, because they are at most quadratic in the actions as the initial H (0) . In Figure 5, we have reported the norms of the generating functions X (r) , ξ (r) and χ (r) 2 with the index r ∈ [1, 60] ; more precisely, we have calculated the sum of the absolute values of the coefficients appearing in (30), (29) and (42), respectively. Actually, the precision of the computation is affected by the truncation rules on the expansions, that we arranged so to neglect all the terms having a Fourier harmonic k with l 1 -norm |k| = 2 i=1 |k i | > 122 . In particular, this implies also that the contributions due to many relevant terms independent from the angles are not taken into account for r > 30 ; thus, the plot of ξ (r) has been stopped when the values corresponding to some indexes r have begun to be unrealistically small with respect to the previous ones. Let us emphasize that the geometrical decrease of the generating functions looks quite sharp in the semi-log scale of Figure 5; this behavior is in agreement with the analytical estimates on the algorithm constructing the Kolmogorov's normal form, when it is reformulated according a classical formal scheme (see [32]). Furthermore, also the l 1 -norm of the external forcing frequency vector Ω (r) is reported in Figure 5. In this case, the geometrical decrease is rather sharp until a "saturation threshold value", that is of order 10 −16 ; for r > 30 the value of Ω (r) is approximately constant. This unpleasant phenomenon can be easily explained, by taking into account that the numerical determination of the initial Ω (0) is affected by the unavoidable round-off errors. Such an uncertainity (due to the application of the frequency analysis numerical method) is propagated to all the sequence of the external forcing frequency vectors by the recursive definitions (27) and (37). Thus, the computed plot of r → Ω (r) agrees with the expectation that lim r→∞ Ω (r) = 0 . Finally, we can conclude that Figure 5 makes evident that the constructing procedure is converging to the Kolmogorov's normal form for dissipative systems, which is characterized by equations (22) and (23).
We now consider another test, checking the accuracy of the conjugacy canonical transformation K (r) which is provided after having performed the r-th normalization step according to the definition (48). Some previous works studying the construction of the Kolmogorov's normal form stressed that it can be used to integrate the equations of motion on an invariant torus characterized by a frequency vector ω . In fact, one can refer to the following ideal scheme (see, e.g., [46] and [31]): where Φ t ω·P is nothing but the flow induced by ω · P , that is the only effective part for the normalized equations of motion (22)-(23), when P = 0; moreover, let us recall that ψ (∞) = lim r→∞ K (r) is the conjugacy transformation whose the existence is ensured by a KAM-like statement (under suitable hypotheses). Of course, ψ (∞) can not be explicitly calculated, but this can be done for (a truncated expansion of) K (r) with a possibly large value of the index r . Thus, it is conveninet to limit ourselves to consider a numerical approximation t → p(t), q(t) of the motion law on an invariant torus characterized by a frequency vector ω , so to check if the following relation is satisfied: where K (r) i provides the approximately normalized i-th action for i = 1, . . . , n 1 . We check formula (51) with r = 60 and, again, in the special case of the dissipative forced pendulum, that is defined by the equations (6) and (4) with the values of the parameters reported in (49). In particular, we consider the motion law t → p(t), q(t) on the attracting invariant torus, related to the golden frequency ω * 1 = (3− √ 5)/2 . In Figure 6, we report the values of the approximately normalized action P = K (60) 1 p(2jπ), q(2jπ) as a function of its canonically conjugated angle Q ≃ Q(0) + 2jω * 1 π , when j = 0, 1, . . . , 10 000 . Let us stress that the initial value of Q(0) (and the corresponding P (0)) is unrelevant, when we are interested in checking the accuracy of the normalized canonical coordinates on all the invariant torus, because it is filled by the quasi-periodic orbit. However, the pair (P (0), Q(0)) is determined after having integrated the equations of motion for a relaxation time-span, according to the discussion reported at the end of section 2. Let us also remark that it is natural to completely neglect the second pair of coordinates; in fact, we recall that the dummy action p 2 does not affect the evolution of all other canonical variables; moreover, sinceq 2 =Q 2 = 1 , the second angle just describes the flowing of time. This explains why we have chosen to plot the approximately normalized first action when t = 2jπ , as in a standard Poincaré map. Figure 6 highlights that the action P (making part of the nearly normalized set of canonical variables) is close to zero for all the considered points generated by the Poincaré map of the flow on the attracting golden torus. Actually, the maximum of |P | is a few orders of magnitude bigger than the round-off errors threshold, in agreement with the expectation that their accumulation is unavoidable, while the very large number of computations required by the expansions are explicitly performed.  Let us emphasize that this positive check of the relation (51) does not concerns just with the algorithm constructing the Kolmogorov's normal form. Indeed, we had to select the good system of equations for which the golden torus was invariant and attracting, this required to determine a priori the numerical value of one of the parameters (namely, the external frequency Ω =Ω (5) reported in (49)), by using the frequency analysis. In our opinion, since the two so different techniques developed in the present paper provide results that are in agreement between them, this comparison strongly confirms our implementations of both methods.