ON THE SOLUTION OF REGULARIZED EQUATIONS IN PERTURBATION THEORY

M. S. Yarov-Yarovoi

Formulas are derived that determine small perturbations in coordinates and time, if the regularizing variable is used as the independent variable. The method is applicable to all types of unperturbed orbits, except rectilinear-parabolic and rectilinear-hyperbolic ones. The determination of perturbations in rectilinear motion will be considered in a separate article.

1. Spatial motion of a material point with rectangular coordinates \(x, y, z\) under the action of conservative forces having the force function

\[ U = U(x, y, z, t), \]

and of nonconservative forces with components

\[ X = X(x, y, z, t), \quad Y = Y(x, y, z, t), \quad Z = Z(x, y, z, t) \]

is described by the system of differential equations

\[ \ddot{x}=\frac{\partial U}{\partial x}+X,\quad \ddot{y}=\frac{\partial U}{\partial y}+Y,\quad \ddot{z}=\frac{\partial U}{\partial z}+Z. \tag{1.1} \]

In problems of celestial mechanics and astrodynamics the function \(U\) has the form

\[ U=\frac{\mu}{r}+R \quad (r=\sqrt{x^2+y^2+z^2}), \tag{1.2} \]

where \(R, X, Y, Z\) have small parameters as factors.

The solution of system (1.1) is, in the general case, rather cumbersome, since the right-hand sides have a complicated analytic structure. This structure, and at the same time the analytic properties of the solution, can be simplified by introducing new independent variables (regularization of the equations [1, 2]). In the present work, only one independent variable \(\tau\) is used for regularization.

If differentiation with respect to \(\tau\) is denoted by a prime, then the differential equation for the time \(t\) is chosen in the form

\[ t' = I(x, y, z, t). \tag{1.3} \]

Here the function \(I\) is as yet arbitrary.

Introducing into system (1.3), in place of \(t\), the new independent variable \(\tau\) by means of equation (1.3), we obtain

\[ t'=I,\quad x''=\frac{I'}{I}x' + I^2\left(\frac{\partial U}{\partial x}+X\right), \tag{1.4} \]

\[ y''=\frac{I'}{I}y' + I^2\left(\frac{\partial U}{\partial y}+Y\right),\qquad z''=\frac{I'}{I}z' + I^2\left(\frac{\partial U}{\partial z}+Z\right). \tag{1.4} \]

By choosing \(I\) in a suitable way, one can ensure that the right-hand sides of equations (1.4) have simpler analytic properties. For example, if, following Siegel [1], in problems of celestial mechanics and astrodynamics one uses an energetic approach to the problem, then as \(I\) one may choose

\[ I=\frac{1}{\dfrac{\mu}{r}+R^*},\qquad R^*=z_1R+z_2A, \tag{1.5} \]

where \(z_1, z_2\) are certain constants whose values are fixed in advance. In formula (1.5)

\[ A=\int (Xdx+Ydy+Zdz)=\int (Xx'+Yy'+Zz')\,d\tau \]

is the work of all nonconservative forces (in fact, a new variable).

1. At first glance one might think that the solution of the system of equations (1.4) is represented in the form of entire functions of \(\tau\), as was the case for the solution of the regularized equations of the two-body problem [2]. However, in the general case this is not so. For example, the regularized equations of the problem of two fixed centers no longer yield a solution in the form of entire functions of \(\tau\) [3].

Therefore we seek the solution of equations (1.4) in the form

\[ t=t_0+\delta t,\qquad x=x_0+\delta x,\qquad y=y_0+\delta y,\qquad z=z_0+\delta z, \tag{2.1} \]

where \(x_0, y_0, z_0, t_0\) are the coordinates and time in the unperturbed motion, satisfying the equations

\[ t_0'=\frac{r_0}{\mu},\qquad x_0''=\frac{r_0'}{r_0}x_0' - \frac{1}{\mu}\frac{x_0}{r_0}, \]

\[ y_0''=\frac{r_0'}{r_0}y_0' - \frac{1}{\mu}\frac{y_0}{r_0},\qquad z_0''=\frac{r_0'}{r_0}z_0' - \frac{1}{\mu}\frac{z_0}{r_0}. \tag{2.2} \]

If the perturbations \(\delta x, \delta y, \delta z, \delta t\) are sufficiently small, then \(r_0/\mu\) and \(I\) also differ little from one another. Then

\[ I=\frac{r}{\mu}-\frac{r}{\mu}\frac{rR^*}{\mu+rR^*},\qquad \frac{I'}{I}x'=\frac{r'}{r}x'-\frac{(rR^*)'}{\mu+rR^*}x' \tag{2.3} \]

(similarly for \(y, z\)),

\[ -I^2\frac{\mu x}{r^3} = -\frac{1}{\mu}\frac{x}{r} + \frac{2R^*+\dfrac{r}{\mu}R^{*2}}{(\mu+rR^*)^2}\,x \]

(similarly for \(y, z\)). Taking into account equations (1.4), (2.1), (2.2), (2.3), we obtain

\[ \delta t'=\frac{r}{\mu}-\frac{r_0}{\mu}+T^*, \]

\[ \tag{2.4} \]

\[ \delta x''=\frac{r'}{r}x' - \frac{r_0'}{r_0}x_0' - \frac{1}{\mu}\left(\frac{x}{r}-\frac{x_0}{r_0}\right)+X^* \]

(similarly for \(y, z\)).

Here we have

\[ T^*=-\frac{r}{\mu}-\frac{rR^*}{\mu+rR^*}, \]

\[ X^*=-\frac{(rR^*)'}{\mu+rR^*}\,x' +\frac{2R^*+\dfrac{r}{\mu}R^{*2}}{(\mu+rR^*)^2}\,x +\frac{r^2}{(\mu+rR^*)^2}\left(\frac{\partial R}{\partial x}+X\right). \tag{2.5} \]

(similarly for \(y,z\)).

We shall assume that \(R^*, R\) and its partial derivatives with respect to \(x,y,z\), as well as \(X,Y,Z\), are quantities of the same order of smallness, and we shall single out, in the right-hand sides of equations (2.4), the terms of first order with respect to \(\delta x,\delta y,\delta z\). We have

\[ r-r_0=\frac{x_0}{r_0}\,\delta x+\frac{y_0}{r_0}\,\delta y+\frac{z_0}{r_0}\,\delta z+\varepsilon(r), \tag{2.6} \]

\[ \frac{r'}{r}\,x'-\frac{r_0'}{r_0}\,x_0' =x_0'\left(\frac{x_0}{r_0^2}\,\delta x\right)' +x_0'\left(\frac{y_0}{r_0^2}\,\delta y\right)' + \]

\[ +\,x_0'\left(\frac{z_0}{r_0^2}\,\delta z\right)' +\frac{r_0'}{r_0}\,\delta x' +\varepsilon\left(\frac{r'}{r}x'\right) \tag{2.7} \]

(similarly for \(y,z\)),

\[ \frac{x}{r}-\frac{x_0}{r_0} =\frac{r_0^2-x_0^2}{r_0^3}\,\delta x -\frac{x_0y_0}{r_0^3}\,\delta y -\frac{x_0z_0}{r_0^3}\,\delta z +\varepsilon\left(\frac{x}{r}\right) \tag{2.8} \]

(similarly for \(y,z\)).

Here

\[ \varepsilon(r)=r_0\left\{ \left[ 1+2\left( \frac{x_0}{r_0}\frac{\delta x}{r_0} +\frac{y_0}{r_0}\frac{\delta y}{r_0} +\frac{z_0}{r_0}\frac{\delta z}{r_0} \right) +\left(\frac{\delta x}{r_0}\right)^2 +\right. \]

\[ \left. +\left(\frac{\delta y}{r_0}\right)^2 +\left(\frac{\delta z}{r_0}\right)^2 \right]^{1/2} -1 -\frac{x_0}{r_0}\frac{\delta x}{r_0} -\frac{y_0}{r_0}\frac{\delta y}{r_0} -\frac{z_0}{r_0}\frac{\delta z}{r_0} \right\}, \]

\[ \varepsilon\left(\frac{r'}{r}x'\right) = \frac{1}{2}(x_0+\delta x)' \left\{ \ln\left[ 1+2\left( \frac{x_0}{r_0}\frac{\delta x}{r_0} +\right. \right. \]

\[ \left. \left. +\frac{y_0}{r_0}\frac{\delta y}{r_0} +\frac{z_0}{r_0}\frac{\delta z}{r_0} \right) +\left(\frac{\delta x}{r_0}\right)^2 +\left(\frac{\delta y}{r_0}\right)^2 +\left(\frac{\delta z}{r_0}\right)^2 \right]\right\}' - \tag{2.9} \]

\[ -\,x_0'\left(\frac{x_0}{r_0^2}\,\delta x\right)' -x_0'\left(\frac{y_0}{r_0^2}\,\delta y\right)' -x_0'\left(\frac{z_0}{r_0^2}\,\delta z\right)', \]

\[ \varepsilon\left(\frac{x}{r}\right) = -\frac{1}{r_0}(x_0+\delta x) \left[ 1+2\left( \frac{x_0}{r_0}\frac{\delta x}{r_0} +\frac{y_0}{r_0}\frac{\delta y}{r_0} + \right. \]

\[ +\frac{z_0}{r_0}\frac{\delta z}{r_0}\right) +\left(\frac{\delta x}{r_0}\right)^2 +\left(\frac{\delta y}{r_0}\right)^2 +\left(\frac{\delta z}{r_0}\right)^2 \right]^{-1/2} - \]

\[ -\frac{1}{r_0}(x_0+\delta x) +\frac{x_0^2}{r_0^3}\,\delta x +\frac{x_0y_0}{r_0^3}\,\delta y +\frac{x_0z_0}{r_0^3}\,\delta z \]

(similarly for \(y,z\)) is the set of terms of second and higher orders with respect to \(\delta x,\delta y,\delta z\).

In the functions \(T^*, X^*, Y^*, Z^*\) one can also separate out the terms of first order of smallness \(T_0^*, X_0^*, Y_0^*, Z_0^*\):

\[ T^*=T_0^*+\varepsilon(T^*),\qquad X^*=X_0^*+\varepsilon(X^*) \tag{2.10} \]

(similarly for \(y,z\)).

Substituting expressions (2.7), (2.8), (2.10) into equations (2.5), we obtain

\[ \delta t' = \frac{1}{\mu}\frac{x_0}{r_0}\,\delta x +\frac{1}{\mu}\frac{y_0}{r_0}\,\delta y +\frac{1}{\mu}\frac{z_0}{r_0}\,\delta z +T_0^*+\varepsilon_t, \]

\[ \begin{aligned} \delta x''={}& \left( -\frac{2r_0'x_0'}{r_0^3}x_0 +\frac{x_0'}{r_0^2}x_0' +\frac{1}{\mu}\frac{x_0}{r_0^3}x_0 -\frac{1}{\mu}\frac{1}{r_0} \right)\delta x \\ &+ \left( -\frac{2r_0'x_0'}{r_0^3}y_0 +\frac{x_0'}{r_0^2}y_0' +\frac{1}{\mu}\frac{x_0}{r_0^3}y_0 \right)\delta y \\ &+ \left( -\frac{2r_0'x_0'}{r_0^3}z_0 +\frac{x_0'}{r_0^2}z_0' +\frac{1}{\mu}\frac{x_0}{r_0^3}z_0 \right)\delta z \\ &+ \left( \frac{r_0'}{r_0} +\frac{x_0'}{r_0^2}x_0 \right)\delta x' +\frac{x_0'}{r_0^2}y_0\,\delta y' +\frac{x_0'}{r_0^2}z_0\,\delta z' +X_0^*+\varepsilon_x \end{aligned} \tag{2.11} \]

(similarly for \(y,z\)).

Here

\[ \varepsilon_t=\frac{1}{\mu}\varepsilon_r+\varepsilon(T^*),\qquad \varepsilon_x= \varepsilon\left(\frac{r'}{r}x'\right) -\frac{1}{\mu}\varepsilon\left(\frac{x}{r}\right) +\varepsilon(X^*) \tag{2.12} \]

(similarly for \(y,z\))—the sets of terms of second and higher orders of smallness.

1. The system of equations (2.11) is a system of linear nonhomogeneous differential equations with variable coefficients, whose right-hand sides (in \(\varepsilon_t,\varepsilon_x,\varepsilon_y,\varepsilon_z\)) contain perturbations of second and higher orders with respect to \(\delta t,\delta x,\delta y,\delta z\). By means of a known device they can be reduced to integral form as follows. Let

\[ q_1=\delta x,\quad q_2=\delta y,\quad q_3=\delta z,\quad q_4=\delta x', \]

\[ q_5=\delta y',\quad q_6=\delta z',\quad q_7=\delta t. \tag{3.1} \]

Then the system of differential equations (2.11) will take the form

\[ q_i'-q_{i+3}=0\qquad (i=1,2,3), \]

\[ q'_j+\sum_{k=1}^{6} a_{jk}q_k=V_j \quad (j=4,5,6), \tag{3.2} \]

\[ q'_7+a_{71}q_1+a_{72}q_2+a_{73}q_3=V_7. \]

The general solution of the corresponding homogeneous system will have the form

\[ q_i=\sum_{j=1}^{7} K_j q_i^{(j)}. \tag{3.3} \]

We shall find the explicit expressions for the functions \(q_i^{(j)}\) in the next subsection. To solve equations (3.2), let us apply the method of variation of arbitrary constants [4]. The matrix of the system of fundamental solutions has the form

\[ \|D\|= \left\| \begin{array}{ccccc} q_1^{(1)} & q_1^{(2)} & \cdot & \cdot & q_1^{(7)}\\ q_2^{(1)} & q_2^{(2)} & \cdot & \cdot & q_2^{(7)}\\ \cdot & \cdot & \cdot & \cdot & \cdot\\ q_7^{(1)} & q_7^{(2)} & \cdot & \cdot & q_7^{(7)} \end{array} \right\|. \tag{3.4} \]

We denote by \(D_{ik}\) the algebraic complement of the element of this matrix lying in the \(i\)-th row and \(k\)-th column. Then

\[ K_i=\int_{\tau}\frac{1}{D}\left(D_{4i}V_4+D_{5i}V_5+D_{6i}V_6+D_{7i}V_7\right)d\tau \quad (i=1,\ldots,7), \tag{3.5} \]

since the right-hand sides of the first three equations are zeros.

The general solution of the nonhomogeneous system is also represented by formulas (3.3), but in these formulas \(K_i\) consist of two terms, the first of which is arbitrary, while the second satisfies equation (3.5).

1. Let us now find the fundamental system for equations (3.2) by the well-known Poincaré theorem. For this purpose, we differentiate partially with respect to the initial values of the coordinates constituting the regularized velocity and the time the general solution of the regularized equations of the two-body problem (2.2) ([2], equations (6.11)). After such differentiation, the initial values may be set equal to the pericentric (or apocentric, if the initial velocity is less than circular) values. As is known, in this case the initial direction of the velocity is perpendicular to the radius vector. Such a restriction on the initial data is necessary in order to simplify as much as possible the rather complicated calculation of the algebraic complements of the determinant \(D\). The formulas obtained will be suitable for all types of orbits except rectilinear-parabolic and rectilinear-hyperbolic ones, since only on these orbits there is no point at which the velocity would have a zero value, i.e., would be perpendicular to the radius vector. Analogous formulas in the case of all types of orbits of rectilinear motion (when the initial value of the velocity is parallel to the radius vector) will be obtained in a special paper.

Under the restrictions made above, at \(\tau=0\),

\[ x=\rho_0,\quad y=0,\quad z=0,\quad x'=0,\quad y'=V'_0=V_0\frac{\rho_0}{\mu},\quad z'=0,\quad t=t_{00}. \]

Therefore

\[ q_1^{(1)}=1+\frac{\mu}{\rho_0^2}\rho_0''\frac{\partial A}{\partial h} =(1+2\alpha^{-1}+2\alpha^{-2})+ \]

\[ +(-2\alpha^{-1}-2\alpha^{-2})\operatorname{ch}H +(\alpha^{-1}+\alpha^{-2})H\operatorname{sh}H, \]

\[ q_1^{(2)}=q_1^{(3)}=0, \]

\[ q_1^{(4)}=B=\sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2}\operatorname{sh}H, \]

\[ q_1^{(5)}=-\frac{V_0}{\rho_0}\frac{\partial A}{\partial h} =\sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,\alpha^{-2}(-2+2\operatorname{ch}H-H\operatorname{sh}H), \]

\[ q_1^{(6)}=q_1^{(7)}=0, \]

\[ q_2^{(1)}=-\frac{\mu V_0}{\rho_0}\rho_0''\frac{\partial A}{\partial h} =\sqrt{2+\alpha}\,\alpha^{-1/2}\bigl[(1+\alpha^{-1})\operatorname{sh}H+ \]

\[ +(-1-\alpha^{-1})H\operatorname{ch}H\bigr], \]

\[ q_2^{(2)}=1+\frac{\rho_0''}{\rho_0}A =-\alpha^{-1}+(1+\alpha^{-1})\operatorname{ch}H, \]

\[ q_2^{(3)}=0,\qquad q_2^{(4)}=\frac{V_0}{\mu}A =\sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,\alpha^{-1}(-1+\operatorname{ch}H), \]

\[ \tag{4.1} \]

\[ q_2^{(5)}=B-2h\frac{\partial B}{\partial h} +2\frac{\mu^2}{\rho_0^2}\rho_0''\frac{\partial B}{\partial h} =\sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2}\bigl[-2\alpha^{-1}\operatorname{sh}H+ \]

\[ +(1+2\alpha^{-1})H\operatorname{ch}H\bigr], \]

\[ q_2^{(6)}=q_2^{(7)}=q_3^{(1)}=q_3^{(2)}=0, \]

\[ q_3^{(3)}=1-\frac{1}{\mu\rho_0}A =(1+\alpha^{-1})-\alpha^{-1}\operatorname{ch}H, \]

\[ q_3^{(4)}=q_3^{(5)}=0,\qquad q_3^{(6)}=B=\sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2}\operatorname{sh}H,\qquad q_3^{(7)}=0, \]

\[ q_4^{(1)}=\frac{\mu\rho_0''}{\rho_0^2}\frac{\partial B}{\partial h} =-\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\alpha^{-1/2}\bigl[(-1-\alpha^{-1})\operatorname{sh}H+ \]

\[ +(1+\alpha^{-1})H\operatorname{ch}H\bigr], \]

\[ q_4^{(2)}=q_4^{(3)}=0,\qquad q_4^{(4)}=1+\frac{2h}{\mu^2}A=\operatorname{ch}H, \]

\[ q_4^{(5)}=-\frac{V_0}{\rho_0}\frac{\partial B}{\partial h} =\sqrt{2+\alpha}\,\alpha^{-1/2}\alpha^{-1}(\operatorname{sh}H-H\operatorname{ch}H), \]

\[ q_4^{(6)}=q_4^{(7)}=0,\qquad q_5^{(1)}=-\frac{V_0\rho_0''}{\mu\rho_0}\left(2A+2h\frac{\partial A}{\partial h}\right)= \]

\[ =-\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\sqrt{2+\alpha}(-1-\alpha^{-1})H\operatorname{sh}H, \]

\[ q_5^{(2)}=\frac{\rho_0''}{\rho_0}B = \frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\alpha^{-1/2}(1+\alpha)\operatorname{sh}H, \qquad q_5^{(3)}=0, \]

\[ q_5^{(4)}=\frac{V_0}{\mu}B =\sqrt{2+\alpha}\,\alpha^{-1/2}\operatorname{sh}H, \]

\[ q_5^{(5)}=1-\frac{2h}{\mu^2}A-\frac{(2h)^2}{\mu^2}\frac{\partial A}{\partial h} +\frac{2\rho_0''}{\rho^2}\left(2A+2h\frac{\partial A}{\partial h}\right) = \]

\[ =\operatorname{ch}H+(1+2\alpha^{-1})H\operatorname{sh}H, \]

\[ q_5^{(6)}=q_5^{(7)}=q_6^{(1)}=q_6^{(2)}=0, \]

\[ q_6^{(3)}=-\frac{1}{\mu\rho_0}B =-\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\alpha^{-1/2}\operatorname{sh}H, \]

\[ q_6^{(4)}=q_6^{(5)}=0,\qquad q_6^{(6)}=1+\frac{2h}{\mu^2}A=\operatorname{ch}H,\qquad q_6^{(7)}=0, \tag{4.1} \]

\[ q_7^{(1)}=\frac{1}{\mu}\tau+\frac{2h}{\mu^3}C-\frac{2\rho_0''}{\mu\rho_0}C -\frac{\mu(\rho_0'')^2}{\rho_0^2}\frac{\partial C}{\partial h} = \]

\[ =\frac{\sqrt{\rho_0}}{\sqrt{\mu}}\alpha^{-1/2} \left[(2+4\alpha^{-1}+3\alpha^{-2})\operatorname{sh}H +(-2\alpha^{-1}-2\alpha^{-2})H+\right. \]

\[ \left.+(-1-2\alpha^{-1}-\alpha^{-2})H\operatorname{ch}H\right], \]

\[ q_7^{(2)}=\frac{V_0}{\mu}A =\frac{\sqrt{\rho_0}}{\sqrt{\mu}}\sqrt{2+\alpha}\,\alpha^{-1}(-1+\operatorname{ch}H), \qquad q_7^{(3)}=0, \]

\[ q_7^{(4)}=\frac{1}{\mu}A =\rho_0\alpha^{-1}(-1+\operatorname{ch}H), \]

\[ q_7^{(5)}=\frac{2V_0}{\mu^2}C+\frac{V_0\rho_0''}{\rho_0}\frac{\partial C}{\partial h} =\rho_0\sqrt{2+\alpha}\,\alpha^{-1/2} \left[(-\alpha^{-1}-\right. \]

\[ \left.-3\alpha^{-2})\operatorname{sh}H +2\alpha^{-2}H +(\alpha^{-1}+\alpha^{-2})H\operatorname{ch}H\right], \]

\[ q_7^{(6)}=0,\qquad q_7^{(7)}=1. \]

Here \(V_0\) is the initial value of the velocity, \(\rho_0\) is the initial value of the radius vector; \(h\) is the constant of the energy integral in the unperturbed motion

\[ h=\frac{1}{2}V_0^2-\frac{\mu}{\rho_0},\qquad \alpha=\frac{\rho_0\,2h}{\mu},\qquad \rho_0''=\frac{2h}{\mu^2}\rho_0+\frac{1}{\mu} =\frac{1}{\mu}(1+\alpha). \tag{4.2} \]

It is not difficult to see that the fundamental system found is normal. The functions \(A\), \(B\), \(C\) were introduced by us in [2]. They are entire functions of \(\tau\). Obviously, \(\partial A/\partial h\), \(\partial B/\partial h\), \(\partial C/\partial h\) are also entire functions of \(\tau\). From formulas (4.1), (4.2) the theorems follow.

Theorem 4.1. The general solution of the variational equations for the regularized equations of the two-body problem is represented by entire functions of the regularizing variable \(\tau\), provided only that the initial value of the radius vector \(\rho_0\) is different from zero.

Theorem 4.2. The elements of the fourth, fifth, and sixth rows of the fundamental matrix are the derivatives with respect to \(\tau\) of the corresponding elements of the first three rows.

Theorem 4.2 is convenient for checking the computations.

In formulas (4.1) there are also given expressions for particular solutions for hyperbolic motion. In this case

\[ H=\frac{\sqrt{2h}}{\mu}\,\tau . \tag{4.3} \]

Remark 4.1. In order to obtain from expressions (4.1) the formulas for elliptic motion, one must use the formulas

\[ H=-iE,\qquad \alpha=e-1,\qquad \rho_0=a(1-e), \tag{4.4} \]

\[ \alpha^{-1/2}H=\frac{E}{\sqrt{1-e}},\qquad \alpha^{-1/2}\operatorname{sh}H=\frac{\sin E}{\sqrt{1-e}}, \]

\[ \alpha^{-1/2}H\operatorname{ch}H=\frac{E\cos E}{\sqrt{1-e}}, \tag{4.5} \]

\[ \operatorname{ch}H=\cos E,\qquad H\operatorname{sh}H=-E\sin E. \]

The formulas for elliptic rectilinear motion are obtained from this for \(e=-1\).

To obtain the formulas for parabolic motion, one must introduce the dimensionless variable [2]

\[ \beta=\frac{\tau}{\sqrt{\mu}\sqrt{\rho_0}} . \tag{4.6} \]

Since

\[ \beta=\alpha^{-1/2}H, \]

in formulas (4.1) it is sufficient to expand the hyperbolic functions of \(H\) in a series in powers of \(\alpha,\beta\). The negative powers of \(\alpha\) then vanish. It follows from Theorem 4.1 that the functions \(q_i^{(j)}\) will be integral with respect to \(\beta\). The coefficients of the expansions of \(q_i^{(j)}\) in powers of \(\beta\) will be polynomials in \(\alpha\). The formulas for parabolic motion are obtained from these expansions for \(\alpha=0\).

The fundamental matrix \(\|D\|\) has, in the parabolic case, the form

\[ \begin{array}{ccc} 1+\dfrac{1}{12}\beta^4 & 0 & 0 \\[1.2em] \sqrt{2}\left(-\dfrac{1}{3}\beta^3\right) & 1+\dfrac{1}{2}\beta^2 & 0 \\[1.2em] 0 & 0 & 1-\dfrac{1}{2}\beta^2 \\[1.2em] \dfrac{1}{\sqrt{\mu}\sqrt{\rho_0}}\,\dfrac{1}{3}\beta^3 & 0 & 0\to \\[1.2em] \dfrac{1}{\sqrt{\mu}\sqrt{\rho_0}}\sqrt{2}\,(-\beta^2) & \dfrac{1}{\sqrt{\mu}\sqrt{\rho_0}}\,\beta & 0 \\[1.2em] 0 & 0 & \dfrac{1}{\sqrt{\mu}\sqrt{\rho_0}}\,(-\beta) \\[1.2em] \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}\left(\beta-\dfrac{1}{3}\beta^3-\dfrac{1}{60}\beta^5\right) & \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}\sqrt{2}\,\dfrac{1}{2}\beta^2 & 0 \end{array} \]

\[ \begin{array}{cccc} \sqrt{\mu}\sqrt{\rho_0}\,\beta & \sqrt{\mu}\sqrt{\rho_0}\sqrt{2}\left(-\dfrac{1}{12}\beta^4\right) & 0 & 0 \\[1.2ex] \sqrt{\mu}\sqrt{\rho_0}\sqrt{2}\,\dfrac{1}{2}\beta^2 & \sqrt{\mu}\sqrt{\rho_0}\left(\beta+\dfrac{2}{3}\beta^3\right) & 0 & 0 \\[1.2ex] 0 & 0 & \sqrt{\mu}\sqrt{\rho_0}\,\beta & 0 \\[1.2ex] \leftarrow\quad 1 & \sqrt{2}\left(-\dfrac{1}{3}\beta^3\right) & 0 & 0 \\[1.2ex] \sqrt{2}\,\beta & 1+2\beta^2 & 0 & 0 \\[1.2ex] 0 & 0 & 1 & 0 \\[1.2ex] \rho_0\,\dfrac{1}{2}\beta^2 & \rho_0\sqrt{2}\left(\dfrac{1}{3}\beta^3+\dfrac{1}{60}\beta^5\right) & 0 & 1 \end{array} \tag{4.7} \]

1. There is no need to discuss in detail how difficult it is to invert the fundamental matrix (3.4), having the concrete form (4.1) or (4.7). In the general case, when the direction of the axes \(x, y, z\) is chosen arbitrarily, the inversion reduces to the computation of determinants of the sixth order with very complicated elements. If the unperturbed orbit lies in the plane \(x, y\), then at most determinants of the fourth order have to be computed. However, even such a computation in fact proves exceptionally cumbersome because of the complexity of the elements of these determinants. In order to simplify these calculations still further, the restriction on the initial conditions of the problem mentioned above is introduced, i.e. the inversion of the matrix (3.4) is carried out under the assumption that the initial velocity is either equal to zero or perpendicular to the initial direction of the radius vector. As has already been said, the case of parallelism of these directions (rectilinear motion in the general case) will be considered in a special paper. The intermediate case, by a choice of the initial value of \(\tau\), can be reduced to the first case, i.e. to the case of perpendicularity of the directions of the initial velocity and the initial radius vector.

The computation of the algebraic complements even under the indicated restrictions involves carrying out an enormous number of calculations. Although, of the 28 algebraic complements necessary for obtaining formulas (3.5), only 13 are nonzero, nevertheless all the computational work together with checking took many weeks. The computations were carried out in parallel and independently for two cases: the general and the parabolic. The final formulas in both cases were compared using the remark (4.1). A general check on all the computations is provided by Liouville’s theorem, which determines the dependence on \(\tau\) of the determinant \(D\) of the fundamental matrix. Indeed, for the system (3.2), in view of \(a_{11}=a_{22}=a_{33}=a_{77}=0\),

\[ \sum_{i=1}^{7} a_{ii} = a_{44}+a_{55}+a_{66} = -4\,\dfrac{r_0'}{r_0} = -4(\ln r_0)'. \]

Hence

\[ D = D^{(0)} \left(\dfrac{r_0}{\rho_0}\right)^4 = \left(\dfrac{r_0}{\rho_0}\right)^4, \]

since the initial value \(D = D^{(0)}\) for the normal fundamental system is always equal to unity.

As a result of the calculations it turned out that all algebraic complements are divisible by \(r_0\). Therefore the nonzero algebraic complements will have the form

\[ \frac{D_{41}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2} \left[ 2\alpha^{-1}\operatorname{sh} H + \left(-\frac{1}{2}-\alpha^{-1}\right)\operatorname{sh} 2H \right], \]

\[ \frac{D_{51}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,\alpha^{-2} \left[ -\frac{3}{2} + 2\operatorname{ch} H - \frac{1}{2}\operatorname{ch} 2H \right], \]

\[ \frac{D_{42}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,\alpha^{-1} \left[ \left(\frac{1}{2}+2\alpha^{-1}\right) + \right. \]

\[ \left. + (-1-2\alpha^{-1})\operatorname{ch} H + \frac{1}{2}\operatorname{ch} 2H + \alpha^{-1}H\operatorname{sh} H \right], \]

\[ \frac{D_{52}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2} \left[ \alpha^{-1}\operatorname{sh} H + \right. \]

\[ \left. + \frac{1}{2}\alpha^{-1}\operatorname{sh} 2H + (-1-2\alpha^{-1})H\operatorname{ch} H \right], \]

\[ \frac{D_{63}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2} \left[ \alpha^{-1}\operatorname{sh} H + \left(-\frac{1}{2}-\frac{1}{2}\alpha^{-1}\right)\operatorname{sh} 2H \right] = \]

\[ = \frac{\rho_0}{r_0}\sqrt{\mu}\sqrt{\rho_0}\alpha^{-1/2}(-\operatorname{sh} H), \]

\[ \frac{D_{44}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \left[ (-2\alpha^{-1}-2\alpha^{-2}) + \right. \]

\[ \left. + (1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch} H - \alpha^{-2}H\operatorname{sh} H \right], \tag{5.1} \]

\[ \frac{D_{54}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{2+\alpha}\,\alpha^{-3/2} [-\operatorname{sh} H+H\operatorname{ch} H], \]

\[ \frac{D_{45}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{2+\alpha}\,\alpha^{-1/2} \left[ \alpha^{-1}\operatorname{sh} H + \left(-\frac{1}{2}-\frac{1}{2}\alpha^{-1}\right)\operatorname{sh} 2H \right] = \]

\[ = \frac{\rho_0}{r_0}\sqrt{2+\alpha}\,\alpha^{-1/2}(-\operatorname{sh} H), \]

\[ \frac{D_{55}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \left[ \left(-\frac{3}{2}\alpha^{-1}-\frac{3}{2}\alpha^{-2}\right) + (1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch} H + \right. \]

\[ \left. + \left(-\frac{1}{2}\alpha^{-1}-\frac{1}{2}\alpha^{-2}\right)\operatorname{ch} 2H \right] = \frac{\rho_0}{r_0} \left[ (1+\alpha^{-1})-\alpha^{-1}\operatorname{ch} H \right], \]

\[ \frac{D_{66}}{D}=\frac{D_{55}}{D}, \]

\[ \begin{gathered} \frac{D_{47}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \rho_0 \alpha^{-1} \Bigg[ \left(\frac{1}{2}+\frac{7}{2}\alpha^{-1}+\frac{9}{2}\alpha^{-2}\right) \\ {}+\left(-1-4\alpha^{-1}-4\alpha^{-2}\right)\operatorname{ch} H +\left(\frac{1}{2}+\frac{1}{2}\alpha^{-1} \right. \\ \left. {}-\frac{1}{2}\alpha^{-2}\right)\operatorname{ch} 2H +\left(\alpha^{-1}+3\alpha^{-2}\right)H\operatorname{sh} H \Bigg], \tag{5.1} \\[1em] \frac{D_{57}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \rho_0 \sqrt{2+\alpha}\,\alpha^{-3/2} \Big[ \alpha^{-1}\operatorname{sh} H \\ {}+\left(\frac{1}{2}+\alpha^{-1}\right)\operatorname{sh} 2H +\left(-1-3\alpha^{-1}\right)H\operatorname{ch} H \Big], \\[1em] \frac{D_{77}}{D}=1. \end{gathered} \]

Formulas for unperturbed parabolic motion were obtained independently.

\[ \begin{gathered} \frac{D_{41}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\left(-\beta-\beta^3\right), \\[0.7em] \frac{D_{51}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\sqrt{2}\left(-\frac{1}{4}\beta^4\right), \\[0.7em] \frac{D_{42}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\sqrt{2} \left(\frac{1}{2}\beta^2+\frac{1}{12}\beta^4\right), \\[0.7em] \frac{D_{52}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0} \left(-\beta-\frac{1}{6}\beta^3\right), \\[0.7em] \frac{D_{63}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0} \left(-\beta-\frac{1}{2}\beta^3\right) = \frac{\rho_0}{r_0}\sqrt{\mu}\sqrt{\rho_0}(-\beta), \tag{5.2} \\[0.7em] \frac{D_{44}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \left(1+\beta^2-\frac{1}{12}\beta^4\right), \qquad \frac{D_{54}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{2}\,\frac{1}{3}\beta^3, \\[0.7em] \frac{D_{45}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \sqrt{2}\left(-\beta-\frac{1}{2}\beta^3\right) = \frac{\rho_0}{r_0}\sqrt{2}(-\beta), \\[0.7em] \frac{D_{55}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \left(1-\frac{1}{4}\beta^4\right) = \frac{\rho_0}{r_0}\left(1-\frac{1}{2}\beta^2\right), \\[0.7em] \frac{D_{66}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \left(1-\frac{1}{4}\beta^4\right) = \frac{\rho_0}{r_0}\left(1-\frac{1}{2}\beta^2\right), \\[0.7em] \frac{D_{47}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \rho_0 \left(\frac{1}{2}\beta^2+\frac{1}{3}\beta^4-\frac{1}{40}\beta^6\right), \\[0.7em] \frac{D_{57}}{D} = \left(\frac{\rho_0}{r_0}\right)^2 \rho_0\sqrt{2} \left(\frac{1}{6}\beta^3+\frac{3}{20}\beta^5\right), \qquad \frac{D_{77}}{D}=1. \end{gathered} \]

To verify these expressions, we used remark (4.1) and the rule: the scalar product of the \(i\)-th row of the matrix \(\|D\|\) by the elements of the \(j\)-th row of the matrix \(\|D\|^{-1}\) is equal to 1 for \(i=j\) and to 0 for \(i\ne j\).

6. Theorem 6.1. If the equations of the perturbed motion (2.11) are brought to the form (3.2) and the method of variation of the arbitrary constants \(K_i\) of the general solution of the corresponding homogeneous system (3.3) with fundamental system (4.1) or (4.7) is applied to them, then the equations for \(K_i\) will have the form

\[ K_1=\int\left(\frac{\rho_0}{r_0}\right)^2 \mu\rho_0\left[E_{1x}(X^{*}+\varepsilon_x^{*})+ E_{1y}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ K_2=\int\left(\frac{\rho_0}{r_0}\right)^2 \mu\rho_0\left[E_{2x}(X^{*}+\varepsilon_x^{*})+ E_{2y}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ K_3=\int\left(\frac{\rho_0}{r_0}\right)^2 \mu\rho_0 E_{3z}(Z^{*}+\varepsilon_z^{*})\,d\beta, \]

\[ K_4=\int\left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\left[E_{4x}(X^{*}+\varepsilon_x^{*})+ E_{4y}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ \tag{6.1} K_5=\int\left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\left[E_{5x}(X^{*}+\varepsilon_x^{*})+ E_{5y}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ K_6=\int\left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\sqrt{\rho_0}\,E_{6z}(Z^{*}+\varepsilon_z^{*})\,d\beta, \]

\[ K_7=\int\left\{\left(\frac{\rho_0}{r_0}\right)^2 \sqrt{\mu}\rho_0\sqrt{\rho_0} \left[E_{7x}(X^{*}+\varepsilon_x^{*})+ E_{7y}(Y^{*}+\varepsilon_y^{*})\right]+ \right. \]

\[ \left. +\sqrt{\mu}\sqrt{\rho_0}(T+\varepsilon_t^{*})\right\}\,d\beta. \]

Here

\[ \varepsilon_t^{*}=\frac{1}{\mu}\varepsilon_r,\qquad \varepsilon_x^{*}=\varepsilon\left(\frac{r'}{r}-x'\right)-\frac{1}{\mu}\varepsilon\left(\frac{x}{r}\right) \]

(and analogously for \(y,z\)).

In equations (6.1), the functions \(E\) with subscripts are dimensionless factors in \(D_{ik}/D\). In the general case

\[ E_{1x}=\alpha^{-1/2}\left[2\alpha^{-1}\operatorname{sh}H+ \left(-\frac{1}{2}-\alpha^{-1}\right)\operatorname{sh}2H\right], \]

\[ E_{1y}=\sqrt{2+\alpha}\,\alpha^{-2}\left[-\frac{3}{2} +2\operatorname{ch}H-\frac{1}{2}\operatorname{ch}2H\right], \]

\[ \tag{6.2} E_{2x}=\sqrt{2+\alpha}\,\alpha^{-1}\left[ \left(\frac{1}{2}+2\alpha^{-1}\right) +\left(-1-2\alpha^{-1}\right)\operatorname{ch}H+ \frac{1}{2}\operatorname{ch}2H+\alpha^{-1}H\operatorname{sh}H \right], \]

\[ E_{2y}=\alpha^{-1/2}\left[ \alpha^{-1}\operatorname{sh}H+\frac{1}{2}\alpha^{-1}\operatorname{sh}2H+ \left(-1-2\alpha^{-1}\right)H\operatorname{ch}H \right], \]

\[ \begin{gathered} E_{3z}=\alpha^{-1/2}\left[\alpha^{-1}\operatorname{sh}H+ \left(-\frac12-\frac12\alpha^{-1}\right)\operatorname{sh}2H\right],\\ E_{4x}=(-2\alpha^{-1}-2\alpha^{-2})+(1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H-\alpha^{-2}H\operatorname{sh}H,\\ E_{4y}=\sqrt{2+\alpha}\,\alpha^{-3/2}[-\operatorname{sh}H+H\operatorname{ch}H],\\ E_{5x}=\sqrt{2+\alpha}\,\alpha^{-1/2}\left[\alpha^{-1}\operatorname{sh}H+ \left(-\frac12-\frac12\alpha^{-1}\right)\operatorname{sh}2H\right],\\ E_{5y}=\left(-\frac32\alpha^{-1}-\frac32\alpha^{-2}\right) +(1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H+\\ \qquad\qquad +\left(-\frac12\alpha^{-1}-\frac12\alpha^{-2}\right)\operatorname{ch}2H,\\ E_{6z}=\left(-\frac32\alpha^{-1}-\frac32\alpha^{-2}\right) +(1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H+\\ \qquad\qquad +\left(-\frac12\alpha^{-1}-\frac12\alpha^{-2}\right)\operatorname{ch}2H,\\ E_{7x}=\left(\frac12\alpha^{-1}+\frac72\alpha^{-2}+\frac92\alpha^{-3}\right) +(-\alpha^{-1}-4\alpha^{-2}-\\ \qquad\qquad -4\alpha^{-3})\operatorname{ch}H+ \left(\frac12\alpha^{-1}+\frac12\alpha^{-2}-\frac12\alpha^{-3}\right)\operatorname{ch}2H+\\ \qquad\qquad +(\alpha^{-2}+3\alpha^{-3})H\operatorname{sh}H,\\ E_{7y}=\sqrt{2+\alpha}\,\alpha^{-3/2} \left[\alpha^{-1}\operatorname{sh}H+ \left(\frac12+\alpha^{-1}\right)\operatorname{sh}2H+\right.\\ \qquad\qquad \left.+(-1-3\alpha^{-1})H\operatorname{ch}H\right]. \end{gathered} \tag{6.2} \]

In the parabolic case we have

\[ \begin{gathered} E_{1x}=-\beta-\beta^3,\qquad E_{1y}=\sqrt2\left(-\frac14\beta^4\right),\\ E_{2x}=\sqrt2\left(\frac12\beta^2+\frac1{12}\beta^4\right),\qquad E_{2y}=-\beta-\frac16\beta^3,\\ E_{3z}=-\beta-\frac12\beta^3,\\ E_{4x}=1+\beta^2-\frac1{12}\beta^4,\qquad E_{4y}=\sqrt2\,\frac13\beta^3,\\ E_{5x}=\sqrt2\left(-\beta-\frac12\beta^3\right),\qquad E_{5y}=1-\frac14\beta^4,\\ E_{6z}=1-\frac14\beta^4, \end{gathered} \tag{6.3} \]

\[ E_{7x}=\frac{1}{2}\beta^2+\frac{1}{3}\beta^4-\frac{1}{40}\beta^6,\qquad E_{7y}=\sqrt{2}\left(\frac{1}{6}\beta^3+\frac{3}{20}\beta^5\right). \]

Theorem 6.1 is easily proved with the aid of equations (3.5), taking into account the expressions (5.1), (5.2).

In equations (6.1) one may pass from integration with respect to \(\beta\) to integration with respect to \(\tau\) or \(H\) by means of the known relations [2]

\[ d\beta=\frac{d\tau}{\sqrt{\mu}\sqrt{\rho_0}}=\alpha^{-1/2}\,dH, \tag{6.4} \]

which follow from the basic relations

\[ \tau=\sqrt{\mu}\sqrt{\rho_0}\,\beta=\sqrt{\mu}\sqrt{\rho_0}\,\alpha^{-1/2}H. \]

The system of equations (5.1) may be used both for constructing analytical theories and for numerical integration for small values of \(H\) (or \(\beta\)), since in it no perturbations of higher orders are discarded. For large values of \(H\) or \(\beta\), in practice it is more advantageous to use the formulas given in Sec. 7.

Let us note that the constants \(K_i\) have a simple geometrical meaning. \(K_1, K_2, K_3\) are, in the first approximation, corrections at periapsis (or apoapsis) to the initial values of the coordinates; \(K_4, K_5, K_6\) are corrections to the initial values of the components of the regularized velocity; \(K_7\) is the correction to the initial value of the time. Thus equations (6.1) will be equations for the corrections to the initial values at periapsis (or apoapsis) of the coordinates, the components of the regularized velocity, and the time, which is also more advantageous in practice than equations for the elements of Keplerian motion [5].

Theorem 6.2. If in equations (6.1) one confines oneself to calculating only perturbations of the first order, then they can be reduced to the form

\[ K_1=\int\left\{R^{(1)}E_{11}+R^{(2)}E_{12}+R^{(0)}\left[E_{1x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+\right.\right. \]

\[ \left.\left.+E_{1y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\right]\right\}\mu\rho_0\,d\beta, \]

\[ K_2=\int\left\{R^{(1)}E_{21}+R^{(2)}E_{22}+R^{(0)}\left[E_{2x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+\right.\right. \]

\[ \left.\left.+E_{2y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\right]\right\}\mu\rho_0\,d\beta, \]

\[ K_3=\int\left\{R^{(0)}E_{3z}\left(\frac{\partial R_0}{\partial x_0}+Z_0\right)\right\}\mu\rho_0\,d\beta, \]

\[ K_4=\int\left\{R^{(1)}E_{41}+R^{(2)}E_{42}+R^{(0)}\left[E_{4x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+\right.\right. \]

\[ \left.\left.+E_{4y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\right]\right\}\sqrt{\mu}\sqrt{\rho_0}\,d\beta, \tag{6.5} \]

\[ K_5=\int\left\{R^{(1)}E_{51}+R^{(2)}E_{52}+R^{(0)}\left[E_{5x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+\right.\right. \]

\[ + E_{5y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\biggr]\biggr\}\sqrt{\mu}\sqrt{\rho_0}\,d\beta, \]

\[ K_6=\int\left\{R^{(0)}E_{6z}\left(\frac{\partial R_0}{\partial z_0}+Z_0\right)\right\}\sqrt{\mu}\sqrt{\rho_0}\,d\beta, \tag{6.5} \]

\[ \begin{aligned} K_7={}&\int\biggl\{R^{(1)}E_{71}+R^{(2)}E_{72} +R^{(0)}\biggl[E_{7x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+{}\\ &\qquad\qquad\qquad\qquad +E_{7y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\biggr]\biggr\}\sqrt{\mu}\rho_0\sqrt{\rho_0}\,d\beta +\int T_0^*\sqrt{\mu}\sqrt{\rho_0}\,d\beta . \end{aligned} \]

Here

\[ R^{(1)}=-\frac{(r_0R_0^*)'}{\mu+r_0R_0^*}\, \frac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}\,r_0}, \tag{6.6} \]

\[ R^{(2)}=\frac{2R_0^*+\dfrac{r_0}{\mu}R_0^{*2}}{(\mu+r_0R_0^*)^2}-\frac{\rho_0^2}{r_0}, \qquad R^{(0)}=\frac{\rho_0^2}{(\mu+r_0R_0^*)^2} \]

and \(\partial R_0/\partial x_0,\ \partial R_0/\partial y_0,\ \partial R_0/\partial z_0,\ X_0,\ Y_0,\ Z_0,\ T_0^*\) are the result of substituting into the corresponding expressions of the unperturbed coordinates and time as functions of the regularizing variable.

The functions \(E_{11},\ldots,E_{72}\) have, in the general case, the form

\[ E_{11}=2\alpha^{-1}(-1+\operatorname{ch}H),\qquad E_{12}=\alpha^{-1/2}(-\operatorname{sh}H), \]

\[ E_{21}=\sqrt{2+\alpha}\,\alpha^{-1/2}\left[2\alpha^{-1}\operatorname{sh}H+(-\alpha^{-1})H+(-1-\alpha^{-1})H\operatorname{ch}H\right], \]

\[ E_{22}=\sqrt{2+\alpha}\left[(-1-2\alpha^{-1})+(1+2\alpha^{-1})\operatorname{ch}H+ \right. \]

\[ \left. +(-1-\alpha^{-1})H\operatorname{sh}H\right], \]

\[ E_{41}=\alpha^{-1/2}\left[(-2-2\alpha^{-1})\operatorname{sh}H+\alpha^{-1}H+(1+\alpha^{-1})H\operatorname{ch}H\right], \]

\[ \tag{6.7} E_{42}=(1+2\alpha^{-1}+2\alpha^{-2})+(-2\alpha^{-1}-2\alpha^{-2})\operatorname{ch}H +(\alpha^{-1}+\alpha^{-2})H\operatorname{sh}H, \]

\[ E_{51}=\sqrt{2+\alpha}\left[-\alpha^{-1}+(1+\alpha^{-1})\operatorname{ch}H\right], \qquad E_{52}=0, \]

\[ \begin{aligned} E_{71}={}&\alpha^{-1/2}\biggl[\left(3\alpha^{-1}+5\alpha^{-2}\right)\operatorname{sh}H +\left(\frac12+\alpha^{-1}+\frac12\alpha^{-2}\right)\operatorname{sh}2H+{}\\ &\qquad\qquad +\left(\alpha^{-1}-3\alpha^{-2}\right)H +\left(-1-4\alpha^{-1}-3\alpha^{-2}\right)H\operatorname{ch}H\biggr], \end{aligned} \]

\[ \begin{aligned} E_{72}={}&\left(-\frac32\alpha^{-1}-5\alpha^{-2}-\frac92\alpha^{-3}\right) +\left(\alpha^{-1}+4\alpha^{-2}+4\alpha^{-3}\right)\operatorname{ch}H+{}\\ &+\left(\frac12\alpha^{-1}+\alpha^{-2}+\frac12\alpha^{-3}\right)\operatorname{ch}2H +\left(-\alpha^{-1}-4\alpha^{-2}-3\alpha^{-3}\right)H\operatorname{sh}H . \end{aligned} \]

If the unperturbed motion is parabolic, then

\[ \begin{gathered} E_{11}=\beta^2,\qquad E_{12}=-\beta,\\ E_{21}=\sqrt{2}\left(-\beta-\frac{1}{6}\beta^3\right),\qquad E_{22}=\sqrt{2}\left(-\frac{1}{2}\beta^2-\frac{1}{12}\beta^4\right),\\ E_{41}=-\beta+\frac{1}{6}\beta^3,\qquad E_{42}=1+\frac{1}{12}\beta^4, \tag{6.8}\\ E_{51}=\sqrt{2}\left(1+\frac{1}{2}\beta^2\right),\qquad E_{52}=0,\\ E_{71}=-\frac{1}{6}\beta^3+\frac{1}{20}\beta^5,\qquad E_{72}=\frac{1}{2}\beta^2+\frac{1}{6}\beta^4+\frac{1}{40}\beta^6 . \end{gathered} \]

The functions \(E_{1x}, \ldots, E_{7y}\) are determined in the general case by formulas (6.2), and in the parabolic case by (6.3).

Proof. Equations (6.5) are obtained from equations (6.1) if one substitutes there the expressions for \(X_0^*, Y_0^*, Z_0^*, T_0^*\) from formulas (2.5). In doing so it is necessary to compute expressions of the type \(D^{-1}(D_{4i}x_0' + D_{5i}y_0')\), \(D^{-1}(D_{4i}x_0 + D_{5i}y_0)\) \((i=1,2,4,5,7)\), whose numerators are divisible by \(r_0\). Therefore, for example,

\[ \frac{1}{D}(D_{41}x_0' + D_{51}y_0')=\frac{\rho_0^2}{r_0}E_{11}, \]

\[ \frac{1}{D}(D_{41}x_0 + D_{51}y_0) = \frac{\rho_0^2}{r_0}\sqrt{\mu}\sqrt{\rho_0}\,E_{12},\ldots \]

This leads to equations (6.5).

Remark 6.1. The first term in the equation for \(K_5\) is easily integrated, since there \(r_0\) cancels. Therefore

\[ K_5=-\frac{\sqrt{\rho_0}}{\sqrt{\mu}}\ln(\mu+r_0R_0^*)+ \]

\[ +\int\left\{R^{(0)}\left[ E_{5x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{5y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\right\}\sqrt{\mu}\sqrt{\rho_0}\,d\beta . \tag{6.9} \]

From this equation one can draw certain conclusions. For example, if in the problem of the motion of a body in the Earth–Moon system one chooses an unperturbed geocentric orbit as the zero approximation, then, at collision with the Moon, \(K_5\) will assume an infinitely large value. Therefore the method of successive approximations cannot be used here. If here one applies the exact equations (6.1), then at collision with the Moon \(x',y'\to 0\) (just as at collision with the Earth [2]), and therefore there should be no singularity at collision with the Moon. In the named problem, in order to construct a unified analytic theory without introducing spheres of action, it is necessary to choose as the zero approximation a more complicated one than the geocentric unperturbed orbit.

The equations given for the quantities \(K_i\) are substantially simplified if one chooses

\[ I=\frac{r}{\mu},\qquad R^*=0. \tag{6.10} \]

Theorem 6.3. If the quantity \(r/\mu\) is taken as the regularizing multiplier and the fundamental system of solutions has the form (4.1) or (4.7), then the quantities \(K_i\) satisfy the equations

\[ K_1=\frac{\rho_0^3}{\mu}\int\left\{ \frac{r^2}{r_0^2}\left[ E_{1x}\left(\frac{\partial R}{\partial x}+X\right)+ E_{1y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+ \frac{\mu^2}{r_0^2}\left(E_{1x}\varepsilon_x^*+E_{1y}\varepsilon_y^*\right) \right\}\,d\beta, \]

\[ K_2=\frac{\rho_0^3}{\mu}\int\left\{ \frac{r^2}{r_0^2}\left[ E_{2x}\left(\frac{\partial R}{\partial x}+X\right)+ E_{2y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+ \frac{\mu^2}{r_0^2}\left(E_{2x}\varepsilon_x^*+E_{2y}\varepsilon_y^*\right) \right\}\,d\beta, \]

\[ K_3=\frac{\rho_0^3}{\mu}\int\left\{ \frac{r^2}{r_0^2}E_{3z}\left(\frac{\partial R}{\partial z}+Z\right)+ \frac{\mu^2}{r_0^2}E_{3z}\varepsilon_z^* \right\}\,d\beta, \]

\[ \begin{aligned} K_4&=\frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int\left\{ \frac{r^2}{r_0^2}\left[ E_{4x}\left(\frac{\partial R}{\partial x}+X\right)+ E_{4y}\left(\frac{\partial R}{\partial y}+Y\right)\right]\right.\\ &\qquad\left. +\frac{\mu^2}{r_0^2}\left(E_{4x}\varepsilon_x^*+E_{4y}\varepsilon_y^*\right) \right\}\,d\beta, \end{aligned} \tag{6.11} \]

\[ \begin{aligned} K_5&=\frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int\left\{ \frac{r^2}{r_0^2}\left[ E_{5x}\left(\frac{\partial R}{\partial x}+X\right)+ E_{5y}\left(\frac{\partial R}{\partial x}+Y\right)\right]\right.\\ &\qquad\left. +\frac{\mu^2}{r_0^2}\left(E_{5x}\varepsilon_x^*+E_{5y}\varepsilon_y^*\right) \right\}\,d\beta, \end{aligned} \]

\[ K_6=\frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int\left\{ \frac{r^2}{r_0^2}E_{6z}\left(\frac{\partial R}{\partial z}+Z\right)+ \frac{\mu^2}{r_0^2}+E_{6z}\varepsilon_z^* \right\}\,d\beta, \]

\[ \begin{aligned} K_7&=\frac{\rho_0^3}{\mu}\,\frac{\sqrt{\rho_0}}{\sqrt{\mu}} \int\left\{ \frac{r^2}{r_0^2}\left[ E_{7x}\left(\frac{\partial R}{\partial x}+X\right)+ E_{7y}\left(\frac{\partial R}{\partial y}+Y\right)\right]\right.\\ &\qquad\left. +\frac{\mu^2}{r_0^2}\left(E_{7x}\varepsilon_x^*+E_{7y}\varepsilon_y^*\right) +\frac{\mu^2}{\rho_0^3}\varepsilon_t^* \right\}\,d\beta . \end{aligned} \]

Proof. Indeed, if \(R^*=0\), then formulas (2.5) give

\[ T^*=0,\qquad X^*=\frac{r^2}{\mu^2}\left(\frac{\partial R}{\partial x}+X\right) \]

(and similarly \(Y^*, Z^*\)). Substituting these expressions into equations (6.1), we obtain the system (6.11).

Equations (6.11), as well as equations (6.1), are completely exact. Therefore they may be used both for constructing analytic theories of the motion of celestial bodies and for numerical integration for small values of \(H\) or \(\beta\). In the latter case their applica-

is no more advantageous than the use of equations (1.1), since the integration step in the regularizing variable will, generally speaking, be larger than in time.

Remark 6.2. If equations (6.1) or (6.11) are differentiated with respect to $\tau$ [or $\beta$, or $H$, formula (6.4)], then differential equations for the quantities $K_i$ are easily obtained. These equations will be more advantageous than the usual equations for Keplerian elements (Lagrange’s or Newton’s equations [5]), since they are applicable to all types of orbits, except rectilinear-parabolic and rectilinear-hyperbolic ones (see § 4). The latter circumstance is of great importance in astrodynamics.

Theorem 6.4. The first-order perturbations in the quantities $K_i$ are determined by the equations

\[ K_1=\frac{\rho_0^3}{\mu}\int \left[ E_{1x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{1y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \]

\[ K_2=\frac{\rho_0^3}{\mu}\int \left[ E_{2x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{2y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \]

\[ K_3=\frac{\rho_0^3}{\mu}\int \left[ E_{3z}\left(\frac{\partial R_0}{\partial z_0}+Z_0\right) \right]\,d\beta, \]

\[ K_4=\frac{\rho_0^3}{\mu}\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int \left[ E_{4x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{4y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \tag{6.12} \]

\[ K_5=\frac{\rho_0^3}{\mu}\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int \left[ E_{5x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{5y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \]

\[ K_6=\frac{\rho_0^3}{\mu}\frac{1}{\sqrt{\mu}\sqrt{\rho_0}} \int \left[ E_{6z}\left(\frac{\partial R_0}{\partial z_0}+Z_0\right) \right]\,d\beta, \]

\[ K_7=\frac{\rho_0^3}{\mu}\frac{\sqrt{\rho_0}}{\sqrt{\mu}} \int \left[ E_{7x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) + E_{7y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta. \]

Proof. If first-order perturbations are calculated, then in the integrands all variables take their unperturbed values, which gives $r=r_0$, and perturbations of the second and higher orders are discarded, i.e. $\varepsilon_t^*=\varepsilon_x^*=\varepsilon_y^*=\varepsilon_z^*=0$. This leads to equations (6.12).

Equations (6.12) are sufficiently simple, and they will be used for constructing an analytic theory of long-range space flights. From equations (6.11) one can obtain estimates for the influence of perturbations of the second and higher orders.

Remark 6.3. To obtain the functions $E_{1x},\ldots,E_{7y}, E_{11},\ldots,E_{72}$ in elliptic motion, it is necessary to use Remark 4.1.

Remark 6.4. The functions $E_{1x},\ldots,E_{7y}, E_{11},\ldots,E_{72}$ are entire analytic functions of $\beta$, whose expansion coefficients in powers of $\beta$ are polynomials in $\alpha$. The indicated functions are also entire with respect to $\tau$.

7. In the preceding section, without special explanation, it was stated that the equations derived would be practically suitable only for un-

large values of \(\beta\) (or \(H\)). The explanation is as follows. If one takes into account the expressions for \(E_{1x}, \ldots, E_{7y}\), then for large values of \(\beta\) (or \(H\)) and for small \(\alpha\) we shall have

\[ K_1 \simeq \frac{\rho_0^3}{\mu}\int \sqrt{2}\left(-\frac{1}{4}\beta^4\right)\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\,d\beta, \]

\[ K_2 \simeq \frac{\rho_0^3}{\mu}\int \sqrt{2}\,\frac{1}{12}\beta^4\left(\frac{\partial R_0}{\partial x_0}+X_0\right)\,d\beta, \]

\[ K_3 \simeq \frac{\rho_0^3}{\mu}\int\left(-\frac{1}{2}\beta^3\right)\left(\frac{\partial R_0}{\partial z_0}+Z_0\right)\,d\beta, \]

\[ K_4 \simeq \frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\int\left(-\frac{1}{12}\beta^4\right)\left(\frac{\partial R_0}{\partial x_0}+X_0\right)\,d\beta, \tag{7.1} \]

\[ K_5 \simeq \frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\int\left(-\frac{1}{4}\beta^4\right)\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\,d\beta, \]

\[ K_6 \simeq \frac{\rho_0^3}{\mu}\,\frac{1}{\sqrt{\mu}\sqrt{\rho_0}}\int\left(-\frac{1}{4}\beta^4\right)\left(\frac{\partial R_0}{\partial z_0}+Z_0\right)\,d\beta, \]

\[ K_7 \simeq \frac{\rho_0^3}{\mu}\,\frac{\sqrt{\rho_0}}{\sqrt{\mu}}\int\left(-\frac{1}{40}\beta^6\right)\left(\frac{\partial R_0}{\partial x_0}+X_0\right)\,d\beta, \]

provided the partial derivatives of the perturbing function and \(X_0, Y_0, Z_0\) are of the same order. Hence

\[ K_1 \simeq \sqrt{\mu}\sqrt{\rho_0}\sqrt{2}\,K_5,\qquad K_2 \simeq -\sqrt{\mu}\sqrt{\rho_0}\sqrt{2}\,K_4. \tag{7.2} \]

If these quantities, close to one another, are substituted into equations (3.3), where \(q_i^{(j)}\) are defined by (4.1) or (4.7), then in computing \(q_i\) it will be necessary to subtract numbers close in magnitude, which is undesirable in practice because of loss of accuracy. Therefore, instead of the quantities \(K_i\), which have the same dimension, we introduce the dimensionless quantities \(C_i\):

\[ C_1=\frac{1}{\rho_0}K_1-\frac{1}{\rho_0}\sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,K_5, \]

\[ C_2=\frac{1}{\rho_0}K_2+\frac{1}{\rho_0}\sqrt{\mu}\sqrt{\rho_0}\sqrt{2+\alpha}\,K_4, \]

\[ C_3=\frac{1}{\rho_0}K_3,\qquad C_4=\frac{\sqrt{\mu}}{\sqrt{\rho_0}}K_4,\qquad C_5=\frac{\sqrt{\mu}}{\sqrt{\rho_0}}K_5, \tag{7.3} \]

\[ C_6=\frac{\sqrt{\mu}}{\sqrt{\rho_0}}K_6,\qquad C_7=\frac{\sqrt{\mu}}{\rho_0\sqrt{\rho_0}}K_7. \]

The formulas for parabolic motion are obtained from these by setting \(\alpha=0\).

Then the fundamental system of particular solutions ceases to be normal and takes the form (zero solutions are omitted)

\[ \begin{aligned} q_1^{*(1)}={}&\rho_0\bigl[(1+2\alpha^{-1}+2\alpha^{-2})+(-2\alpha^{-1}-2\alpha^{-2})\operatorname{ch}H\\ &\qquad\qquad\qquad\qquad\quad+(\alpha^{-1}+\alpha^{-2})H\operatorname{sh}H\bigr], \end{aligned} \]

\[ q_1^{*(4)}=\rho_0\alpha^{-1/2}\operatorname{sh}H, \]

\[ q_1^{*(5)}=\rho_0\sqrt{2+\alpha}\bigl[(1+2\alpha^{-1})-2\alpha^{-1}\operatorname{ch}H+\alpha^{-1}H\operatorname{sh}H\bigr], \]

\[ q_2^{*(1)}=\rho_0\sqrt{2+\alpha}\,\alpha^{-1/2} \bigl[(1+\alpha^{-1})\operatorname{sh}H+(-1-\alpha^{-1})H\operatorname{ch}H\bigr], \]

\[ q_2^{*(2)}=\rho_0\bigl[-\alpha^{-1}+(1+\alpha^{-1})\operatorname{ch}H\bigr], \qquad q_2^{*(4)}=\rho_0\sqrt{2+\alpha}\,(-\operatorname{ch}H), \]

\[ q_2^{*(5)}=\rho_0\alpha^{-1/2}\bigl[(3+\alpha)\operatorname{sh}H+(-2-\alpha)H\operatorname{ch}H\bigr], \]

\[ q_3^{*(3)}=\rho_0\bigl[(1+\alpha^{-1})-\alpha^{-1}\operatorname{ch}H\bigr], \qquad q_3^{*(6)}=\rho_0\alpha^{-1/2}\operatorname{sh}H, \]

\[ q_4^{*(1)}=(\sqrt{\rho_0}/\sqrt{\mu})\,\alpha^{-1/2} \bigl[(-1-\alpha^{-1})\operatorname{sh}H+(1+\alpha^{-1})H\operatorname{ch}H\bigr], \]

\[ q_4^{*(4)}=(\sqrt{\rho_0}/\sqrt{\mu})\,\operatorname{ch}H, \]

\[ q_4^{*(5)}=(\sqrt{\rho_0}/\sqrt{\mu})\sqrt{2+\alpha}\,\alpha^{-1/2} \bigl[-\operatorname{sh}H+H\operatorname{ch}H\bigr], \tag{7.4} \]

\[ q_5^{*(1)}=(\sqrt{\rho_0}/\sqrt{\mu})\sqrt{2+\alpha}\,(-1-\alpha^{-1})H\operatorname{sh}H, \]

\[ q_5^{*(2)}=(\sqrt{\rho_0}/\sqrt{\mu})\alpha^{-1/2}(1+\alpha)\operatorname{sh}H, \]

\[ q_5^{*(4)}=(\sqrt{\rho_0}/\sqrt{\mu})\sqrt{2+\alpha}\,\alpha^{-1/2}(-\alpha\operatorname{sh}H), \]

\[ q_5^{*(5)}=(\sqrt{\rho_0}/\sqrt{\mu})\bigl[\operatorname{ch}H+(-2-\alpha)H\operatorname{sh}H\bigr], \]

\[ q_6^{*(3)}=(\sqrt{\rho_0}/\sqrt{\mu})\alpha^{-1/2}(-\operatorname{sh}H), \qquad q_6^{*(6)}=(\sqrt{\rho_0}/\sqrt{\mu})\operatorname{ch}H, \]

\[ \begin{aligned} q_7^{*(1)}={}&(\rho_0\sqrt{\rho_0}/\sqrt{\mu})\,\alpha^{-1/2} \bigl[(2+4\alpha^{-1}+3\alpha^{-2})\operatorname{sh}H\\ &\qquad+(-2\alpha^{-1}-2\alpha^{-2})H +(-1-2\alpha^{-1}-\alpha^{-2})H\operatorname{ch}H\bigr], \end{aligned} \]

\[ q_7^{*(2)}=(\rho_0\sqrt{\rho_0}/\sqrt{\mu})\sqrt{2+\alpha}\,\alpha^{-1}(-1+\operatorname{ch}H), \]

\[ q_7^{*(4)}=(\rho_0\sqrt{\rho_0}/\sqrt{\mu})\bigl[(1+\alpha^{-1})+(-1-\alpha^{-1})\operatorname{ch}H\bigr], \]

\[ \begin{aligned} q_7^{*(5)}={}&(\rho_0\sqrt{\rho_0}/\sqrt{\mu})\sqrt{2+\alpha}\,\alpha^{-1/2} \bigl[(2+3\alpha^{-1})\operatorname{sh}H\\ &\qquad+(-2\alpha^{-1})H+(-1-\alpha^{-1})H\operatorname{ch}H\bigr], \end{aligned} \]

\[ q_7^{*(7)}=(\rho_0\sqrt{\rho_0}/\sqrt{\mu}). \]

In the parabolic case the fundamental matrix has the especially simple form:

\[ \|D\| = \left\| \begin{array}{cccc} \rho_0\left(1+\dfrac{1}{12}\beta^4\right) & 0 & 0 \\[1.0em] \rho_0\sqrt{2}\left(-\dfrac{1}{3}\beta^3\right) & \rho_0\left(1+\dfrac{1}{2}\beta^2\right) & 0 \\[1.0em] 0 & 0 & \rho_0\left(1-\dfrac{1}{2}\beta^2\right) \\[1.0em] \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}\,\dfrac{1}{3}\beta^3 & 0 & 0 \\[1.0em] \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}\,\sqrt{2}(-\beta^2) & \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}\,\beta & 0 \\[1.0em] 0 & 0 & \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}}(-\beta) \\[1.0em] \dfrac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}} \left(\beta-\dfrac{1}{3}\beta^3-\dfrac{1}{60}\beta^5\right) & \dfrac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}}\,\sqrt{2}\,\dfrac{1}{2}\beta^2 & 0 \\[1.0em] \rho_0\beta & \rho_0\sqrt{2} & 0 & 0 \\[0.8em] \rho_0\sqrt{2}(-1) & \rho_0\beta & 0 & 0 \\[0.8em] 0 & 0 & \rho_0\beta & 0 \\[0.8em] \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}} & 0 & 0 & 0 \\[1.0em] 0 & \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}} & 0 & 0 \\[1.0em] 0 & 0 & \dfrac{\sqrt{\rho_0}}{\sqrt{\mu}} & 0 \\[1.0em] \dfrac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}}\left(-\dfrac{1}{2}\beta^2\right) & \dfrac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}}\,\sqrt{2}\beta & 0 & \dfrac{\rho_0\sqrt{\rho_0}}{\sqrt{\mu}} \end{array} \right\|. \tag{7.5} \]

Using Remark 4.1, in the elliptic case \([a=e-1,\ H=-iE,\ \rho_0=a(1-e)]\) one can obtain

\[ q_1^{*(1)} = \dfrac{a}{1-e} \left[(1+e^2)-2e\cos E-eE\sin E\right], \]

\[ q_1^{*(4)}=a\sqrt{1-e}\sin E, \]

\[ q_1^{*(5)} = a\sqrt{1+e}\left[(-1-e)+2\cos E+E\sin E\right], \]

\[ q_2^{*(1)} = a\sqrt{\dfrac{1+e}{1-e}}\,e(-\sin E+E\cos E), \]

\[ q_2^{*(2)}=a(1-e\cos E),\qquad q_2^{*(4)}=a\sqrt{1+e}\,(1-e)(-\cos E), \]

\[ q_2^{*(5)}=a\sqrt{1-e}\,[(2+e)\sin E+(-1-e)E\cos E], \]

\[ q_3^{*(3)}=a(\cos E-e),\qquad q_3^{*(6)}=a\sqrt{1-e}\,\sin E, \]

\[ q_4^{*(1)}=\frac{\sqrt a}{\sqrt\mu}\,\frac{e}{1-e}(\sin E-E\cos E),\qquad q_4^{*(4)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{1-e}\cos E, \]

\[ q_4^{*(5)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{1+e}\,(-\sin E+E\cos E), \]

\[ q_5^{*(1)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{\frac{1+e}{1-e}}\,e(-E\sin E),\qquad q_5^{*(2)}=\frac{\sqrt a}{\sqrt\mu}\,e\sin E, \]

\[ q_5^{*(4)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{1+e}\,(1-e)\sin E, \]

\[ q_5^{*(5)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{1-e}\,[\cos E+(1+e)E\sin E], \tag{7.6} \]

\[ q_6^{*(3)}=\frac{\sqrt a}{\sqrt\mu}(-\sin E),\qquad q_6^{*(6)}=\frac{\sqrt a}{\sqrt\mu}\sqrt{1-e}\cos E, \]

\[ q_7^{*(1)}=\frac{a\sqrt a}{\sqrt\mu}\,\frac{1}{1-e}\,[(1+2e^2)\sin E-2eE-e^2E\cos E], \]

\[ q_7^{*(2)}=\frac{a\sqrt a}{\sqrt\mu}\sqrt{1-e^2}\,(1-\cos E), \]

\[ q_7^{*(4)}=\frac{a\sqrt a}{\sqrt\mu}\sqrt{1-e}\,e(-1+\cos E), \]

\[ q_7^{*(5)}=\frac{a\sqrt a}{\sqrt\mu}\sqrt{1+e}\,[(-1-2e)\sin E+2E+eE\cos E], \]

\[ q_7^{*(7)}=\frac{a\sqrt a}{\sqrt\mu}(1-e)\sqrt{1-e}. \]

Remark 7.1. Hence, just as from formulas (7.4), one can obtain expansions of \(q_i^{*(j)}\) in powers of \(\beta\). For this it is necessary to make the substitutions, respectively, \(E=\sqrt{1-e}\,\beta\), \(e=1+\alpha\), or \(H=\alpha^{1/2}\beta\), and expand \(\sin E\), \(\cos E\), \(\operatorname{sh}H\), \(\operatorname{ch}H\) in powers of \(\beta\). This gives the same expansions, which for \(\alpha=0\) (\(e=1\)) pass into formulas (7.5). We note that all functions (7.4), (7.5), (7.6) will be entire with respect to \(\beta\) or \(\tau\), and the coefficients of their expansions in \(\beta\) or \(\tau\) will be polynomials in \(\alpha\) or \(e\).

Theorem 7.1. If, as solutions of the equations in variations, one chooses (7.4) or (7.5), then the equations for the quantities \(C_i\) take the form

\[ C_1=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{1x}^{*}(X^{*}+\varepsilon_x^{*})+E_{1y}^{*}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ C_2=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{2x}^{*}(X^{*}+\varepsilon_x^{*})+E_{2y}^{*}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ C_3=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{3z}^{*}(Z^{*}+\varepsilon_z^{*})\right]\,d\beta, \]

\[ C_4=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{4x}^{*}(X^{*}+\varepsilon_x^{*})+E_{4y}^{*}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \tag{7.7} \]

\[ C_5=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{5x}^{*}(X^{*}+\varepsilon_x^{*})+E_{5y}^{*}(Y^{*}+\varepsilon_y^{*})\right]\,d\beta, \]

\[ C_6=\mu \int \left(\frac{\rho_0}{r_0}\right)^2 \left[E_{6z}^{*}(Z^{*}+\varepsilon_z^{*})\right]\,d\beta, \]

\[ C_7=\mu \int \left\{\left(\frac{\rho_0}{r_0}\right)^2 \left[E_{7x}^{*}(X^{*}+\varepsilon_x^{*})+E_{7y}^{*}(Y^{*}+\varepsilon_y^{*})\right] +\frac{1}{\rho_0}(T^{*}+\varepsilon_t^{*})\right\}\,d\beta. \]

Here, in the general case,

\[ E_{1x}^{*}=\alpha^{-1/2}\left[-\operatorname{sh}H+ \left(1+\frac{1}{2}\alpha\right)\operatorname{sh}2H\right], \]

\[ E_{1y}^{*}=\sqrt{2+\alpha}\left[ \frac{3}{2}\alpha^{-1}+(-1-2\alpha^{-1})\operatorname{ch}H+ \frac{1}{2}\alpha^{-1}\operatorname{ch}2H\right], \]

\[ E_{2x}^{*}=\sqrt{2+\alpha}\left[ -\frac{3}{2}\alpha^{-1}+(1+\alpha^{-1})\operatorname{ch}H+ \frac{1}{2}\alpha^{-1}\operatorname{ch}2H\right], \]

\[ E_{2y}^{*}=\alpha^{-1/2}\left[ (-1-\alpha^{-1})\operatorname{sh}H+ \frac{1}{2}\alpha^{-1}\operatorname{sh}2H\right], \]

\[ E_{3z}^{*}=\alpha^{-1/2}\left[ \alpha^{-1}\operatorname{sh}H+ \left(-\frac{1}{2}-\frac{1}{2}\alpha^{-1}\right)\operatorname{sh}2H\right], \]

\[ E_{4x}^{*}=(-2\alpha^{-1}-2\alpha^{-2})+ (1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H -\alpha^{-2}H\operatorname{sh}H, \]

\[ E_{4y}^{*}=\sqrt{2+\alpha}\,\alpha^{-3/2}(-\operatorname{sh}H+H\operatorname{ch}H), \tag{7.8} \]

\[ E_{5x}^{*}=\sqrt{2+\alpha}\,\alpha^{-1/2}\left[ \alpha^{-1}\operatorname{sh}H+ \left(-\frac{1}{2}-\frac{1}{2}\alpha^{-1}\right)\operatorname{sh}2H\right], \]

\[ E_{5y}^{*}= \left(-\frac{3}{2}\alpha^{-1}-\frac{3}{2}\alpha^{-2}\right) +(1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H+ \]

\[ +\left(-\frac{1}{2}\alpha^{-1}-\frac{1}{2}\alpha^{-2}\right)\operatorname{ch}2H, \]

\[ E_{6z}^{*}=\left(-\frac{3}{2}\alpha^{-1}-\frac{3}{2}\alpha^{-2}\right) +(1+2\alpha^{-1}+2\alpha^{-2})\operatorname{ch}H+ \]
\[ +\left(-\frac{1}{2}\alpha^{-1}-\frac{1}{2}\alpha^{-2}\right)\operatorname{ch}2H, \]

\[ E_{7x}^{*}=\left(\frac{1}{2}\alpha^{-1}+\frac{7}{2}\alpha^{-2} +\frac{9}{2}\alpha^{-3}\right) +(-\alpha^{-1}-4\alpha^{-2}-4\alpha^{-3})\operatorname{ch}H+ \]
\[ +\left(\frac{1}{2}\alpha^{-1}+\frac{1}{2}\alpha^{-2} -\frac{1}{2}\alpha^{-3}\right)\operatorname{ch}2H +(\alpha^{-2}+3\alpha^{-3})H\operatorname{sh}H, \]

\[ E_{7y}^{*}=\sqrt{2+\alpha}\,\alpha^{-1/2} \left[\alpha^{-2}\operatorname{sh}H+ \left(\frac{1}{2}\alpha^{-1}+\alpha^{-2}\right)\operatorname{sh}2H+\right. \]
\[ \left.+(-\alpha^{-1}-3\alpha^{-2})H\operatorname{ch}H\right]. \]

In the parabolic case we have

\[ E_{1x}^{*}=\beta,\qquad E_{1y}^{*}=\sqrt{2}(-1), \]

\[ E_{2x}^{*}=\sqrt{2}\left(1+\frac{3}{2}\beta^{2}\right),\qquad E_{2y}^{*}=-\beta+\frac{1}{2}\beta^{3}, \]

\[ E_{3z}^{*}=-\beta-\frac{1}{2}\beta^{3}, \]

\[ E_{4x}^{*}=1+\beta^{2}-\frac{1}{12}\beta^{4},\qquad E_{4y}^{*}=\frac{\sqrt{2}}{3}\beta^{3}, \tag{7.9} \]

\[ E_{5x}^{*}=\sqrt{2}\left(-\beta-\frac{1}{2}\beta^{3}\right),\qquad E_{5y}^{*}=1-\frac{1}{4}\beta^{4}, \]

\[ E_{6z}^{*}=1-\frac{1}{4}\beta^{4}, \]

\[ E_{7x}^{*}=\frac{1}{2}\beta^{2}+\frac{1}{3}\beta^{4}-\frac{1}{40}\beta^{6},\qquad E_{7y}^{*}=\sqrt{2}\left(\frac{1}{6}\beta^{3}+\frac{3}{20}\beta^{5}\right). \]

Theorem 7.1 is easily proved if formulas (7.3) are used, taking into account equations (6.1).

Remark 7.2. If in formulas (7.8) one makes the substitution \(H=\alpha^{1/2}\beta\) and expands the functions \(E_i^{*}\) with lower indices in series in powers of \(\beta\), then all negative powers of \(\alpha\) cancel. In this case the functions \(E_{1x}^{*},\ldots,E_{7y}^{*}\) are represented by absolutely convergent series in \(\beta\) (or \(\tau\)) for arbitrary values of the latter (i.e., they will be entire functions of \(\beta\) and \(\tau\)), and the coefficients of the expansions in \(\beta\) will be polynomials in \(\alpha\). For \(\alpha=0\) the resulting series pass into the polynomials (7.9).

Theorem 7.2. If \(I=r/\mu\) is chosen as the regularizing factor, i.e. \(R^{*}=0\), then equations (7.7) for the quantities \(C_i\) pass into the equations

\[ \begin{aligned} C_1 &= \frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{1x}\left(\frac{\partial R}{\partial x}+X\right)+E^*_{1y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+\right.\\ &\hspace{45mm}\left.+\frac{\mu^2}{r_0^2}\left(E^*_{1x}\varepsilon^*_x+E^*_{1y}\varepsilon^*_y\right)\right\}\,d\beta, \end{aligned} \]

\[ \begin{aligned} C_2 &= \frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{2x}\left(\frac{\partial R}{\partial x}+X\right)+E^*_{2y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+\right.\\ &\hspace{45mm}\left.+\frac{\mu^2}{r_0^2}\left(E^*_{2x}\varepsilon^*_x+E^*_{2y}\varepsilon^*_y\right)\right\}\,d\beta, \end{aligned} \]

\[ C_3=\frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{3z}\left(\frac{\partial R}{\partial z}+Z\right)\right]+\frac{\mu^2}{r_0^2}E^*_{3z}\varepsilon^*_z\right\}\,d\beta, \]

\[ \begin{aligned} C_4 &= \frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{4x}\left(\frac{\partial R}{\partial x}+X\right)+E^*_{4y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+\right.\\ &\hspace{45mm}\left.+\frac{\mu^2}{r_0^2}\left(E^*_{4x}\varepsilon^*_x+E^*_{4y}\varepsilon^*_y\right)\right\}\,d\beta, \end{aligned} \tag{7.10} \]

\[ \begin{aligned} C_5 &= \frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{5x}\left(\frac{\partial R}{\partial x}+X\right)+E^*_{5y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+\right.\\ &\hspace{45mm}\left.+\frac{\mu^2}{r_0^2}\left(E^*_{5x}\varepsilon^*_x+E^*_{5y}\varepsilon^*_y\right)\right\}\,d\beta, \end{aligned} \]

\[ C_6=\frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{6z}\left(\frac{\partial R}{\partial z}+Z\right)\right]+\frac{\mu^2}{r_0^2}E^*_{3z}\varepsilon^*_z\right\}\,d\beta, \]

\[ \begin{aligned} C_7 &= \frac{\rho_0^2}{\mu}\int\left\{\frac{r^2}{r_0^2}\left[E^*_{7x}\left(\frac{\partial R}{\partial x}+X\right)+E^*_{7y}\left(\frac{\partial R}{\partial y}+Y\right)\right]+\right.\\ &\hspace{35mm}\left.+\frac{\mu^2}{r_0^2}\left(E^*_{7x}\varepsilon^*_x+E^*_{7y}\varepsilon^*_y\right)+\frac{\mu^2}{\rho_0^3}\varepsilon^*_t\right\}\,d\beta, \end{aligned} \]

where the functions \(E^*\) with subscripts are determined by formulas (7.8), (7.9). To obtain the equations in the elliptic case, we use Remark 4.1.

The equations for first-order perturbations will have the form

\[ C_1=\frac{\rho_0^2}{\mu}\int\left[E^*_{1x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+E^*_{1y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\right]\,d\beta, \]

\[ C_2=\frac{\rho_0^2}{\mu}\int\left[E^*_{2x}\left(\frac{\partial R_0}{\partial x_0}+X_0\right)+E^*_{2y}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right)\right]\,d\beta, \tag{7.11} \]

\[ C_3=\frac{\rho_0^2}{\mu}\int\left[E^*_{3z}\left(\frac{\partial R_0}{\partial z_0}+Z_0\right)\right]\,d\beta, \]

\[ C_4=\frac{\rho_0^2}{\mu}\int\left[ E_{4x}^{*}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) +E_{4y}^{*}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \]

\[ C_5=\frac{\rho_0^2}{\mu}\int\left[ E_{5x}^{*}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) +E_{5y}^{*}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta, \]

\[ C_6=\frac{\rho_0^2}{\mu}\int\left[ E_{6z}^{*}\left(\frac{\partial R_0}{\partial z_0}+Z_0\right) \right]\,d\beta, \]

\[ C_7=\frac{\rho_0^2}{\mu}\int\left[ E_{7x}^{*}\left(\frac{\partial R_0}{\partial x_0}+X_0\right) +E_{7y}^{*}\left(\frac{\partial R_0}{\partial y_0}+Y_0\right) \right]\,d\beta. \tag{7.11} \]

If the equations for first-order perturbations in elliptic motion are derived with the aid of Remark 4.1, then instead of equations (7.11) we obtain

\[ C_1=\frac{a^2}{\mu}\int\left\{ (1-e)\left[-\sin E+\left(\frac{1}{2}+\frac{1}{2}e\right)\sin 2E\right]\frac{\partial R_0}{\partial x_0} +\right. \]

\[ \left. +\sqrt{1-e^2}\left[-\frac{3}{2}+(1+e)\cos E-\frac{1}{2}\cos 2E\right]\frac{\partial R_0}{\partial y_0} \right\}\,dE, \]

\[ C_2=\frac{a^2}{\mu}\int\left\{ \sqrt{1-e^2}\left[\frac{3}{2}-e\cos E-\frac{1}{2}\cos 2E\right]\frac{\partial R_0}{\partial x_0} +\right. \]

\[ \left. +\left[e\sin E-\frac{1}{2}\sin 2E\right]\frac{\partial R_0}{\partial y_0} \right\}\,dE, \]

\[ C_3=\frac{a^2}{\mu}\int\left\{ \left[-\sin E+\frac{1}{2}e\sin 2E\right]\frac{\partial R_0}{\partial z_0} \right\}\,dE, \]

\[ C_4=\frac{a^2}{\mu}\int\left\{ \frac{1}{\sqrt{1-e}}\left[-2e+(1+e^2)\cos E+E\sin E\right]\frac{\partial R_0}{\partial x_0} +\right. \]

\[ \left. +\sqrt{1+e}\,[\sin E-E\cos E]\frac{\partial R_0}{\partial y_0} \right\}\,dE, \]

\[ C_5=\frac{a^2}{\mu}\int\left\{ \sqrt{1+e}\left[-\sin E+\frac{1}{2}e\sin 2E\right]\frac{\partial R_0}{\partial x_0} +\right. \tag{7.12} \]

\[ \left. +\frac{1}{\sqrt{1-e}}\left[-\frac{3}{2}e+(1+e^2)\cos E-\frac{1}{2}e\cos 2E\right]\frac{\partial R_0}{\partial y_0} \right\}\,dE, \]

\[ C_6=\frac{a^2}{\mu}\int\left\{ \frac{1}{\sqrt{1-e}}\left[-\frac{3}{2}e+(1+e^2)\cos E -\frac{1}{2}e\cos 2E\right]\frac{\partial R_0}{\partial z_0} \right\}\,dE, \]

\[ C_7=\frac{a^2}{\mu}\int\left\{ \frac{1}{(1-e)\sqrt{1-e}}\left[ \left(-\frac{3}{2}-\frac{5}{2}e-\frac{1}{2}e^2\right)+ \right.\right. \]

\[ \begin{aligned} &+ (1 + 2e + e^2)\cos E + \left(\frac{1}{2} + \frac{1}{2}e - \frac{1}{2}e^2\right)\cos 2E + {}\\ &\quad + (2 + e)E\sin E \bigg]\frac{\partial R_0}{\partial x_0} + \frac{\sqrt{1+e}}{1-e}\bigg[\sin E + {}\\ &\quad + \left(\frac{1}{2} + \frac{1}{2}e\right)\sin 2E + (-2-e)E\cos E\bigg]\frac{\partial R_0}{\partial y_0}\bigg\}\,dE . \end{aligned} \]

Hence, from formulas (7.6), for

\[ e=-1 \tag{7.13} \]

one obtains formulas for the perturbations in the elliptic case of rectilinear motion.

The equations derived in this section will also be used as working equations for constructing an analytic theory of distant space flights (to the Moon, to the major planets of the Solar System, etc.).

The equations for the perturbations in the general case of rectilinear motion will be presented in a separate paper.

In conclusion I express my deep gratitude to Prof. G. N. Duboshin and to the staff of the Department of Celestial Mechanics and Gravimetry of Moscow State University, where I reported this work, for their valuable comments. In addition, this work was part of a report delivered by me in Riga at the conference on the theory of motion of artificial celestial bodies in May 1964. I also express my appreciation to the participants of that conference.

References

1. Siegel C. L. Vorlesungen über Himmelsmechanik. Springer—Verlag, 1956. (A Russian translation is available: Siegel K. L. Lectures on Celestial Mechanics. IL, 1959).

2. Yarov-Yarovoi M. S. Differential Equations, 1, No. 7, 962—976, 1965.

3. Aksenov E. P. Intermediate orbits of artificial Earth satellites. Cosmic Research, 2, issue 1, 6, 1964.

4. Stepanov V. V. A Course of Differential Equations. Fizmatgiz, 1958, p. 282.

5. Duboshin G. N. Celestial Mechanics. Basic Problems and Methods. Fizmatgiz, 1963.

Received by the editors
18 January 1965

P. K. Shternberg State Astronomical Institute