by Joachim Schubart, Astron. Rechen-Institut,
University of Heidelberg, Heidelberg, Germany
schubart@ari.uni-heidelberg.de
Within my earlier Web page Results on Fictitious Orbits at the Inner Jovian 2/1 Resonance, hereafter called Part 1, in Section 4 I have reported on numerical studies of three orbits that approximate the evolution of the orbits of real asteroids. These orbits show libration of the critical argument, crarg, of the 2/1 resonance. Here I describe a way to approximate the long-period evolution of orbital elements of such orbits by means of a Fourier series that depends on trigonometric functions of two independent angles. The terms of the series are multiplied by constant coefficients. The angles change with the time. One of the angles rotates much faster than the other one. It turns out that the use of a modified scale of time is useful for such cases. In the present Part 2 the integration of the orbits depends on simple versions of the three-body problem. Short-period effects do not appear or can be removed by a suitable filter.
Again TL designates the period of libration of crarg and TP refers
to the period of the revolution of Dlp that is equal to the difference of the
longitudes of perihelion of asteroid and Jupiter.
If values of the
semi-major axis, a, of one of the orbits are plotted versus the time, T,
oscillations about a mean value appear and correspond to the cycles of TL.
However, the length of TL varies during a cycle of the much longer period TP.
According to this the oscillations of a are not equally spaced in the plot.
They are closer to each other at Dlp = 0 than at Dlp = 180 deg. Now I introduce
another measurement of time by a variable T2 that proceeds faster than T at
Dlp = 0 and slower at T = 180 deg. If T2 - T is suitably chosen, the
oscillations of a are equally spaced in a plot versus T2. The idea of using a
modified scale of time is not new. In Sections 2 and 3
I apply it to a numerical study of one of the three
orbits that evolves in a nearly quasi-periodic way and depends on two long
periods of quite different length. The new scale is valuable for
an approximation of the elements by a limited number of terms of the
Fourier series. In Section 4 I study two other examples by this method.
a and aj are the semi-major axes of asteroid and Jupiter
e = eccentricity, lp = longitude of perihelion, lm = mean longitude
ej, lpj, and lmj are the respective quantities of Jupiter
Then: Dlp = lp - lpj and: crarg = 2*lmj - lm - lp
* is the multiplication sign, ** indicates exponentiation
In Section 4 of Part 1 the third orbit depends on starting values of the
real asteroid (78801). Both the values of e and the ratio TP/TL are large
in this case. Here I study the long-period effects of this object by a
comparatively simple model. I neglect the attraction of Saturn and turn the
orbital plane of the asteroid into that of Jupiter, the only attracting planet.
The study proceeds on the basis of the planar elliptic restricted three-body
problem, but the equations of motion are further simplified by an averaging
process that removes the effects of the only short-period argument.
I have mentioned this simplified model in Section 6 of Part 1, and that
I have used the model to integrate several fictitious orbits and even the orbit
studied now. For references on earlier work done by this model go to Part 1.
I call the orbit based on the simplified model 78801A. It shows
a nearly quasi-periodic evolution with time. Since there are no effects of
short period, the periods TP and TL rule the evolution. The periods are about
equal to 9408.2 and 421.66 yr, respectively. TP/TL is close to 22.312 .
I have started the orbit with the values a = 3.2655179 au, e = 0.3749632
and with the angular variables (lm-2*lmj) = 272.38 deg , lp = 129.9714 deg ,
lpj = 0 . The orbit of Jupiter depends on aj = 5.201768 au, ej =0.048974 .
The ratio of the masses of sun and Jupiter is equal to 1047.36126 .
In Part 1 I have mentioned the large shifts that correspond to phi in case
of (78801). In the present case they can reach about +/-300 deg. To
approximate phi I use the expression, here with phi in deg and alpha
in radians,
phi = 301.9*sin(alpha-0.091*sin alpha) + 3*sin(alpha-1.1*sin alpha)
+ 2*sin3*alpha - 0.5*sin4*alpha
According to these large shifts an attempt to approximate the orbital elements
by a Fourier series with respect to the angles alpha and b1 will produce terms
of higher order with large coefficients. However a transformation to other
basic angles that includes the change from T to a modified scale of time can
lead to a more favourable situation.
Let T2-T be equal to phi/c4 and b2 = b1 + phi with phi in the above
approximation, so that b2 is nearly equal to beta. Then: c4*T = c4*T2 - phi
and: b2 = c3 + c4*T2 , alpha = c1 + c2*T2 - phi*c2/c4 .
If: mu = alpha + phi*c2/c4 , then mu = c1 + c2*T2 and b2 = c3 + c4*T2
are linear functions of T2. b2 = 0 corresponds to a moment of T2, when
crarg passes zero and increases during an assumed mean libration of this angle.
Since TP/TL is large, c2/c4 is small so that the term phi*c2/c4 only is a
comparatively small addition to alpha. Now beta, a measure of the cycles of
libration,
is close to a linear function of T2. This will allow a study of the orbital
elements by a Fourier series that is based on mu and b2 as functions of T2,
and that shows a fast convergence.
Since the results of the integration of orbit 78801A are stored at equally
spaced intervals of T, a first step consists in a process of interpolation
between the stored values to provide a new tabulation that proceeds with
equally spaced intervals of T2. For each value of T2 the corresponding value
of T is implicitly given and determined at first.
The new tabulation lists values of the orbital elements that are ruled by
the two long periods and that are nearly
quasi-periodic functions of T2. Therefore a Fourier approximation of these
functions is possible in a limited interval. According to a special quality
of the equations of motion that is preserved by the transformation, I expect
that a trigonometric expansion of an orbital element will contain only
a part of the terms that represent a Fourier series in general,
if the constants c1, c2, c3, c4 are well adjusted. In case of a
and of e the series proceeds with terms that depend on
cos(m*b2)*cos(n*mu) and sin(m*b2)*sin(n*mu) with m,n = 0,1,2,... , but in case
of crarg and of lp the terms depend on cos(m*b2)*sin(n*mu) and
sin(m*b2)*cos(n*mu) .
These trigonometric products are multiplied by coefficients to be determined. For this the usual way consists in multiplying each value of the studied orbital element by the value of the trigonometric product in view that results for the corresponding T2. For example, I assume that a is given by a sum of terms of the kind mentioned before. To get the coefficient of a term cos(m*b2)*cos(n*mu) with m and n greater than zero, I multiply all the stored values of a by four times the value of this product that results for the respective T2. The resulting sequence of values represents a function of T2 that is equal to the wanted coefficient plus a sum of periodic terms. Some of these terms have very long periods. I tabulate this sequence.
It is possible to remove the periodic terms from the values of the sequence with a suitable filter by a method mentioned in Part 1. However, terms with a very long period can pass the filter together with the coefficient. Here I prefer a method of averaging that only uses the values of a restricted interval of T2. To separate the coefficient from the sum of periodic terms by averaging, I make use of the near equality of 3*TP and 67*TL. I derive mean values of the tabulated values of the sequence in intervals that nearly cover 67*TL. I shift the beginning of these intervals in an interval that nearly equals 3*TP. Then the mean of all the means of the intervals is nearly free from an influence of most of the periodic terms and I can use it as an approximation of the wanted coefficient. Some long-period terms are less well reduced in amplitude, but they enter the process with small coefficients. It is possible to reduce these coefficients, if the members of the series with the largest determined coefficients are subtracted from the studied element. Then another process of averaging follows. Such a process makes use of stored values from an interval of about 56500 yr.
Here I collect some results from an application of the methods described in the preceding section to the orbit 78801A. I approximate the difference a - 3.27486 , in au, by the sum of terms ca*cos(m*b2)*cos(n*mu) and of terms cb*sin(m*b2)*sin(n*mu) with m,n = 0,1,2,... , with coefficients ca and cb that are functions of m and n . In Table 1 I list the largest coefficients in units of the fifth decimal (in au). The values of ca with m=1 describe an oscillation by TL with changes of amplitude by TP. The cb with m=2 indicate small shifts in the phase of this oscillation. With more accurate values and additional terms the sum does not differ from a by more than three units of the seventh decimal in an interval of about 60000 yr.
Table 1.
10**5 * ca 10**5 * cb
m = 0 1 2 3 m = 1 2 3 4
n = 0 0 1640 0 -12 n = 1 3 136 1 -1
1 -18 185 7 12 2 -3 -40 1
2 -13 -43 -5 -1 3 0 19 -1
3 1 14 3 4 -3 -9
4 -1 -6 -2 5 4
5 2 1 6 -2
6 -1 -1 7 1
8 -1
lp revolves in the backward direction, but lp+mu oscillates. I approximate lp+mu, in radians, by the sum of terms cc*cos(m*b2)*sin(n*mu) and of terms cd*sin(m*b2)*cos(n*mu) with m,n = 0,1,2,... , with coefficients cc and cd that are functions of m and n . In Table 2 the values of cc with m=0 demonstrate an oscillation by TP. The cd with m=1 show a comparatively small oscillation by TL. If a more accurate and extended computation is compared with the studied argument, a small effect by the nearly resonant term sin(b2-22*mu) appears. When the computation includes this effect, the sum does not deviate from lp+mu by more than 6 units of the seventh decimal.
Table 2.
10**5 * cc 10**5 * cd
m = 0 1 2 3 m = 1 2 3
n = 1 -17271 19 10 2 n = 0 -230 -26 1
2 -2318 136 1 -1 1 -62 9 -1
3 480 -7 1 1 2 -36 -11
4 -183 14 3 6 -1
5 64 -5 4 -6 -1
6 -25 3 5 3
7 9 -1 6 -2
8 -4 1 7 1
9 2 8 -1
10 -1
In analogy to the preceding example the coefficients determined for e show an oscillation due to TP and a smaller one due to TL. The coefficients that correspond to the cb of a are small in case of e. The accuracy of the approximation by the series found for e is well again, but in case of crarg the approximation differs from the tabulated values by up to one unit of the fifth decimal, in radians. This is no surprise, since crarg depends on the mean longitude and that accumulates any small effects in a.
In all the studied trigonometric series the terms of higher order show a satisfactory decrease of the coefficients with increasing values of m and n. Especially the decrease with n for a given value of m that is greater than zero demonstrates the success of the change to new angular variables and to a modified scale of time.
In Section 4 of Part 1 the second orbit depends on starting values of the real asteroid (28459). Again I study this orbit by use of a simple model of the forces and of averaged equations of motion in the way described at the beginning of Section 2. This study leads to the values TP = 7843.4 and TL = 375.30 yr of the basic periods and to the ratio TP/TL = 20.899 . I retain the designations and I apply the methods of Section 2. Now the shifts that correspond to phi are smaller. I put phi = 140.6 * sin alpha - 1.6 * sin 2*alpha , in deg. This allows the transformation to T2 and to the new variables mu and b2. Then I determine the coefficients of the Fourier series of the element a in the former way, making use of the near equality of TP and 21*TL. In this process I need the stored values of an interval of 15700 yr for a derivation of the coefficients. Since the stored values cover 60000 yr, I derive the coefficients from two well separated intervals and use the mean value of the two determinations.
A special consideration is necessary with respect to the number of terms ca*cos(m*b2)*cos(n*mu) and cb*sin(m*b2)*sin(n*mu) that is useful for an approximation of a in the present case. Since the amplitude of libration of crarg is comparatively large, now I have to include terms with larger values of m in the series. However, I have to omit terms with n greater than 10. This is due to the similarity of the frequencies of b2 and 21*mu. The slow argument 21*mu-b2 is restricted to a comparatively small interval during the process of determination of the coefficients. Terms of the types cos(n*mu+m*b2) and cos(n*mu-m*b2) can replace the terms with ca and cb in the series. Now the relation (10*mu+j*b2) + (11*mu-b2-j*b2) = (21*mu-b2) with j = m or -m shows that the two arguments on the left side approximately revolve in opposite directions with nearly the same frequency. Therefore the method cannot separately determine the coefficients of the respective terms. The determination of ca and cb of the original terms with n = 10 will include a part of the effects by terms with n = 11. The same is true for n=9 and n=12, and so on. An analogous consideration refers to crarg and to the other orbital elements. However, a step of iteration is especially valuable in this case: All the terms with coefficients determined in a first step are subtracted from the studied element and the method is applied to the residuals. Corrections to the coefficients result in this way, and the final representation of the elements is similar in accuracy to that of orbit 78801A.
At the beginning of Section 2 I have mentioned the simple model of the forces that produces results without short-period effects and allows a large step of integration. I have used it in case of orbit 78801A and of the preceding example. Now I continue with a study that is based on the rigorous equations of motion of the planar elliptic restricted three-body problem sun, Jupiter, (78801). Here I need 5*10**6 steps to cover an interval of about 60000 yr by integration. I use the starting values of 78801A given above together with lm = 226.6539 and lmj = 157.137 deg, lp = Dlp. These values equal the corresponding starting values of orbit 78801 of Part 1. In the present case the results show short-period effects that depend on combinations of the mean frequencies of the revolutions of Jupiter and asteroid and of other frequencies. To apply the procedure used before I remove these effects from the results on a, e, crarg, and (lp + alpha) by a filter that does not change the first four or five decimals of the amplitudes of the long-period effects. I store the results of this process at equally spaced intervals of time. The study depends on these results. The long periods are similar in length to those of 78801A: TP=9388 yr, TL = 421.9 yr, TP/TL = 22.2519 .
For the transformation from T to T2 I use the formula of phi of 78801A with minor changes. Mainly I change the large coefficient to 306 deg. In the determination of the coefficients of the Fourier series I apply the methods of Section 2, again using the approximate equality of 3*TP and 67 * TL. To get an accuracy in the representation of the studied variables by the series as before, I have to use terms with n up to 19 and to apply a step of iteration in analogy to that described for the preceding example. In case of (lp+alpha) I have included small effects by the terms sin(b2-22*mu) and sin(b2-23*mu) in the computation. If the stored results of the filtering process are then compared with the corresponding values given by the series, it turns out that the accuracy of the approximation is even slightly better than that mentioned for 78801A in Section 3.