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

An analytical method for investigation of the evolution of dynamical systems {\it 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 the 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 the General Relativity [1], orbits and precession of planets of the Solar system [2,3], galactic dynamics and gravitating N-body problem [4], etc. Therefore, development of more advanced methods, especially of those enabling noniterative, and hence with negligible error storage at each iteration, for description of 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 the 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: mail@xxx.xx time error-free schemes developed even for particular models such as [6], can be efficient.
Methods of theory of dynamical systems are rather efficient for numerical studies of complex systems, including chaotic, 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 to characterise 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 at resolvent computation. Hence, the Hamiltonian and the initial function can be represented in more convenient variables depending on the particular dynamical system and aims, e.g. unbounded and non-periodic, differing from action-angle variables.
Hamiltonian formalism introduced for dissipative systems needs specific steps, and the transformation to evendimensional Liouville system usually mounts additional constraints, e.g. the dual Hamiltonian description [13,14] for arXiv:1611.08165v1 [astro-ph.IM] 24 Nov 2016 damping systems requires 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 the error-free 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 Hamiltonian are considered in [10,11] in a unified forṁ 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 Considering f (·) as a function of generalized coordinate, conjugate momenta and time, the equation describing the evolution of this function in the phase space changing according to (1), gets the form Then, according to 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 By means of the Laplace transform one can rewrite equation (5) 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 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 Per contra indeed, due to (9), (11) f (x, s) ≡ L e tL g(x) = R s (L)g(x) .
Finally, if A = L(β ) and B = L 0 by means of (16) -(17) we get the desired formula, f (x,t) = e tL g(x) = e tL 0 +β where i.e. at fixed perturbation β , by variation of N, we can reach a given accuracy. In order to reach higher accuracies, the exp(tL) in (19) should be calculated recursively substituting the equation itself under the integral.
Notably, for N = 1 f (x,t) = e t(L 0 +β L 1 ) g(x) Hence, the evolution for 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.
Presented method can be slightly modified (by adding an extra variable t withṫ = 1 equation) for the possible explicit time-dependent Hamiltonian systems.
The developed scheme, besides providing 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 analytical solution for 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 3 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 upon the convenience 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 duality principle [15], one obtains rescaled Hamiltonian and by "Dual" KAM theorem perturbed invariant tori stay preserved. Consider a dynamical system given by a two-dimensional Hamiltonian of two oscillators in action-angle form with different combination of a small perturbation V 1 (I, ϑ ) = I 1 cos ϑ 2 , V 2 (I, ϑ ) = I 1 sin ϑ 2 + I 2 cos ϑ 1 .
We have the following exact form for the operator e tL 0 g(I, ϑ ) = g(I, ϑ + tω).
(27)  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 proportional, at the accuracy 10 −6 .

Conclusions
We described a method of numerical investigation of evolution of dynamical systems independent on time accuracy, i.e. not based on iterations as the typical numerical codes but based on exact computation of the evolution function. In such formulation the method enables the use of computer algebra.
Essential point here is that the avoidance of the Laplace transform and spatial integration eliminates possible imaginary component appearance during computer algebra calculation, as could happen at 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 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 gravitating N-body problem, i.e. to the long term planetary system evolution, galaxy, galaxy cluster dynamics.
Another intriguing area of application can be the evolution of cosmological perturbations vs 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 for-mation and evolution, back reaction effects can be revealed in more details.