CC BY
10
УДК 521.1
DOI 10.19110/1994-5655-2021-6-20-35
С.М. ПОЛЕЩИКОВ
ПОТЕНЦИАЛЫ, ДОПУСКАЮЩИЕ ИНТЕГРИРОВАНИЕ ВОЗМУЩЕННОЙ ЗАДАЧИ ДВУХ ТЕЛ В РЕГУЛЯРНЫХ
КООРДИНАТАХ
г. Сыктывкар polsm@list.ru
Аннотация
Изучается задача разделения переменных в некоторых системах координат, полученных с помощью L-преобразований. Даны потенциалы, допускающие разделение регулярных переменных в возмущенной задаче двух тел. Потенциалы содержат произвольные гладкие функции. Рассмотрен пример потенциала, приводящий к построению явного решения задачи в эллиптических функциях. Выделены случаи ограниченного и неограниченного движения. Приведены результаты численных экспериментов.
Ключевые слова:
возмущенная задача двух тел, L-матрицы, интегрируемость, эллиптические функции
Abstract
The problem of separation of variables in some coordinate systems obtained with the use of L-transfor-mations is studied. Potentials are shown that allow separation of regular variables in a perturbed two-body problem. The potential contains two arbitrary smooth functions. An example of a potential is considered allowing explicit solution of the problem in terms of elliptic functions. The cases of bounded and unbounded motion are shown. The results of numerical experiments are given.
Keywords:
perturbed two-body problem, L-matrices, integrabil-ity, elliptic functions
Introduction
S.M. POLESHCHIKOV
POTENTIALS ALLOWING INTEGRATION OF THE PERTURBED TWO-BODY PROBLEM IN REGULAR COORDINATES
Syktyvkar
The integrable cases of motion equations have great practical value. Their significance is determined by the fact that with the help of their solutions one can analyze the motion. In a number of cases integrable problems are used to construct intermediate orbits [1,2]. One of non-trivial examples of integrated systems is the particle motion in a Newtonian field with additional constant acceleration vector. This had been investigated earlier by a number of authors [3-5] and applied to analysis of space flights with constant jet acceleration. In 1970, this problem was studied using regular coordinates obtained from the KS-matrix [6]. In contrast to [6], in [7] integration of the same problem was performed in regular coordinates obtained with the use of ¿-transformations.
In the present work we considera problem of constructing potentials allowing integration of the equations of motion. The idea of our approach consists in the following. First, a new dynamic system is constructed, having more degrees of freedom than the original one. To do this, an ¿-transformation is applied. The theory of L-matrices and their applications is given in [8, 9]. Using new coordinates, a general potential is selected, allowing separation of variables in the Hamilton - Jacobi equation. After this, an inverse transform to original coordinates is performed, using explicit formulas. As a basis for selecting general potential with the required inte-grability property, a well known Stackel theorem is used [10]. This theorem gives necessary and sufficient conditions for separation of variables for orthogonal Hamilton systems, i.e. systems whose Hamiltonian contains only squares of generalized momentums.
Note that separation of variables depends on a choice of a coordinate system. We consider here three
kinds of coordinate systems: regular, bipolar and spherical. The last two systems are introduced in regular coordinates. Canonical equations in regular coordinates are constructed using arbitrary ¿-transformations from the initial canonical motion equations of the perturbed two-body problem. The new equations have also orthogonal form and are invariant with respect to ¿-similarity transforms. In the nonperturbed case these equations do not have singularity at the attracting center. Due to invariance with respect to some perturbing potentials allowing integrability, one can introduce two additional angular parameters.
As a result of this approach the general solution of original system is represented in parametric form, where fictitious time plays the role of parameter, while the physical time depends on this fictitious time and initial data. This sort of integrability is sometimes called 'Sundman integrability' [11].
As an example of integrable case of the perturbed two-body problem the special kind of potential is given. In this example the explicit solution of the problem in terms of elliptic functions is expressed, and the criterion of bounded motion is formulated.
Notation. Everywhere below vectors are regarded as column vectors, and are given in bold letters. The sign T placed over the vector or matrix symbol denotes transposition. A quantity evaluated at the initial moment of physical or fictitious time is denoted by zero superscript: f (0) = f0.
1. The separation of variables
Let us consider the Hamiltonian function of the perturbed two-body problem
H = H(x, y) = i|y|2 - ^ + V,
¡л = j(m + mg), r = |x|,
(1)
where x = (x1;x2,x3)T is the position vector of the point of mass m with respect to the point of mass m0; y = (yi, y2, y3)T is the generalized impulses (yi = Xi, i = 1,2,3); y is the gravitational constant; V = V(x) is the perturbed potential.
For construction of the equations of motion in regular coordinates we will need the ¿-transformation z = ¿(q)q generated by the ¿-matrix of the fourth order that has the following properties:
¿(q^T(q) = ¿T№(q) = |q|2E V q G R4, (2)
(L(q)p)i = (L(p)q)i, i =1,...,p,
(3)
(¿(q)p)i = -(¿(p)q)i, i = p +1,..., 4 (4) V q, p G R4.
Here E is the unitary matrix. The conditions (2) - (4) simultaneously hold only for p = 1 or p = 3. The quantity p is the rank of ¿-transformation. The following theorem can be proved [8,9]. Theorem 1
An arbitrary ¿-matrix generating ¿-transformation of rank three, has the form
L(q)
( qT K1K4 \ qT K2K4 qT K3K4 qT K4
(5)
/
where orthogonal skew-symmetric matrices Ki,K2, K3, K4 are equal to either
or
Ki = aiiU + a2i V + a3i W,
K4 = aiX + а2У + аз2,
Ki = aiiX + a2iY + a3i Z, K4 = ai U + a2V + a3W.
i = 1,2,3,
1, 2, 3,
(6)
(7)
The triplet of vectors ei = (a1i, a2i, a3i) 1,2, 3, forms an orthonormal basis in R3, and e = (a1,a2,a3)T is an arbitrary unitary vector.
Conversely, the arbitrary four skew-symmetric matrices in the form (6) or (7) define the ¿-matrix by the formula (5).
In the formulae (6) and (7) there are the so-called basic skew-symmetric orthogonal matrices
U :
W
У
0 -1 0 0'
10 0 0
0 0 0 -1
0 0 10
0 0 0 -1 0 0 -1 0
0 1 0 0
1 0 0 0
0 -1 0 0 10 0 0 0 0 0 1 0 0-10
V
X =
z
0 0 -1 0
0 0 0 1
10 0 0
0 -1 0 0
0 0 -1 00 0 0 1
10 01
00 00
0 0 0 -1
0 0 10
0 -1 0 0
10 0 0
The matrices Ki are called generators of the ¿-matrix. If K1,K2, K3,K4 are calculated by the formulae (6) then ¿(q) is called the ¿-matrix of the first type, otherwise the ¿-matrix of the second type.
We transfer from variables t, xi, yi to the new variables t, qj, pj by the formulae
dt
• dr,
x У
Mq)q,
2|q|
^(q)p, q, p e R4
(8)
where the matrix A(q) is found from (5) by rejection of the fourth line:
ЛЫ
qT KiK
qT K2K4 qT K3K4
Considerthe equations of motion in newvariables
qi, Pi dj dr
dK dpj
^ = ^, ^ = , j = 0,1, 2, 3,4 (9) dr
dK
dqj
with the Hamiltonian
1
K =8 |p|2 + po|q|2 + |q|2^c(q), yc(q) = V (x(q)).
In this system the first equation with j = 0 corresponds to transformation of time: dq° = \q\2dr. The variable p° is conjugate to q° and has a constant value. If
x(0)
x
y(0) = y0
(11)
are initial conditions for the variables of the system with the Hamiltonian (1), then, as it is proved in [12,13], with the initial values defined by formulae
= A(q° )q0
x
?o(0) = 0,
Po(0) = -H (x0, y0), p0 = 2ЛТ (q0)y0
(12)
the solution of (9) becomes, under the transformation (8), a solution of the system with the Hamiltonian (1) satisfying the initial conditions (11). The function qTK4p preserves a constant value along solutions of (9), and with the initial conditions from (12), this value is zero [13]. Hence, the equality qTK4p = 0 is the first integral of this system. The variable q° coincides with physical time t.
Note that the systems with Hamiltonian (1) and (10) have different orders. The choice of initial values by the formulae (12) means that there is a special construction of the system (9) for each trajectory of the system with Hamiltonian (1).
Let's pick up the form of potential V, admitting division of variables. For this purpose we shall take advantage of the theorem proved by Stackel [10]. Theorem 2
The system with Hamiltonian
H
E'
i=1
i(qi,■■ ■,qn)(K7)P'i + vi(qi)
admits separation of variables in the Hamilton - Jacobi equation if and only if there is a nonspecial matrix $ of order n which elements psi depend only on qi, such as
(1,0,..., 0f,
Фс
(13)
where c = (01,02,..., cn)T.
In this case the integrals of motion will be
t - ßi = £
Pii(Çi)dqi
i=1
Vfi(qi) '
-ßs = Z
i=1
Vsi(qi)dqi
Vfi(qi) '
(14)
pi
= Vfi(qi), i = 1,.
where fi(qi) = 2(ai^u(qi)+.. .+an<fni(qi)-Vi(qi))\ ai, (i = 1,..., n) is constant. As q° a simple root of the function fi(qi) is taken.
Consider again the separation of variables in regular coordinates qi. The Hamiltonian looks like (10). In this case we have
1
01 = 02 = 03 = 04 = 4,
14
|q|2(p° + Vc(q)) = 4J] Vs(qs).
s=1
The solution of system (13) will be, for example, the matrix
/ 4 0 0 0 \
Ф
-1 1 0 0 0 -1 10 0 0 -1 1
The potential is defined up to a constant. As p° is a constant, we obtain
Vc(q)
1
4|q|
r(V1(qi) + V2(q2) +
+ Vs(qs) + V4(q4 ))■ (15)
Let's find expression for the potential Vc(q) in original coordinates x. We notice that variables xi and r are quadratic forms of the variables q1, q2, q3, q4. Using the ¿-similarity transformation it is possible to choose an L-matrix such as a linear combination B1x1 + B2x2 + B3x3 that will be equal to the sum of squares of qi with some coefficients. Note that for any L-matrix we have r = \q\2. As Vc (q) is to be of the form (15), the required potential in x-coordinates will be the function of the form
V (x) = ^(Ar + B1x1 + B2x2 + B3x3). (16) r
Let's specify a choice of L-matrix with the required property. Introduce the notation
в = yjb2 + B2 + B2, bi
в
B '
-i, i = 1,2,3.
Suppose that the L-matrix is ofthe first type. That is, K1, K2, K3 are calculated by the formula (6); for simplicity we assume that K4 = -Y. Then
Ar + B(b1x1 + 62x2 + 63x3) =
= Ar - BqT
(b1ü11 + b2 a12 + b3a,13)U+
+ (b1a21 + b2a22 + b3a23)V+
+ (b1a31 + b2a32 + b3a33)W
Y q.
Choose the parameters aij of L-matrix in such a way that the following equalities hold:
b1an + b2a12 + b3a13 = 1, b1a21 + b2a22 + b3a23 = 0,
b1a31 + b2a32 + b3a33 = 0.
(17)
Geometrically, the vector i1
(b1.b2.b3 )T,
i3 = (a31,a32,a33)
the solution to this system means that
= (an,ai2,ai3)T coincides with b = and the vectors i2 = (a2i,a22,a23)T, T are orthogonal to b. Moreover,
it follows from the structure of the L-matrix that vectors il, i2, and i3 form a frame. It is evident that the system (17) has infinite number of solutions. We write its general solution. For the first vector we have
il = (bi,b2,b3)T.
2
s
For i2 and i3 we assume, in the case h2 + h2 = 0, that We obtain
12
vbf+bf
(b2 cos a + Ь\Ь3 sin a, —hi cos a+
+ h2h3 sin a, — (h + h2)sin
i3 = — 1 (—b2 sin a + b1b3 cos a, b1 sin a+
Vb2 + b22v
+ b2b3 cos a, —(b2 + b2)cos aj .
If b1 + b2 = 0, then b = (0, 0, b3)T, b3 = ±1. Therefore, we can take the following vectors as the general solution of the system (17):
дФ ,— qi = TT" = V Qi cos Q2, dpi
дФ 1— q2 = = VQi sinQ2, dp2
дФ I-
q3 = = V Q3 cos Q4,
дрз
дФ ,— q4 = = V Q3 sin Q4,
ii = (0,0,Ьз)т, 12 = (
cos a, sin a, 0 T
Pi P2 P3
Pz
1
dp4
дФ =
ßQi = 2VÖT
дФ
(19)
(pi cos Q2 + P2 sin Q2),
dQ2
дФ
vQi(—Pi sin Q2 + P2 cos Q2), 1 :(P3 cos Q4 + P4 sin Q4),
дQз
дФ I— 4 = = V Q?(—P3 sin Q4 + P4 cos Q4).
дQ
4
i3 = h3^ — sin a, cos a, 0^ .
The coordinates Q, Q2, Q3, Q4, obtained from (19), will be called bipolar. From the last four equations we The quantity a £ [0, 2n] plays the role of an arbitrary find pi, p2, p3, p4:
parameter of the general solution.
After choosing the parameters aij, the matrix A(q) is determined uniquely. The solution of (17) gives
1
Й2 1
Vc(q) = 7^0 (Ajqj2 — BqTUYq)
,--P2
Pi = 2Pi\f Qi cosQ2--= sinQ2,
Q1
,— P2
P2 = 2Pi\/Qi sinQ2 + cosQ2,
Q1
P4
--2 ((A + B)q2 + (A + B)q2 +
jqj
+ (A — B )q3 + (A — B)q42).
p3
p4
2P3VQ3 cos Q4---= sin Q4,
Q3
P
2P3VQ3 sin Q4 +—-= cos Q4.
Q3
(20)
Hamiltonian in q-coordinates corresponding to this po- In the new variables the Hamiltonian K becomes tential becomes
к = E4(1 p2+4p°q2+4^2))
i=i
where D1 = D2 = A + B, D3 = D4 = A — B. The canonical system of the equations falls into foursubsys-tems
CdT = 1 Pi, ^T = —2(Po+Di)qi,i = 1,2,3,4. (18) dT 4 dT
_ 1 / p 2 P 2 N
K =^4QiP!2 + Qi +4Q3P32 + _
+ Po(Qi + Q3) + (Qi + Q3)V,
where function V is expressed in terms of Qi.
Similar to the above, consider separation of variables in bipolar coordinates. In the notations of theorem 2 we now have
Ci = Qi, C2
1
4Q
C3 = Q3, C4
1
These systems are equivalent to four harmonious oscillators. Integrals of motion are obtained either from (14), As a solution to (13) one can take the matrix or directly from solving (18). Thus, separation of variables for potential (16) is carried out.
For regular q-coordinates, we introduce a new coordinate system. To preserve the canonical form of equations of motion, we use the canonical transformation with generating function
4Q3
/
Ф
Ф = Pi vQ cos Q2 + P2 vQ
sin Q2 +
+ P3VQ cos Q4 + P4VQ sin Q4.
1 Q
4Q
00 00
0 1 0
Qi
V
00
4QI
1
(21)
/
For the potential V admitting separation of variables, we
1
a
find
the form
V = Ö1TÖ3 i^ + Q^н
+ Q3V s(Qs) + V 4Q4 ) 4Q3
In g-coordinates we obtain the form
1
V2(arctan ^)
^ = |qj2 ^ + g2 )V ^ + g2 )+ 4(g2 + g2 )
V4(arctan Ц)
4(g2 + g2)
I2 V(gi
+ (g2 + g2)V з(д2 + g2) +
+
Passing to x-coordinates, we use the concrete L-trans-formation
xi = 2gig4 + 2д2дз, x2 = — 2gig3 + 2g2g4,
x3 = g2 + g2 — g2 — g2,
(22)
which follows from (5), (6) with K1 = V, K2 = W, K3 = U, K4 = -Y. Taking into account that for any L-matrix the equality r = q2 + q2 + q2 + q2 holds, we obtain
q2 + q2 = 2(r - x3), ql + qi = 2(r + x3)-
The general solution of the first equation is
g3
r — X3
cosф, g4
22 ф e [0, 2п].
r — X3 .
sin ф,
Then
x1 sin ф — x2 cos ф x1 cos ф + x2 sin ф
gi =-^ , —, g2 =
V2y/r — x3
V2y/r—x3
In a similar way we may introduce a parameter, using the second equation,
lr + x3 gi = —2—cos фl, g2 =
фl e [0, 2п].
r + x3
2
sin ф1
As is well known [12], with L-transformation for a point in R3 at a distance r from the origin, there corresponds a point of some circle of radius y/r in R4. The variables qi contain an arbitrary parameter p (or p1), giving parametrization of the given circle. In the original coordinates xi this parameter disappears. Note that
q2 q4
— = tan p1, — = tan p. q1 q3
We therefore assume functions V2, V4 to be constant. Then we arrive at a potential of the form
V(x) = - Gi((r + x3)/2) + G2((r — x3)/2) r
, (23)
where G1, G2 are arbitrary smooth functions. The Hamiltonian in bipolar coordinates for this potential takes
Gi(Qi)
K = Qi(pi + Po + I + ——
- P22
-Pf
Л2 ' Qi J ' 4Qi 2
+ Q3( f + P0 + ^1 + 1 P2
Q3
4Q3 2
In view of the solution (21) for fi from the theorem 2 we have
f (ГЛ \ П ( ai a2 «3 Gi(Qi)
/i(Qi) = 4 Q— 4Q — Q— Po —
f2(Q2) = 2«2,
/3(Q3) = 2(- - pn- G2(Q3) 1 /3(Q3) = 2l Q3 4Q2 P0 Q3 J
Q3
/4(Q4) = 2«4.
Then integrals of motion are obtained by formulas (14).
Let's consider one more case of separation of variables. Introduce the spherical coordinates in q-co-ordinates
q1 = vQ1cos Q2 cos Q4, q2 = VQ1 sin Q2 cos Q4, q3 = v01 cos Q3 sin Q4, q4 = v01 sin Q3 sin Q4.
(24)
We supplement the transformation (24) to obtain a canonical transformation of impulses
P1 =2^/~Q~1 cosQ2 cosQ4P1-
sin Q^ cos Q2 sin Q4
"2---P 4,
"v/QT cos Q4 VQ7
P2 = sin Q2 cos Q4Pi +
+ cos Q2 p sin Q2 sin Q4 p VQi cos Q4 vQi '
P3 = 2vQ cos Q3 sin Q4Pi —
sin Q3 cos Q3 cos Q4
P3 +---P4,
(25)
vQi sin Q4 vQi
P4 = sin Q3 sin Q4Pi +
+ cos Q3 p + sin Q3 cos Q4 p
y/QisinQ4 3 y/Ql 4'
Then in new variables the Hamiltonian will be
K = ^4QiPi2 + -
P22 P32 P42 A 2 + ^ Д ^ +Д +
1 ' Qi cos2 Q4 ' Q1 sin2 Q4 ' QW + P0Q1 + Q1F. In the notations of Stackel theorem we have
C1 = Qi, C2 = C3 =
4Qi cos2 Q4'
1 1
C4 =
4Qi sin2 Q4
4Q
In this case the solution of (13) will be the matrix
Ф
/ Q 00
0 1 0 0 0 1 1 0 0
0
\
V 4Q1
cosl Q4
sin2 Q4 1
(26)
The potential V, admitting separation of variables, can be written as
V = Q QVi(QI) + 4Q-2"Q_
Qi V 4Qi cos2 Q4
1 V 2 (Q2)+
+
4Qi sin2 Q4 In view of relations
1 V 3(Q3) + 4Q- V 4(Q4 )) .
Qi cos2 Q4 = qi + q2 Qi sin2 Q4 = q2 + q2
Qi = |q|2 =
2 = r + X3 2 ' r - X3
tan Q4
q2 + q2 _ 1 - Х3/Г
q2 + q2 1 + x3/r'
Q2 = arctan —, Q3 = arctan — q1 q3
following from (22), (24), and the remarks above, we obtain the required form of potential in x-coordinates
1 r 2 A
V (x) =1 [6'i(r)+rrx;+
+ + G )
r - x3 r r
(27)
where Gi, G2 are arbitrary smooth functions and A, B arbitrary constants.
Now assume that a Hamiltonian (1) with the potential (27) is given. Applying L-transformation (22), we write the new Hamiltonian in q-coordinates as
1
K =8 |p|2 + _Pq | q|2 + Gi(|q|2)+
+
A
+
B
+
q2 + q! q2 + q2
, 1 r (qi + q!- q2 - q2 + G2
|q
2
2
Fulfilling canonical transformation (24), (25), we have
+
k = qi( p2 + pq + Gi(Qi)
2
Qi
1 (P 2 \
+ 77i-+4A +
4Qi cos2 Q4 V 2 J
+
-^ (P2 + 4B) +
4Qi sin2 Q4 V 2 )
1 ( P2 \ + p4 + 4G2(cos2Q4 ^ .
Taking into consideration matrix (26), we then obtain
л(ад = 2( Qi +
«4
4Q
■Pq
— 1 in —
Gi(Qi)
Qi 4Qi Qi
/2(Q2) = 2(a2 - 4A), f3(Q3) = 2(a3 - 4B),
/4 (Q4) = ^ -
a2
a3
cos2 Q4 sin2 Q4
- a4 - 4G2(cos2Q4)j.
The integrals of motion follow from (14).
Note that using an arbitrary ¿-transformations allows one to introduce two parameters into the potentials obtained. Tthese two parameters are determined by some constant unit vector b. For example, instead of (27) one can write
1 r 2 A
V (x) = 1 G1(r)+ 2AT +
v ; r . 1 w r + bTx
2B
1 ( bT x
+-iTF" + " G2 -
r — bT x r V r
2. Integration of the system of equations in a special case
In this section we perform straightforward integration of a system with potential of the form (23) having additional parameters. Namely, consider the potential
V = V (x) = - Ц Gi((r + bT x)/2)+ r\
+ G2((r - bTx)/2)),
(28)
where Gi, G2 are some smooth functions, and b = (bi,b2,b3)T an arbitrary unit vector. Note that the vector b provides two parameters in explicit form. Having in mind only theoretical investigation (integrability problem), one can take b to be the ort along the xi-axis. On the other hand, from the more practical point of view, introducing vector b gives us additional degree of freedom necessary for applied problems of celestial mechanics. In such problems, the axes are usually connected with some special directions (equinox or zenith). Therefore the presence of the vector b in potential (28) allows one to turn the coordinate system at one's will.
As Gi, G2, one can take, for example, functions of the form
V + bT x)k, V - bT x)k, k = 1,2,...
rr
We consider a finite linear combination
N
V =--
(Afc(r + bTx)k + Bk(r - bTx)fc). (29)
fc=1
Here Ak, Bk are constants. Such a potential was considered in [16]. This case leads in general to hyperelliptic integrals.
For an interested reader there is a problem: find a real perturbing potential which can be approximated by functions of the form (29). Note that the combination
— B (r + bT x)2 + B (r — bT x)2 = —B bT x 4r 4r
1
r
gives potential corresponding to a constant force. Applications of such potential were considered in [3-5].
The canonical equations of motion have the form
dxi dt
=У:,
d = — xi — x (Gi((r H bTx)/2)+
+ G-((r — bT x)/2)) + + 1 (G/((r + bTx)/2)(r H bi)+
(30)
+ G3(r — bT x)/2)( ^ — bi)) r
where i = 1, 2, 3 and the sign prime indicates the derivative.
This system is the same as the equation of the perturbed two-body problem
x+ x =¿(G!((r H bTx)/2) —
- G'2((r - bTx)/2)) b+ + ¿ (Gi((r + bTx)/2)+
+ G2((r - bTx)/2))x-
- {°i((r + bTx)/2)+
+ G2((r - bTx)/2)) x.
From this one can see that the perturbation is defined by two forces. The first force is collineartothe fixed vector b, and its module varies in dependence on vector x. The second force is the central one.
We are going to show that the system (30) is integrable in regular variables found by ¿-transformations. Transformation (8) contains an arbitrary L-matrix. A special choice of this matrix allows one to separate the variables in the case of an arbitrary constant unitary vector b.
Let us consider the term in (10) containing Vc(q). In the new variables this becomes
|q|Vc(q) = -Gi((|q|2 +
+ qT (biKi + № + b3K3)K4q)/2)-
-G2((|q|2 - qT(biKi + I2K2 + b3K3)KAq)/2).
We assume that the L-matrix has the first type and K4 = -Y. Then
|q|Vc(q) = -Gi((|q|2 - C)/2)-
where
С = qT
— G3((|q|2 H С)/2), (31)
(b/a// H b-a/! + bзalз)U+ H (bia2i H b2a22 + bзa2з)V+
H (bia3i H b2a32 H b3a33)W У q.
Let's select parameters of ¿-matrixes aij from a system (17). Then
С = qT иУq = qT
( —q/ \ —q-
q3 q4
2 2,2,2 q- — q- h q;- H q-
Substituting the found value C in (31), we obtain
|q|2Vc(q) = -C1(q21 + q2) - (q2 + q2).
It follows that the Hamiltonian (10) is represented in the form of the sum
K = Ki + K2,
where
Ki = 8tel + P)+ Po(q2 + q!) - Gi(q2 + q2),
K2 = °(p4 + p; )+ po(q4 + q3) — G 2 (q- + q3 ).
As the value of p0 is constant, the system (9) splits into two independent subsystems
dqi д K1 dт dpi ' dqi _ дK; dт dpi '
dpi dт
dpi dт
dKi
dqi ' OK; dqi '
1 ,2 ,
З 4.
(32)
(33)
We integrate the system (32) again. In the bipolar coordinates Hamiltonian Ki, and accordingly the system, have the form
_ 1 / p 2 x
Ki = g(4QiP2 + Q-) + poQi - Gi(Qi) ,
dQi = Q Q = dт Ql 1' dт 4Q1,
f = - 2P! + § - po + cl(Q/)-
dP;
dт
(34)
0.
Since the Hamiltonian Ki does not explicitly depend on т and Q2, the system (34) has two integrals,
1 p 2 E
2Qip2 +8007 + P0Qi - ^i^f , (35)
P2 = ci.
Here, E7 and c7 are the constants of integration. Taking these integrals into account, the equation for P1 may be written in the following form
dPi =
1Г = 4Q3 8Q1
El + G/(Qi) - Gl(Ql)
Q
Eliminating dT from equations for P1, Q1 and integrating the resulting equation, we find
P/
A
2Q/
V^i(Qi), Si = ±1,
where
Фl(Ql ) = -c/ + E/Qi + C2Q/ + 8Q/Gi(Qi)
and c2 is integration constant defined by
The formulae of inverse transformation
C2 = 4(P0)2 +
(Q1)2
_E1_ G1 (Q1)
Q0 8 Q0 ■
Due to nonnegativity of Q1, from the first equation of the system (34) it follows that
¿1 = sign Q1
Substituting the derived P1 to the first equation of (34) we find
Qi
T + C3 = 2^1
dQ1
Z
-УФТШ
(36)
Using the continuity principle, the sign before the integral (36) cannot change when $1(Q1) is non-zero. Therefore, the function t(Q1) in this case behaves monotoni-cally. Inverting the integral (36), we obtain Q1 as a function of t; we substitute this function in the second equation of the system (34). Then we get
Q2 = -f
C1 I dT
4 J Q1(T)
0
+ C4, C4 = Q0.
Here 03 and 04 are the integration constants. Thus, the values Q1, Q2, P1 are represented as functions of t. If $1(Q1) is a polynomial, the integral (36) is, in general, hyperelliptic.
The integration of the system (33) is done similarly. As a result we find
Q4 = -f
C5 I dT
4 i Q3 (t)
0
+ C8, C8 = Q°0,
P3
where
2Q3
^2(Q3), ¿2 = sign Q'3, P4 = C5,
$2 (Q3) = - 05 + E2Q3 + 06 Q3 + 8Q3G2(Q3).
Here 06 and E2 are the integration constants defined by the equalities
ce = 4(Pj0 )2 +
5 8G2(Q0)
E2
1
(Q0)2 Q0 Q0 P2
-y = 2Q3P2 + 8Q3 + P0Q3 - G2(Q3). (37) The function Q3(t ) is found by a reversion of the integral
Q3
T + C7 = 202
dQs
V^QT
(38)
Thus, the values Q3, Q4, and P3 are also determined as functions of the variable t. The lower limits £ and n in integrals (36) and (38) are chosen according to the location of Q1 and Q3 with respect to the roots of functions $1(Q1) and $2(Q3), respectively.
q2 + q2, q2 + qi
Q1
Q3
p1 = q1P1 + q2P2
tan Q2 = ^, tan Q4 = ||,
2(q2 + q2) '
3 + q4P
2(q3 + q2) '
P3 = q3P3 + q4P4
P2
P4
-q2P1 + q1P2,
-q4P3 + q3P4
allow to define initial values of the variables Qi and Pi
(i = 1, 2, 3,4).
The values of integration constants c1, c2, c3, c4, and E1 are determined by the initial values of Q1, Q2, P1, P2. These five constant values are connected with each other by the integral (35). In the same way, the constant values c5, c6, c7, c8, and E2 are connected by the integral (37) and are defined by initial values of Q3, Q4, P3, and P4. From p0 = —H(x0, y0) we also find relation E1 + E2 = 8^. One has to add the above relations for c2 and c6 to these connections. Besides, as the bilinear relation qTK4p = 0 is the integral of (9), in our case we have
qT(-Y)p = -q2P1 + qvp2 + q4P3 - q3P4 = 0.
Therefore the equality P2 also holds.
P4, or equivalently c1 = c5,
Applying further the first four formulas (19) and (20), we find qi, pi (i =1,2,3,4) as functions of t. Finally, integrating the two remaining equations of (9), we obtain p° = —H(x°, y°) and physical time expressed through t,
T
t = q0 = J |q|2dT + C9 = t1 + t2,
(39)
where
t1 = Q1(T )dT, t2 = Q3(T )dT, C9 = 0.
Thus, the system (9) is completely integrated and we can, at least locally, find a required trajectory. Here it is necessary to note, that if perturbing potentials G1, G2 in (30) are analytic, then, as it is known from a course of the differential equations, the solution of the problem will also be analytic. Let us suppose that the local inversion of integrals (36), (38) appeared to be a globally determined function. In this case we can conclude, by uniqueness of analytic continuation, that this inversion gives not only local, but also global solution of the problem (30). This is the case when functions G1, G2 are polynomials of degree two or three. In this case (36) and (38) are the elliptic integrals, for which inversion we have the well developed technique of elliptic functions at our disposal; thus, we have found the solution of (30) in explicit form.
2
C
1
T
à
2
T
T
3. Inversion of the integral in elliptic case
In this section we consider one case of functions Gi and G2, which reduces to elliptic integrals. Other cases have been studied in [7, 14,16,17]. Let us take the functions G and G2
G i =
G2 =
A_
+ bT x
B_ i
— bT x
+ Ai (r + bTx) + A2(r + bTx)2, (40)
+ Bi (r - bTx) + B2 (r- bTx)2, (41)
where A-i, Ai, A2, B-i, Bi, and B2 are the parameters of the potential. Then for the functions $ i(Qi) and $2(Q3) in (36) and (38) we have the expressions
$ i(Q i) = ci + EiQ i + C2 Q2 + 32A2Q3,
$2(^3) = С + E2Q3 + c6Qi + 32B2Q3.
Here ci
-c2 +4B_ i
-c2 + 4A_ i
- i, C2 = 16A i - 8po, C5 = - i, C6 = 16Bi - 8po. Firstly, let us note that the variables Qi and Q3 are non-negative by definition, and that from integrals (36) and (38) it follows that the ranges of these variables are determined by the inequalities
фl(QI) > 0, Ф2^3) > 0.
(42)
Let us reverse the integral (36). The number of roots of the polynomial $ i and their positions depend on the value of A2. With A2 = 0 the degree of $ i(Qi) equals to two. The integral (36) is found in elementary functions, so this case is not considered here. We distinguish two cases: A2 < 0, A2 > 0. Let's note the roots of $ i(Qi) as £i, £2, £3. The cases under consideration will be sequentially numbered by parameter iA.
I. Assume that A2 < 0. In this case $ i(-^) > 0, $ i(+ro) < 0. The value $ ^0) = ci may be both positive and negative. For actual motion there should be at least one positive root. The qualitatively different cases of the graph of $ i(Qi) are shown in Fig. 1 and 2. In the case of three real roots (Fig. 2), the axis of ordinates goes between £i, £2 if ci < 0, and left with respect to £i or between £2, £3 if ci > 0.
Fig. 1. The graph of Фх(фх). The case A2 < 0. Рис. 1. График Фх(^х). Случай A2 < 0.
The case iA = 1. Suppose that Ф 1 has one real root £ 1, and that Q0 e (0,£ 1) (Fig. 1). Let's write the integral (36) in the form
T + C3
¿1
Qi
dz
2^—mJ - z)(z2 + bz + c)'
£1
where the square trinomial z2 + bz + c has no real roots and is positive for all z, and
b = £' + 32A? c = b£' + 3EAJ (c> 0)- (43)
Apply the substitution
z = ,i - a
l-cos^, a = 4 /,2 1 + cos ф
a = J ,2 + b,i + c
in the integral and put the notations
о , ,i - Qi ь2 1 ,i + b/2
ф1 =2 arctan \ -, k = - 1 +--
w a 2 V a
li = 2 у/-2aA2.
Fig. 2. The graph of Фх(фх). The case A2 < 0. Рис. 2. График Фх(^х). Случай A2 < 0.
Then we derive
T + c3 = -1.
Ф1
dф
уJl - k2 sin2 ф
(44)
Putting here t = 0, we find an integration constant c3:
c3
sign P0
dф
li 0 ф - k2 sin2 ф '
ф0 = 2 arctan i Qi
Check that k2 < 1. As z2 + bz + c has no real
roots, we have b2 - 4c < 0. Therefore,
(,i + 2)2 <,2 + b,i + c = a2 ^
, i + b/2
< 1.
Hence, \ki| < 1. Reversing the integral (44) derived above, we come to the function
Q i = , i + a -
2a
1 + cn(l i(t + c3); ki)'
It is easy to see that for Qi e (0,£i) the denominator cn(w) + 1 = 0. Calculating the derivative of Qi, we get
6i = — signsn (li (t + c3); ki). For the variable Q2 we
1
a
find
Q2 = +
4£ + a)
+
aC1
2h(£2 - a2)
liO+сз)
du
0
Oc3
1 + П1 cn(u; k1 )
du
1 + n1 cn(u; k1)
+Q
n1 = 1 +
2a
£1 - a
Note that
2
П1 k2 = C
2 1 - k1 = 1
4a£1
> 0.
Therefore for calculating the integral of (1 + n1 cn(u; k1))-1 we apply the formula (341.03) [15]
du
1 + n cn(u; k)
1
1 n2
П u
n2 1
; k) - ng1
n2 = 1, (45)
with
1 n2 - 1
g1(u) = ö1
X ln
2 V k2 + k'2n2 \/k2 + k'2n2 sn(u; k) + Vn2 - 1 dn(u; k)
\/k2 + k'2n2 sn(u; k) - Vn2 - 1 dn(u; k) k'2 = 1- k2.
The case iA — 2. Suppose that $1(Q1) has three real roots 0 < £1 < £2 < £3, and Q° e (0,£1). Let's write (36) as
T + C3
ài
Zi
dz
2V-2A2J y/(£1 - z)£ - z)(£3 - z) '
Qi
Making the substitution p = arcsin ^J(£1 — z)/(£2 — z) and reversing the resulting integral, we find
Q1 = £2 -
(£2 - £1)
cn2(I1 (t + C3); k1)'
where
£1, h=V-2A2(£3 - £1 ),
C3
sign P0
л
dф
1 k2 sin2
if Sin ф
£1 - Q0
ф1
arcsin -
£2 - QÏ
Now we calculate ¿1. We differentiate Q1 and use the formula of double argument for elliptic sine. We have
Q1 = 2M6—6)cn-3(u; k1)(—1)sn(u; k1)dn(u; k^ = = — ¿1(^3—6)cn-4(u; k1 )(1 —^2sn4(u; k1))sn(2u; k{),
where the notation u brevity. Therefore,
11(t + c3) is introduced for
¿1 = sign Q'1 = - signsn (2^1 (t + C3); k1).
Here we note the elliptic integral of the third kind Now we find Q2
as
u
n(u, n; k) = J
dv
1 - nsn2(v; k)
Q2 = Cf + C1(£2 - £1)
For t1 we have
2a
t1 = (£1 + a)T - —
11
4£2 4l1£1£2
- n(^1C3,n1; k1)
П(Ь(т + C3),щ; k1)-
+ Q0, n1 = ^. £1
li (o+C3)
du
1 + cn(u; k1)
11C3
du
For the value of physical time, corresponding to the variable Q1, we have
1 + cn(u; k1)
t1 = £2T +
£1- £2 h
li(o+C3)
du
OiC3
du
J cn2(u;k1) J cn2(u;k1)
00
The integral of (1 + cn(u; k\)) 1 is calculated by the mu|a (313 02H15] formula (341.53) [15]
where the integral of cn 2(u; k1) is calculated by the for-
f dv . dn(u; k)sn(u; k)
J =u-E(u)± 1 ±k) ,(46)
°
where E(u) = E(p; k) is incomplete elliptic integral of the second kind (p = am u).
dv
1
cn2(v; k) 1 - k2
(1 - k2)u-
- E(u) +
dn(u; k) sn(u; k) cn(u; k)
k
1
u
l
1
2
n
u
The case iA
3. The polynomial Ф^^х)
has three real roots $i < $2 < $3, and Q°
(тах{0,(2},(з). Let us write (36) as
T + Сз
¿1
Qi
dz
2V-2Â2J y/(z - $i)(z - 6)(6 - z) '
Î3
The reduction of this integral to the standard form (44) is carried out by the substitution p = arcsin \/(£3 — z)/(6 — £2)■ The result of reversion can be presented in the form
Qi = 6 + (6 — 6)sn2(ii(t + C3); h), where the following notations are used
kl = А /73-) ll
Сз
$3 - $i '
sign P0
y/-2A2 ($3 - ),
Ф0
dф
^Jl - k"2 sin2 ф
where the square trinomial z2 + bz + c > 0 for all z. The coefficients b and c are found by the formulae (43). Applying the substitution
z = + a
l - cos ф l + cos ф '
= \j $ 2 + ь$ 1 + С
and reversing the resulting integral, we come to the function
Qi = $1 - a +
2a
1 + cn(li(T + С3); ki)'
where
k2 = 1 ii-$i±b/2
У li = V2aA2.
С3
sign P0
Ф1
d^
0 l - k 2 sin2 ф
$3 - Qi
ф
arcsin -
$3 - $2
For 6i we find 6i = — signsn (2li(t + c3); ki). Substitute Qi in the formulae for Q2 and ti. We find
Q2 =
С
4l i $3
П(1 i (t + С3),ni ; ki)- П(1 i C3,ni; ki)
+ q2 ,
ф° = 2arctan* ' Qi $i
As above, one can show that k2 < 1. The resulting function Q 1 (t) is unbounded, as it has an infinite number of poles on real straight line, which are found by the formula
4m +2^., . t = ---K (ki) — C3, m G Z.
t i = $3T +
$2 - $3 l i
ll (r+C3)
У sn2(u; ki )du-
llC3
sn2(u;ki)du
The integral of squared elliptic sine is calculated by the formula [15]
u
J sn2(v; k)dv = (u — E(u)).
0
II. Assume further that A2 > 0. Now we have $ ^—ro) < 0, $ > 0, and $ ^0) = ci. The
qualitatively different cases of the graph $ 1(Q 1) are shown in Fig. 3 and 4.
The case iA = 4. The polynomial $1(Q1) has one real root ^ and, accordingly, Q1(0) > max{0,^1}. The graph of $1 (Q1) in this case is shown in Fig. 3. Let us write the integral (36) as
T + С3
Si
Qi
dz
2V2Ä2J y/(z - $i)(z2 + bz + с) ! ii
Fig. 3. The graph Фх(фх). Case A2 > 0. Рис. 3. График Фх(^х). Случай A2 > 0.
Fig. 4. The graph Фх(фх). Case A2 > 0. Рис. 4. График Фх(фх). Случай A2 > 0.
G
a
a
l
i
l
i
i
Furtherwe find that ¿1 = signsn (l 1 (т + c3); k 1 ). For variable Q2 we have
Q2 =
ci т
Reversing (47) and using the inverse substitution, we find the required function
Qi = £ i + (£2 — £ i)sn2(l i(t + C3); ki),
4(£1 - a)
where
ac
2l 1 (W-a2)
ll (t+C3)
du
1 + n 1 cn(u; k 1 )
C3
sign P0
Ф°1
dф
I-1C3
du
1 + n 1 cn(u; k 1 )
+ Q
ф? = arcsin -
^Jl — kf sin2 ф
Q - £1
£2 — £1 '
П1 = 1 —
2a
£1 + a
One can show that ö1 = signsn (211(т + c3); k1). For Q2 we find
Note that
Q2
c1
n
2
n1 — 1
1 +
n1 — 1
(£1 — a)2 4a£1
П(11(т + C3),m; k1) —
< 0 < kf.
4l1£1
— П(11С3,n1; k1)
+ Q0, n1 = 1 — .
£1
Therefore for c^lculating the integral of the function (1 + For the first summand of physical time t in (39) we have ni cn(u; ki))-i the formula (45) is to be applied with
91 =
1 — n2
k2
+ k'2n2
arctan
k2 + k'2n2 sn(u; k) 1 — n2 dn(u; k)
Î1 = £1т +
£2 — £1 li
r ll (t+C3)
k'2 = 1 — k2.
If £i = 0 then ni = —1, and for calculating Q2 the formula (46) should be used. For ti we have
J sn2(u; k1)du— I1C3
—J sn2(u;k1)du
2a
t1 = (£1 — а)т + — l1
ll (t + C3)
du
Let us consider the case iA = 6. Suppose Qo e (max{0, £3}, to). The integral (36) has the form
1 + cn(u; k1)
т + C3
¿1
Qi
dz
I1C3
du
2V2Ä2J V(z — £1)(z — £2)(z — £3) '
(3
1 + cn(u; k1)
Suppose that $i has three real roots £i < £2 < £3. The graph of the function $i(Qi) in this case is given in Fig. 4. This case also splits into two subcases:
£i < Q0 < £2 and £3 < Qi.
The case iA = 5. Suppose that Qi e
Using the substitution p = arcsin y/(z — £3)/(z — £2) we transform this integral to the standard form (47). The resulting reversion of the integral in this case is the following
Qi = £2 + £3 — £2
cn2(l1 (т + C3); k1)'
(max{0,£1},£2). We write (36) as
where
т + C3 =
¿1
Qi
dz
k1 =
2V2Ä2J V(z — £1)(z — £2)(z — £3) ' (l
, l1 = V2A2(£3 — £1),
We apply the substitution ф = arcsin ^J(z — £1)/(£2 — £1) to this integral and use the notations
c3 =
sign P0
Ф1
dф
1 — k12 sin2 ф
k1
h = ^2A2(£3 — £1 ).
ф1 = arcsin
Q — £3 Q? — £2 '
Then our integral has the standard form
+ ¿1 т +С3 = T
Ф1
d^
1 0 \A — k2 sin2 ф
Now the function Qi (t) has an infinite number of poles of the second order, hence it is unbounded. The poles are found by the formula
2m + 1 . .
t = —:-K(ki) — C3, m e Z.
l1
l
1
1
1
l
1
Further we find the values ¿1, Q2, ti
¿1 = signsn (21i(t + C3); ki),
Q2 = £1! + Ci(£2 - £3)
Q2 46 + 4li66
4. Numerical examples and analysis of motions
In the examples below we consider the motion of a particle in perturbed gravitational field of a planet with
spherical density distribution, which gravitational param-
32
П(/1(т + сз) n\-k1)- eter is taken to be / = 398601.3 km3/s 2. The perturb-
0 £2
- n(liC3,rai; ki) + Q2, ni = —,
£3
ti = Ьт +
£3 - £2 li
h h+cs)
du
cn2(u;ki)
hcs
du
cn2(u;ki)
An inversion of the integral (38) is fulfilled in a similar way. This integral and the function $2 differ only by notations from the integral (36) and the function $1. Therefore, after some evident renaming, we find the expressions for Q3, Q4, and t2. We number these cases sequentially by the parameter is. Then we have:
is = 1 o B2 < 0, n2,n3 G C, 0 < Q3 <rn
(Q3 is bounded).
is = 2 o B2 < 0, 0 < Q3 < n1 <V2 <m
(Q3 is bounded).
is = 3 o B2 < 0, <n2 < Q3 <m
(Q3 is bounded).
is =4 o B2 > 0, n2,V3 G C,m < Q3
(Q3 is unbounded).
is = 5 o B2 > 0, n1 < Q3 <V2 <m
(Q3 is bounded).
iB = 6 o B2 > 0, m <^2 <m <Q3
(Q3 is unbounded).
The study above yield the following theorem. Theorem 3
The motion of the particle is bounded if and only if at the initial moment both variables Q1 and Q3 are restricted on the right by the roots of the polynomials $1 and $2, correspondingly.
Now we give a definition of retaining potential, introduced in [14]. Definition 1
A potential is named as retaining, if for arbitrary initial conditions the motion of a particle in a perturbed field corresponding to this potential is bounded.
Thus, potential (28), where G1, G2 are defined by the formulae (40), (41) for A2 < 0 and B2 < 0, is retaining. Generally, the formulae (36), (38) are not elliptic integrals, and we cannot present a solution in explicit form. Nevertheless, the above-stated qualitative result remains true [16].
ing force is defined by the potential (28), with G1 and G2 calculated by the formulae (40) and (41). For convenience (to have no fractions), a dimensionless direction vector b for the constant force is used. While doing calculations, this direction vector is assumed to be normalized. The dimensions of parameters A_1 and B-1 are [km4/s2], A1 and B1 are [km2/s2], A2 and B2 are [km/s2]. Calculations and construction of orbits were performed using the Maple system with 32 digits. In each example, for convenience of its analysis, the values of
circular and parabolic velocities vc
upar
of Keplerian
motion are given. The perturbations under consideration are great, they are non-typical for the Earth's satellites. For this reason, we do not give Keplerian elements of osculating orbits for the corresponding initial values. The initial position of a particle is marked by a point on the corresponding figure.
Example 1. Initial values of coordinates and velocities of a particle:
x1 = 8200 km, x2 = 0 km, x3 = 6000 km,
X2 = 8.6 km/s, X1 = X3 = 0 km/s,
(vcir & 6.26 km/s, vpar & 8.86 km/s).
In an unperturbed case these values define an elliptic motion.
Parameters of potential are as follows:
A_1 = 0.004 km4/s2, A1 = 0.06 km2/s2,
A2 = 0.2 • 10_7 km/s2, B_1 = 0.0001 km4/s2,
B1 = 0.008 km2/s2, B2 = —0.3 • 10_4 km/s2.
Coordinates of direction vector are b = (—1, 2,1)T. In the case under consideration the roots of polynomials $1 and $2 are
£1 & 1478, 6 & 115346, Q1 & 4631 ^ Q1 G (£1,6), n2 & 1707, n3 & 31031,
Q3 & 5529 ^ Q3 G (V2,V3).
Therefore, the motion is bounded. This is the case iA = 5, is = 3.
The coordinates and velocities have been calculated during a time range, corresponding to two revolutions of the particle around the attracting centre without perturbations, that is t g [0, 2T], where T is calculated by the formula
T = -hkk , hk
|x 0|2 m
J—1---—. (48)
2 |x0| ' ( )
Here hk is the Keplerian energy. Let's remind that L-transformation doubles the angles at the origin of coordinates.
Fig. 5. The case iA = 5, iB = 3. Рис. 5. Случай iA = 5, iB = 3.
Note that in this example the potential is not retaining. Nevertheless, the motion appears to be bounded. The trajectory of the particle is shown in Fig. 5.
Example 2. Initial values of coordinates and velocities of a particle:
xi = 8200 km, = 0 km, x3 = 6000 km,
x2 = 9.9 km/s, X1 = X3 = 0 km/s,
(vcir & 6.26 km/s, vpar & 8.86 km/s).
In unperturbed case the motion belongs to hyperbolic type.
Parameters of potential are as follows:
A-1 = 0.004 km4/s2, A1 = 0.006 km2/s2, A2 = -0.2 • 10-7 km/s2, B-1 = 0.0001 km4/s2, Bi = 0.008 km2/s2, B2 = -0.3 • 10-7 km/s2.
As A2 and B2 are negative we have a retaining potential. Coordinates of direction vector are b = (1, 2, — 1)T. The roots of polynomials:
£2 & 2126, £3 & 122192633,
Qi & 5529 ^ Qi e (&,&),
П2 & 1699, П3 & 81506371,
Q3 & 4631 ^ Q3 e (п2,ПЗ).
Therefore, the motion is bounded. This is the case iA = 3, iB = 3. The integration is carried out during the time range corresponding approximately to t = 1759.74 days. The particle trajectory is shown in Fig. 6.
Fig. 6. The case îa = 3, iB =3. Рис. 6. Случай iA = 3, iB = 3.
Example 3. Initial values of coordinates and velocities of a particle are as follows:
x\ = 6000 km, x2 = 0 km, x3 = -8000 km,
X2 = 7.9 km/s, X1 = X3 = 0 km/s,
(vcir ^ 6.31 km/s, vpar ~ 8.93 km/s).
In an unperturbed case these values define an elliptic motion.
Parameters of a potential are as follows:
A-1 = 0.04 km4/s2, A1 = 0.03 km2/s2, A2 = -0.2 • 10-5 km/s2, B-1 = 0.1 • 10-4 km4/s2, Bi = -0.0003 km2/s2, B2 = 0.3 • 10-4 km/s2.
Here the potential is not retaining. Coordinates of direction vector are b = (1,1,1)T. The roots of polynomials are as follows:
£2 « 2686, £3 « 20699,
Qi « 4423 ^ Qi G (£2^3),
ni « 3256, n2,V3 G C,
Q0 « 5577 ^ Q0 >ni-
Therefore, the motion is unbounded. This is the case îa = 3, iB = 4.
The integration is carried during the time range corresponding approximately to t = 3.23 days. The par-
Fig. 7. The case iA = 3, iB =4. Рис. 7. Случай iA = 3, iB = 4.
In the case of unbounded motion, to define the integration interval firstly one has to find the nearest pole of Q1(r) and/or Q3(t) in the direction of ascending т. Suppose this nearest pole is at т = т1. Then we choose a small positive value e and divide the segment [0, т1 —e] into N equal subsegments. The value N is to be selected from practical reasons. The orbit should be visually a smooth curve. In our examples the value N = 100 was used. After that, the calculations by the formulae derived above are carried out in equidistant nodes.
The following example demonstrates an application of our formulae for testing a numerical integration method. The original system of motion equations (30) is considered. The Runge-Kutta-Fehlberg method of the eighth order with automatic choice of integration step is tested. The step is chosen by a method of the seventh order. The corresponding pair of programs, implemented in FORTRAN, is below noted as RKF8(7). Integration of equations (30) was performed by RKF8(7) with relative local error of the method e = 10-13, and all calculations were carried out with double precision (real*8). The gravity parameter and the units of measurement are the same as above. A hypothetical particle is considered, repeatedly encountering the attracting centre. The trajectory obtained by explicit formulae is taken to be standard (reference). Its coordinates have been obtained using Maple with 32 digits (in FORTRAN this corresponds to quadruple precision (real*16)).
Example 4. Initial values of coordinates and velocities of a particle are as follows:
x1 = 7000 km, x2 = 0 km, x3 = 6000 km,
Table 1
Estimation of the precision of numerical integration
Таблица 1
Оценка точности численного интегрирования
n t(day) SH х 10"12 Sx1 х 10"12 Sx2 х 10"12 Sx3 х 10"12 Sr х 10"12
1 .3382444 1 0.2 1 1 0.4
10 4.9080991 2 6 12 213 10
50 24.1940313 41 729 104 1667 108
100 48.4322508 53 4399798 523748 95154 330606
500 242.7821163 294 77898 31418 151259 77206
1000 485.2955201 556 554500 332688 1067003 330900
x2 = 7.9 km/s, x1 = x3 = 0 km/s,
(Vc
6.58 km/s, vpar & 9.30 km/s).
In an unperturbed case we have an elliptic motion.
Parameters of retaining potential are as follows:
A-i = 0.1 km4/s2, Ai = -0.02 km2/s2,
A2 = -0.2 • 10-5 km/s2, B-1 = -0.004 km4/s2,
Bi = -0.001 km2/s2, B2 = -0.001 km/s2.
Coordinates of direction vector are b = (-1, -3,1)T. The roots of polynomials and $2 are as follows:
6 - 764, - 58639, Qi - 4459 ^ Ql e (6,6),
n2 - 504, n3 - 7209, Q3 - 4761 ^ Q3 e (m,V3).
Therefore, the motion is bounded. The case iA = 3,
is = 3.
The calculations were carried out during the time ranges corresponding to 1, 10, 50, 100, 500, and 1000 revolutions of the particle around attracting centre without perturbations. The trajectory of the particle for three revolutions is shown in Fig. 8.
Fig. 8. The case iA = 3, iB = 3. The motion is bounded. Рис. 8. Случай iA =3, iB = 3. Движение ограничено.
Table 1 contains the values, near the end of the integration interval, of the relative error for the energy constant SH, the coordinates of particle position vector xi, and its absolute value r
SH = Hlzli, Sxi = Xrfil, , = 1,2,3,
H
Xi
| r — rc \
where H0 is the value at the initial moment, H at an arbitrary moment; xc, rc are the values found by exact formulae. In the second column the intervals of physical time t (in days) are given, for which numerical integration of system (30) was carried out.
From these data we can see that if the integration interval increases, the relative errors of H and x3 do not decrease. For coordinates x1, x2, and absolute value r, with n = 100, these errors increase, then they diminish, and then increase again.
The numerical examples show efficiency of the formulae we obtained. Besides, the theorem 3 allows to determine, given the initial position and velocity of a particle, whether its motion is bounded or unbounded.
Conclusion
In this paper we consider three sorts of coordinates (regular g-coordinates, bipolar coordinates, spherical coordinates). For each of the systems, the forms of potentials admitting complete separation of variables are given. Thus, the original equations for such potentials allow integration "in the sense of Sundman". In a similar way one can build, for regular g-coordinates, other coordinate systems for which Hamiltonian has orthogonal form, and with the use of Stackel theorem build potentials allowing the above-mentioned integrability.
Application of these potentials is a separate and independent problem. These potentials are of practical importance, which approximate some real forces.
The author is grateful to professor A.Zhubr for useful comments and discussions.
Литература-References
1. Aksenov E.P. Teorija dvizhenija iskusstvennyh sputnikov zemli [Theory of the motion of the Earth's artificial satellites]. Moscow: Nauka, 1977. 360 p.
Аксенов Е.П. Теория движения искусственных спутников земли. М.: Наука, 1977. 360 с.
2. Ferrandiz J.M., Floria L. Towards a systematic definition of intermediaries in the theory of artificial satellites // Bull. Astron. Inst. Czechosl. 1991. Vol. 42. P. 401 - 407.
3. Beletsky V.V. Traektorii kosmicheskih poletov s postojannym vektorom reaktivnogo uskorenija [Trajectories of Space Flights with a Constant Vector of Reactive Acceleration] // Cosmic Research. 1964. Vol. 2, No. 3. P. 787 - 807. Белецкий В.В. О траекториях космических
полетов с постоянным вектором реактивного ускорения // Космические исследования. 1964. Т. 2, № 3. С. 787 - 807.
4. Kunitsyn A.L. O dvizhenii rakety v central'nom silovom pole s postojannym vektorom reak-tivnogo uskorenija [Rocket Motion in a Central Force Field with a Constant Vector of Reactive Acceleration] // Cosmic Research. 1966. Vol. 4, No. 2. P. 324 - 332.
Куницын А.Л. О движении ракеты в центральном силовом поле с постоянным вектором реактивного ускорения // Космические исследования. 1966. Т, 4. № 2. С. 324 - 332.
5. Demin V.G. Dvizhenie iskusstvennogo sputnika v necentral'nom pole tjagotenija [The motion of an artificial satellite in the eccentric gravitational field]. Moscow: Nauka, 1968. 352 p.
Демин В.Г. Движение искусственного спутника в нецентральном поле тяготения. М.: Наука, 1968. 352 с.
6. Kirchgraber U. A problem of orbital dynamics, which is separable in ^S-variables // Celest. Mech. 1971. Vol. 4. P. 340 - 347.
7. Poleshchikov S.M. One integrable case of the perturbed two-body problem // Cosmic Res. Vol. 42, No. 4. P. 398 - 407.
8. Poleshchikov S.M., Kholopov AA Teorija L-matric i reguljarizacija uravnenij dvizhenija v nebesnoj mehanike [Theory of L-matrices and regularization of motion equations in Celestial Mechanics]. Syktyvkar: Syktyvkar Forest Inst., 1999. 255 p.
Полещиков С.М., Холопов AA Теория L-мат-риц и регуляризация уравнений движения в небесной механике. Сыктывкар: СЛИ, 1999. 255 с.
9. Poleshchikov S.M. Regularization of motion equations with L-transformation and numerical integration of the regular equations // Celest. Mech. and Dyn. Astr. 2003. Vol. 85, No. 4. P. 341 - 393.
10. Pars LA. A treatise on analytical dynamics. NY: Wiley, 1965. 636 p.
11. Kholshevnikov K.V. Ob integriruemosti v nebesnoj mehanike [On the integrability in celestial mechanics] // Analiticheskaja nebesnaja mehanika [Analytical celestial mechanics]. Kazan: Kazan University publ., 1990. P. 5 - 10.
Холшевников К.В. Об интегрируемости в небесной механике // Аналитическая небесная механика. Казань: Изд-во Казан. ун-та, 1990. С. 5 - 10.
12. Stiefel E., Scheifele G. Linear and regular celestial mechanics. Berlin: Springer-Verlag, 1971. 304 p.
13. Poleshchikov S.M. Regularization of canonical equations of the two-body problem using a generalized ^S-matrix // Cosmic Res. 1999. Vol. 37, No. 3. P. 302 - 308.
14. Poleshchikov S.M. Motion of a particle in a perturbed field of the attracting centre // Cosmic Res. 2007. Vol. 45, No. 6. P. 522 - 535.
15. Byrd P.F., Friedman M.D. Handbook of elliptic integrals for engineers and physicists. Berlin: Springer-Verlag, 1954. 355 p.
16. Poleshchikov S.M., Zhubr A.V. A set of potentials allowing integration of the perturbed two-body problem in regular coordinates // Cosmic Res. 2008. Vol. 46, No. 3. P. 202 - 214.
17. Poleshchikov S.M. Integriruemyj sluchaj voz-mushhennoj zadachi dvuh tel, porozhdajushhej elementarnye funkcii [An integrable case of the perturbed two-body problem producing elementary functions] // Proc. of Syktyvkar Forest Inst. 2006. Vol. 6. P. 31 - 57. Полещиков С.М. Интегрируемый случай возмущенной задачи двух тел, порождающей элементарные функции // Труды СЛИ. 2006. Т. 6. С. 31 - 57.
Статья поступила в редакцию 01.11.2021.
