Evolution of perturbed dynamical systems: analytical computation with time independent accuracy

An analytical method for investigation of the evolution of dynamical systems with independent on time accuracy is developed for perturbed Hamiltonian systems. The error-free estimation using of computer algebra enables the application of the method to complex multi-dimensional Hamiltonian and dissipative systems. It also opens principal opportunities for the qualitative study of chaotic trajectories. The performance of the method is demonstrated on perturbed two-oscillator systems. It can be applied to various non-linear physical and astrophysical systems, e.g. to long-term planetary dynamics.


Introduction
The high accuracy investigation of the evolution of dynamical systems is a principal task in various astrophysical and physical problems, e.g. trajectories of artificial satellites testing General Relativity [1], orbits and precession of planets of the Solar System [2,3], galactic dynamics and gravitating N-body problem [4], etc. Therefore, the development of more advanced methods, especially of those enabling noniterative, and hence with negligible error storage at each iteration, for a description of the evolution of dynamical systems or broadening the classes of systems to which they are applicable, are of particular interest.
Among the principal issues is the numerical treatment of systems in chaotic regimes, including the tools for revealing chaos. For example, since, on the one hand, Lyapunov exponents are characteristics in the limit of infinite time, on the other hand, long iterative trajectory computations lead to inevitable accumulation of errors, suggested numerical procedures for n-dimensional non-linear systems (e.g. [5]) can be applied and interpreted with limited efficiency. Hence, a e-mail: gurzadyana@gmail.com time error-free schemes developed even for particular models such as [6] can be efficient.
The methods of the theory of dynamical systems are rather efficient for numerical studies of complex systems, including chaotic and random ones [7][8][9], and that approach is used below. We develop further the resolvent method suggested in [10,11]. The main issues elaborated below are: 1. Avoiding of the explicit use of Laplace and inverse Laplace transforms, and obtaining an essentially simplified evolution function representation, we broaden the classes of dynamical systems to which the method is applicable, namely, to non-linear Hamiltonian systems with many degrees of freedom. This is crucial, particularly, for N-body codes [4]. 2. The evolution function is obtained analytically i.e. with an accuracy independent of time and establishes the asymptotic and time-average characteristics of the system and, hence, enables one to deal also with chaotic trajectories. To this end, the error-free estimation of phase point (state) evolution in time allows one to characterize the behavior and system's accessory at ergodic or mixing phase space regions. This is, again, a principal topic for gravitating N-body systems known to possess chaotic properties [12]. 3. The choice of the initial function, i.e. of the one not affected by perturbations, as a periodic function is justified in various physical contexts; the newly applied procedure does not involve the Fourier expansion used previously in the resolvent computation. Hence, the Hamiltonian and the initial function can be represented by more convenient variables, depending on the particular dynamical system and aims, e.g. unbounded and non-periodic, differing from the action-angle variables.
The Hamiltonian formalism introduced for dissipative systems needs specific steps, and the transformation to even-dimensional Liouville system usually mounts additional constraints, e.g. the dual Hamiltonian description [13,14] for damping systems requires the introduction of auxiliary supplement. Based on this, the developed method allows one to follow the trajectory change of the Hamiltonian's principal term, not paying attention to the additional part. Then the reduction in [11] is again applicable via the scheme introduced here.
The method developed below, as we show also with computations on sample systems, enables one to perform an errorfree study for broad classes of dynamical systems-of perturbed Hamiltonian and dissipative ones-and can be applied to complex physical systems.

The evolution function and the resolvent
Hamiltonian equations for a perturbed system given by the Hamiltonian are considered in [10,11] in a unified form, where Systems (1) are the key ones in Kolmogorov-Arnold-Moser (KAM) theory stating the survival of the perturbed tori when certain conditions are fulfilled [7].
Then for any smooth function f (x) the following equation holds: where Considering f (·) as a function of generalized coordinates, conjugate momenta and time, the equation describing the evolution of this function in phase space changing according to (1) gets the form Then, according to the Hamiltonian equations, The evolution equation can be changed from the canonical variables to action-angle ones used in KAM-analysis to describe the phase portrait (Tor n × R n ) of the trajectories (see [7], §50). Then and where one can rewrite Eq. (1) as follows: where the resolvent R s (L) is equal to For any well-defined operators A and B one can show the resolvent satisfies the following relation: Then, by means of the inverse Laplace transform, the time domain of R s (·) can be defined, e.g. for A where the contour of integration is Thus, one can rewrite (13) as In contrast, indeed, due to (9) and (11) f Finally, if A = L(β) and B = L 0 by means of (16)-(17) we get the desired formula, where i.e. at fixed perturbation β, by variation of N , we can reach a given accuracy. In order to reach higher accuracies, ex p(t L) in (19) should be calculated recursively substituting the equation itself under the integral. Notably, for N = 1 Hence, the evolution for the initial function f (x, 0) = g(x) is found. Then, after selecting the point (state) in the phase space, e.g. the maximum of f ( · , t), for each t we pursue the time evolution of the point (state), thus determining the evolution of the Hamiltonian system in those time instances. The presented method can be slightly modified (by adding an extra variable t with the equationṫ = 1) for the possible explicit time-dependent Hamiltonian systems.
The developed scheme, besides providing a computation accuracy independent on time, also allows one to investigate the qualitative properties of the Hamiltonian system within the ergodic theory approach. Namely, despite the fact that an analytical solution for the Hamiltonian equations may not even exist due to internal resonances between different degrees of freedom, the precise time evolution for any state can be derived using the developed method. Such results can be used to distinguish to which class of dynamical systems the investigated system belongs to, e.g. of ergodic, weak mixing, K-system, etc.
The obtained exact form of f ( · , t) is more convenient for the study of asymptotic and average properties of the system. In a Hilbert space the time average of the evolution (20) of the system is obtained by von Neumann's theorem as

Error-free evolution: sample systems
The elimination of the explicit Laplace transform provides superior feasibility for computer algebra calculation than those described in [10,11]. Note that one is free to choose the initial function g(x) from the set of the smooth functions as seems convenient for the studied system or the aims. For the sample system, we will consider β 1 for simplicity. It is important to note that, in the case of sufficiently large perturbations, i.e. β → ∞, by applying the duality principle [15], one obtains a rescaled Hamiltonian and by the "dual" KAM theorem perturbed invariant tori are preserved. Consider a dynamical system given by a two-dimensional Hamiltonian of two oscillators in action-angle form with a different combination of a small perturbation, We have the following exact form for the operator: The initial function f (I, ϑ, 0)  Then the evolution of the systems by means of the described procedure is obtained, respectively, − ω 2 2β (I 1 − I 10 )I 2 (cos ϑ 1 − cos(ϑ 1 + tω 1 )) +(sin(ϑ 1 + tω 1 ) − sin ϑ 1 ) sin(ϑ 2 − ϑ 20 + tω 2 ) The evolution of the surface determined by the function f 1 (I, ϑ, t), i.e. corresponding to V 1 (I, ϑ) perturbation system, at different instances of time for I 1 = const, I 2 = const is shown in Fig. 1. The surface change and therefore the trajectories are regular in view of t = 0, 20, 50 and t = 10, 30 matching. In Figs. 2 and 3 the variations of (I, ϑ) by the time for different values of β are shown. One can observe deviations compared with the unperturbed (integrable) system's evolution. In Fig. 4 the contour lines are presented for the function f 2 (I, ϑ, t). The phase tori deformed at the increase of β in time are not destroyed, but remain only slightly deformed, as stated by KAM theory.
During computations, at the increase (decrease) of β by an order of magnitude, the calculation time varied with proportionality, at the accuracy 10 −6 .

Conclusions
We described a method of numerical investigation of the evolution of dynamical systems independent on the time accuracy, i.e. not based on iterations as the typical numerical codes but based on an exact computation of the evolution function. In such a formulation the method enables the use of computer algebra.
The essential point here is that the avoidance of the Laplace transform and spatial integration eliminates the possible appearance of an imaginary component during computer algebra calculation, as could happen with the previous procedure. Then the resulting form for the evolution function is far more practical in further study and analysis. We demonstrated the method on an example of the two-dimensional Hamiltonian of two oscillators and the error-free time evolution was shown for various values of the perturbation parameter β; the survival of the perturbed tori was shown in accordance to KAM theory.
The developed method can be applied to various physical and astrophysical many-dimensional non-linear systems, most notably, to the gravitating N-body problem, i.e. to the long-term planetary system evolution, galaxy, and galaxy cluster dynamics.
Another intriguing area of application can be the evolution of cosmological perturbations vs. the N-body formalism, particularly, for scalar perturbations within the discrete gravitating masses approach of [16] due to the possibility of adopting the modified method (dependent on the specific prospect) to dissipative systems (e.g. in quasistatic approximation or via scheme discussed in [13]). Then, after estimation of the error-free expressions, the cosmic structure formation and evolution, and back reaction effects can be revealed in more detail.
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .