## 2 Models and Algorithms
### Evolutionary Equations at Small Eccentricities  
Note that in many problems of celestial mechanics, eccentricities of body orbits are small (see [3; 5]). Let us consider this practically important case of small eccentricity \(e<<1\) and expand the right-hand sides of Equations (17)-(22) and (29)-(34) in Model 1 and 2, respectively, in power series in terms of \(e\). First, we can find an approximate solution to the Kepler Equation (9) and represent the eccentric anomaly \(E\) in the form of a converging series in \(e\) (see [5]), as follows:  
\[E=M+e\sin M+\frac{e^{2}}{2}\sin(2M)+\ldots.\]  
Using this solution, we obtain  
\[\cos E=\cos M-e\sin^{2}M-\frac{3e^{2}}{2}\cos M\sin^{2}M+\ldots,\] \[\sin E=\sin M+\frac{e}{2}\sin(2M)+\frac{e^{2}}{8}(3\sin(3M)-\sin M )+\ldots. \tag{35}\]  
Note that applying the computer algebra system Wolfram Mathematica (see [41]), one can calculate such expansions with any required accuracy, but here we restrict ourselves to computations up to the second order in \(e\). On substituting the expansions (35) into (17)-(21) and expanding their right-hand sides in powers of eccentricity \(e\), in the case of Model 1, we obtain  
\[\frac{da}{dt}=a\Big{(}1+2e\cos M+2e^{2}\cos(2M)\Big{)}\frac{\dot{ \gamma}(t)}{\gamma(t)},\] \[\frac{de}{dt}=\Big{(}\cos M+e\cos(2M)-\frac{9e^{2}}{2}\cos M\sin ^{2}M\Big{)}\frac{\dot{\gamma}(t)}{\gamma(t)},\] \[\frac{d\omega}{dt}=\Big{(}\frac{\sin M}{e}+\sin(2M)+\frac{e}{4}( 1+9\cos(2M))\sin M+\] \[\frac{e^{2}}{6}  
One can readily see that in the right-hand sides of Equation (36), there are some terms that do not depend on the mean anomaly \(M\), while most of other terms are proportional to the sine and cosine functions of \(M\) and its multiple. Such oscillating terms determine only the short-term oscillations of the orbital elements and do not influence their behavior with regard to long time intervals or secular perturbations. To eliminate the short-term oscillations of the orbital elements, it is sufficient to average the Equation (36) over the mean anomaly \(M\), which is an increasing function of time. We assume that the masses and the secular parts of the orbital elements change very slowly with time and the procedure of averaging does not change them.  
Recall that averaging of a function \(f(M)\) over the variable \(M\) and transition to the secular perturbations \(f^{(sec)}\) is reduced to calculating the integral  
\[f^{(sec)}=\frac{1}{2\pi}\int_{0}^{2\pi}f(M)dM.\]  
Substituting the right-hand sides of Equation (36) instead of the function \(f(M)\) and performing the integration, we obtain the following differential equations:  
\[\frac{da^{(sec)}}{dt}=a^{(sec)}\frac{\dot{\gamma}(t)}{\gamma(t)},\ \ \frac{dM^{(sec)}}{dt}=\frac{dE^{(sec)}}{dt}=\frac{\sqrt{\kappa}}{a^{3/2}\sqrt{ \gamma(t)}}, \tag{37}\]  
while the time derivatives of the secular perturbations of the rest of the orbital elements are zeros. Solving Equation (37), we obtain  
\[a^{(sec)}=a_{0}\gamma(t),\ \ M^{(sec)}=E^{(sec)}=\frac{\sqrt{\kappa}}{a_{0}^{ 3/2}}(\Phi(t)-\Phi(t_{0})),\ \ \Phi(t)=\int_{0}^{t}\frac{dt}{\gamma^{2}(t)}\, \tag{38}\]  
where \(a_{0}\) is an initial value of the semi-major axis \(a\) at the instant \(t=t_{0}\).