- Email: [email protected]
Construction of polynomial approximation for Hansen coefficients
Computer Physics Communications 124 (2000) 40–48 www.elsevier.nl/locate/cpc
Construction of polynomial approximation for Hansen coefficients Akmal A.Vakhidov
1
Lohrmann Observatory, Technical University Dresden, Mommsenstrasse 13, Dresden 01062, Germany Received 18 March 1999; received in revised form 17 May 1999
Abstract The problem of constructing efficient approximations of Hansen coefficients using polynomials in terms of the eccentricity is considered. The properties of the following approximating schemes are studied in detail: approximation by fragments of Taylor series, interpolation by Lagrange polynomials, Chebyshev approximation. The accuracy of all these approximating schemes is investigated for different values of eccentricities and for different indices of the Hansen coefficients. 2000 Elsevier Science B.V. All rights reserved. PACS: 95.10 Ce – Celestial Mechanics; 02.30 Mv – Approximations and Expansions Keywords: Hansen coefficients; Polynomial approximation; Computer algebra
1. Introduction If we construct analytical or semianalytical theories of motion of celestial bodies, then in many cases it is very convenient to use the expansions of coordinates into trigonometric series in multiples of the mean anomaly depending on the time linearly. Such an approach gives us the possibility of solving analytically the equations of motion using the method of consequent approximations (see, for example, [1,2]) or to realize the averaging procedure based on von Zeipel or Lie-series methods [3]. Quite generally, the expansion of coordinates in multiples of the mean anomaly in the case of elliptic motion is given by the formula: n ∞ X r exp(j mv) = Xsn,m (e) exp(j sM), a s=−∞ 1 E-mail: [email protected].
where (r, v) are the polar coordinates of the celestial body, a is the semi-major axis of the elliptic orbit, M is√the mean anomaly, n, m, s are some integers, j = −1. The coefficients Xsn,m (called as Hansen coefficients) depend on the eccentricity e of the orbit. The problem of efficient computation of Hansen coefficients was investigated in detail in many papers (see also [4–7]), but this problem is really very complicated and to date we have no standard algorithms for computing these coefficients efficiently for all possible values of indices n, m, s and eccentricity e. A great problem in using the expansions with Hansen coefficients in theories of motion of celestial bodies is connected with the necessity to recompute these coefficients many times for different (very often a little different) values of eccentricity. Indeed, when using a semianalytical theory of motion or theory presented via orbit elements, it is necessary to recompute all Hansen coefficients at each step of the numerical integration of the equations of motion. It
0010-4655/00/$ – see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 9 9 ) 0 0 4 3 1 - 2
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
leads to many inconveniences from the point of view of economy of computing time and decreases the efficiency of the theory itself. Therefore sometimes it seems to be much more efficient to construct some approximations, which give the possibility of computing the Hansen coefficients on the determined interval of values of eccentricity. Such approximations could be obtained before the usage of motion theory and allow one to evaluate the coefficients using compact approximating formulae instead of more complicated direct methods. The accuracy reached in approximations should satisfy to the accuracy of the corresponding analytical or semianalytical theory. For many practical problems of celestial mechanics in theories of the first order the accuracy of computing the Hansen coefficients should be equal to 10−4 . The errors then will be smaller than the terms of the second order, which are not taken into consideration. Such an accuracy is enough for theories of the first order for such problems as the motion of minor planets or the motion of artificial Earth satellites, where the small parameter of the problem is of order 10−3 . If we want to take into consideration the perturbations of the second order then the accuracy of computing the disturbing function (and Hansen coefficients, respectively) should be 10−7 . But of course, the desired accuracy of computing the Hansen coefficients depend on the concrete problem and is different for different problems of celestial mechanics. Therefore it is very important to construct the approximations in a way that could incorporate the simplicity of computing the Hansen coefficients and the desired accuracy simultaneously. In the present paper we investigate three approximating schemes for computing the Hansen coefficients: (1) approximation by fragments of Taylor series; (2) interpolation with the help of Lagrange polynomials; (3) Chebyshev approximation. We present the results of computer experiments to estimate the accuracy of these approximating schemes and advantages and inconveniences of usage of each of them. For numerical experiments we consider the Hansen coefficients both with n < 0 (the case of perturbations from internal body) and with n > 0 (perturbations from external body), different values of indices m and s, and different intervals of eccentricity for approximations. All computer experiments described below were realized on Sun Sparc workstations in the computer alge-
41
bra system MAPLE V (Release 4) [8] with different numbers of digits. The software for constructing the described approximations is available from the author upon request.
2. Approximation of Hansen coefficients by fragments of Taylor series For estimating the accuracy of polynomial approximations of Hansen coefficients we will consider the behaviour of the logarithm of the relative error defined as Xe − Xa , D = log Xe where Xe is the “exact value” of Hansen coefficient and Xa is the value obtained with the help of the approximating polynomial. By the words “exact value” we understand here not the true value computed with absolute accuracy, but the value obtained using direct methods for computing the Hansen coefficients. These “exact values” Xe can be computed, in particular, on the basis of the formulae of Hill’s method [9]. This method is based on the representation of Hansen coefficients via hypergeometric functions and Bessel functions which can be evaluated with arbitrary precision as the standard built-in functions of the computer algebra system MAPLE V [8]. Hill’s method allows also the computation of the Hansen coefficients in symbolic form for the forthcoming expansion of the corresponding analytical expressions of Hansen coefficients into Taylor series. This procedure can also be realized by means of computer algebra. After expanding the Hansen coefficients into Taylor series we take into consideration only terms up to determined degree N with respect to (e − e0 ) and obtain the approximating polynomial in the form: Xsn,m (e) =
N X
Ak (e − e0 )k ,
k=0
where the coefficients Ak are computed numerically for all values of n, m, s and e0 with the help of standard computer-algebraic functions for Taylor expansions. Our computer experiments, realized for different indices n, m, s and in the vicinity of different points e0 , show that the Taylor expansions always converge
42
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
and increasing the degree N of approximating polynomial we always improve also the accuracy of approximation. On the other hand, the construction of Taylor expansions in the neighbourhood of the point e0 6= 0 is very complicated from the point of view of necessary computer resources (large computer memory and computing time) and can be easily realized only for polynomials of low degrees (N 6 15). Let us present now some estimations of the accuracy of Taylor approximations of Hansen coefficients. We would like to point out that all these estimations were obtained only on the basis of numerical experiments, and the aim of the present paper is to help to choose the most reasonable way to construct the efficient polynomial approximations for Hansen coefficients, but not to give the strong analytical estimations on the numerical efficiency of the corresponding approximations. If e0 = 0, n < 0, m 6= | − n − 1|, then the Hansen coefficients Xsn,m (e) are monotonic functions of eccentricity, have a constant sign and increase with increasing the eccentricity e. The logarithm of the relative error of approximation D < 0, but beginning
from a determined boundary value of eccentricity it increases very rapidly up to zero with increasing the eccentricity. This boundary value of eccentricity depends on degree N of approximating polynomial. In particular, for N = 10 the fast increase of parameter D begins already in a small vicinity of the point e0 = 0, and the corresponding approximations are valid only for very small eccentricities. For N = 20 and number of digits 16 the parameter D does not increase and always D < −12 in the domain of eccentricities e < 0.1. In the case N = 50 the logarithm of error D < −12 by digits 16 already for eccentricities e < 0.4, for some values of indices n, m, s even for e < 0.5. It means that the similar approximations (at N = 50) represent the Hansen coefficients with a good precision already on a relative large interval of eccentricities. The behaviour of parameter D as a function of eccentricity for n = −5, m = 2, s = 1 is presented in Fig. 1 for N = 20 (curve a) and N = 50 (curve b). In the case e0 = 0, n < 0, m = |−n−1| the Hansen coefficients Xsn,m (e) are not monotonic functions of eccentricity and show the change of sign at determined
Fig. 1. Accuracy of Taylor approximation of Hansen coefficients as a function of eccentricity for the following cases: (a) e0 = 0, N = 20, n = −5, m = 2, s = 1; (b) e0 = 0, N = 50, n = −5, m = 2, s = 1; (c) e0 = 0, N = 50, n = −5, m = 4, s = 3; (d) e0 = 0.3, N = 15, n = −5, m = 4, s = 1.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
values of e. In the small vicinity of such points the approximation errors becomes very large and sometimes the parameter D becomes even positive. It means that Xe − Xa > Xe , and the approximations do not give reliable results. On the other hand, such large errors appear only in the domain where the parameter D begins to increase with increasing the eccentricity and these errors can be reduced if the degree N of the approximating polynomial is enough large. In particular, for N = 20 the parameter D < −12 in the domain e < 0.1; for N = 50 we have D < −12 in the domain e < 0.4 (see, as an example, Fig. 1, curve c). It means that increasing N we extend the domain of validity of Taylor expansions also in the case n < 0, m = | − n − 1|. The case e0 = 0, n > 0 is similar to the case e0 = 0, n < 0, m = | − n − 1|. The Hansen coefficients here change the sign at determined values of eccentricity and these domains can be approximated by Taylor expansions with low relative accuracy, but increasing N we extend the interval of validity of these approximations: for N = 20 the parameter D < −10 for all eccentricities e < 0.1, and for N = 50 we have D < −10 for all values of eccentricity e < 0.4. We also notice that the accuracy of approximations decreases with increasing the modulus of difference |m − s|, but these effects are very small in comparison with effects obtained by variation of degree of approximating polynomial N . Finally, some remarks concerning the case e0 6= 0. In this case the domain of validity of Taylor expansions becomes very small; it is only a very small vicinity of the point e0 . The behaviour of parameter D is similar to the case e0 = 0 in the domain e > e0 . But if e < e0 , then we observe a very fast increase in the parameter D up to positive values, that makes the Taylor approximations not valid in this region at all (see, for example, Fig. 1, curve d). Taking into consideration that the construction of Taylor expansions in the neighbourhood of the point e0 6= 0 up to high degrees N is very complicated we can recommend: (1) to use the approximation of Hansen coefficients by fragments of Taylor series in the neighbourhood of the point e0 = 0 (up to respective degree N ) in the case of small eccentricities;
43
(2) not to use the Taylor expansions and to choose other types of polynomial approximations of Hansen coefficients in the case of large eccentricities.
3. Approximation of Hansen coefficients with the help of Lagrange polynomials Quite generally, the interpolation of a function X(e) on the determined interval [a, b] by Lagrange polynomial can be made with the help of the following formula: X(e) =
N X i=0
X(ei )
N Y e − ek , ei − ek
k=0 k6=i
where N + 1 is the number of basic points for interpolation (corresponds to the degree of polynomial), ei = a + ih, ek = a + kh and h = (b − a)/N is the step of interpolation. X(ei ) is the value of the function X(e) in the point e = ei . It is very important to notice that by increasing the degree of interpolating polynomial N for Taylor expansions we always increase the accuracy of approximations. In the case of interpolation by Lagrange polynomials it is not so. Indeed, increasing N we increase the number of basic points for interpolation, but decrease the value of divisor ei − ek which leads to substantial loss in the accuracy. The problem of small divisors can always be solved by increasing the number of digits in the process of computation. On the other hand, for each number of digits there is a boundary degree Nmax of interpolating polynomial so that increasing N in the domain N > Nmax we will lose the accuracy of approximation. Therefore it is very efficient to use a large number of digits for construction of interpolating polynomial, although for computing the Hansen coefficients with the help of these interpolating polynomials it is enough to use the smaller number of digits. In all computer experiments described below we have used Digits = 32 for construction of the interpolating polynomial and Digits = 16 for computation of Hansen coefficients on the basis of this polynomial. All experiments were performed using the computer algebra system MAPLE V. Interpolation on large intervals of eccentricity (1e > 0.5) is not very efficient, if we would like to
44
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Fig. 2. Accuracy of interpolation of Hansen coefficients by Lagrange polynomials as a function of eccentricity for the following cases: (a) 0.0 < e < 0.9, N = 30, n = −5, m = 2, s = 1; (b) 0.1 < e < 0.3, N = 20, n = −5, m = 2, s = 1; (c) 0.0 < e < 0.1, N = 15, n = −5, m = 3, s = 1; (d) 0.0 < e < 0.1, N = 15, n = −5, m = 3, s = 3; (f) 0.4 < e < 0.5, N = 10, n = 5, m = 3, s = 1.
approximate Hansen coefficients with a high accuracy. For the interval [0.0; 0.9] the optimal degree of interpolating polynomial Nmax ≈ 30 for above indicated number of digits. On the Fig. 2 (curve a) we can see the behaviour of parameter D (introduced in the previous section) as a function of eccentricity for coefficients with n = −5, m = 2, s = 1. We can see that the maximal accuracy is reached in the center of approximated interval, and the maximal errors are on the borders, especially in the vicinity of the point e = 0, where the absolute values of the Hansen coefficients themselves are also very small. The similar behaviour of parameter D is observed also for other indices n, m, s on this interval of approximation. Interpolation on middle intervals of eccentricity (1e ≈ 0.2–0.4) is quite efficient. Here we have a possibility of reaching the maximal accuracy of interpolation. The optimal degree of approximating polynomial Nmax depends on the interval of interpolation, and, as an example, for the interval of eccentricity [0.1; 0.3] and Digits = 16 we have Nmax = 20. Curve b on the Fig. 2 shows the behaviour of parameter D for the indices: n = −5, m = 2, s = 1. For other indices the
behaviour of function D is similar. For the interval of eccentricity [0.5; 0.7] the optimal degree of polynomial Nmax ≈ 15 and the accuracy of interpolation decreases with increasing N > Nmax . The singularity appears in the small vicinity of the point e = 0, where the Hansen coefficients have small absolute values. In this domain all approximations represent the Hansen coefficients with very large relative errors, although the absolute errors here can be small because of small values of the Hansen coefficients themselves. Interpolation on small intervals of eccentricity (1e ≈ 0.1) gives also quite acceptable accuracy (Fig. 2, curves c and f). The exception is the case of Hansen coefficients Xsn,m with m = s, if the left border of interval of approximation e = 0. The accuracy of approximation of such coefficients is always very low (see, for example, Fig. 2, curve d). This effect takes place also for other indices (also if n < 0), but it does not appear if the left border of interpolation interval e 6= 0. The optimal degree of approximating polynomial in this case is Nmax ≈ 15. With increasing the absolute values of eccentricity Nmax decreases, but not very rapidly.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
On the basis of computer experiments we can conclude: (1) the maximal efficiency is reached by interpolations on middle intervals of eccentricity, the optimal degree of interpolating polynomial in this case is placed between Nmax ≈ 15 and Nmax ≈ 20 for digits = 32; (2) the interpolation in the neighbourhood of the point e = 0 has always large relative errors; (3) the construction of the interpolating polynomial should be made with a large number of digits in order to avoid losing accuracy because of small divisors; the forthcoming computation of Hansen coefficients can be performed with a standard (in particular, double) precision. 4. Chebyshev approximation of Hansen coefficients For constructing the Chebyshev approximations of the Hansen coefficients we can use the standard function chebyshev of the computer algebra system MAPLE V [8]. The accuracy of approximation is a
45
built-in parameter eps of this function. The optimal degree of the approximating polynomial is therefore determined here automatically. We notice that the parameter eps is defined here as the absolute accuracy eps = |Xe − Xa | and not equivalent to the parameter D introduced in Section 2 of the present paper. Here Xe is the “exact value” of the approximated function and Xa is the value obtained from approximating polynomial. Numerical experiments show that the Chebyshev approximation is able to provide in many cases a higher accuracy of approximations than the Taylor and Lagrange polynomials. The relative error D can be very large here only in the domain of very small values of the Hansen coefficients. On Fig. 3 (curves a, b, c) we can see the typical behaviour of parameter D for the Chebyshev approximation on the interval 0.0 < e < 0.9 with the number of digits 16 and parameter eps = 10−16 . The errors of approximation are very large for small eccentricities, where the absolute values of the Hansen coefficients are very small. The behaviour of the Hansen coefficients themselves with the same indices is presented on Fig. 4. The de-
Fig. 3. Accuracy of Chebyshev approximation of Hansen coefficients as a function of eccentricity for the following cases: (a) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = −2; (b) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = 1; (c) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = 3; (d) 0.0 < e < 0.9, eps = 10−16 , n = 2; m = 2; s = 1.
46
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Fig. 4. Logarithm of absolute values of Hansen coefficients on the interval 0.0 < e < 0.9 for the following cases: (a) n = −5, m = 4, s = −2; (b) n = −5, m = 4, s = 1; (c) n = −5, m = 4, s = 3.
gree of the approximating polynomial for the interval 0.0 < e < 0.9 and eps = 10−16 evaluated by means of MAPLE-function chebyshev is between N = 40 and N = 50 for all cases which were under consideration. If we decrease the interval of approximation then the degree of the approximating polynomial is also decreased but slowly. For the interval of approximation 1e = 0.1 this degree is between N = 10 and N = 20, if digits = 16 and eps = 10−16 . It means that it is not reasonable to use approximations on many short intervals of eccentricity instead of one approximation on a large interval. The jumps (as in the case of curve c on the Fig. 3 in the vicinity of the point e = 0.66) appear in the neighbourhood of the points where the Hansen coefficients change sign. In the vicinity of such points the absolute values of the Hansen coefficients become small and the parameter eps is also small, but the relative error D can be very large. In Fig. 3 we also show the relative errors of the Chebyshev approximation of the Hansen coefficients in the case of positive index n (curve d). It is also very important to notice that construction of Chebyshev approximations for large positive values of index n is very complicated and requires a large computer memory and computing time.
5. Conclusion Numerical experiments show that the usage of approximating polynomials gives a possibility of computing the Hansen coefficients some hundreds of times faster than using direct methods. Moreover, such an approach allows us also to utilize computer memory very economically. It is very important for increasing the efficiency of theories of motion of celestial bodies containing the expansions with Hansen coefficients. Concerning the efficiency of different approximation schemes themselves, we can make the following conclusions: Lagrange interpolations are very convenient, if we would like to construct the approximating polynomial quickly, but interpolate the Hansen coefficients with not a very high accuracy. These interpolations are especially efficient for middle and short intervals of eccentricity (1e ≈ 0.1–0.3) and in the domain of moderately large eccentricities (0.2 < e < 0.5). It is also very important to evaluate the optimal degree of interpolating polynomial beforehand in order to avoid the loss of accuracy because of small divisors for polynomials of high degrees.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Chebyshev approximations are always able to provide very small absolute errors of approximation, but the relative errors can be very large in the domain of small eccentricities. In many practical problems such large relative errors do not play any important role at all, because they are connected with effects, the influence of which is not very important (for small eccentricities there are also always the coefficients with not very small absolute values and these coefficients are approximated with a high relative accuracy). On the other hand, in order to study the determined types of perturbations in the motion of a celestial body it can sometimes be necessary to evaluate the Hansen coefficients with the determined indices. In this case we should be sure that the approximations give for these coefficients small relative and absolute errors simultaneously. Although Chebyshev approximations allow us to reach a higher accuracy than the Lagrange interpolations, we should notice that construction of Chebyshev approximations can be very complicated (especially for large positive indices n) and require a large computer memory and computing time. Taylor approximations are efficient only if they are constructed in the neighbourhood of the point e0 = 0. It is the only type of approximations which can provide the small relative errors in the domain of small eccentricities. Increasing the degree of approximating polynomial always improves the accuracy of approximations. The optimal interval for Taylor approximations: e < 0.5. The optimal degree of approximating polynomial (for the case 0.0 < e < 0.5) is about N = 50. Let us make now some remarks concerning the behaviour of approximations in the points where the Hansen coefficients change sign. We can see that in the neighbourhood of such points the relative errors D become large sometimes (curve c on the Figs. 1 and 3). On the other hand, the absolute values of the Hansen coefficients in these points become small (curve c on the Fig. 4) and in our numerical experiments we can see that the absolute errors eps = |Xe − Xa | also remain small. Our estimations show that ∀n, m if we P compute the full sum s Xsn,m (e) exp(j sM), then the effect of changing the sign of Hansen coefficients with the determined index s does not introduce any additional errors because the corresponding coefficients themselves are very small and can be eliminated from the consideration at all. But if we would like to study
47
the influence of the determined harmonical term of the disturbing function on the motion of a celestial body and the corresponding term contains the Hansen coefficient passing through zero, then the use of approximating polynomials can be not reliable, and the errors of approximations should be evaluated beforehand in order to avoid the loss of accuracy. The same remark concerns also the case of very small eccentricities, where only the Taylor approximations in the vicinity of e0 = 0 give both small relative and small absolute errors. Finally, we would like to point out again that all conclusions made above are based on many computer experiments, and it is not the aim of this paper to give the strong analytical estimations of efficiency of different approximating schemes. The results presented in the paper should help one choose the most reasonable way for constructing the polynomial approximations of Hansen coefficients for concrete problems of celestial mechanics and help to increase the efficiency of corresponding theories of motion of celestial bodies. Acknowledgements The author would like to express the great gratitude to Dr. Irina Tupikova for her help in this research and many useful discussions. The author is very grateful also to Prof. M. Soffel for his kind attention to this work. Some discussions with Dr. S. Klioner were also very useful. References [1] W.M. Kaula, Theory of Satellite Geodesy (Blaisdell, 1966). [2] E. Gaposhkin, Satellite dynamics, SAO Special Report, 1973, p. 353. [3] I.V. Tupikova, A.A. Vakhidov, M. Soffel, Explicit semianalytical theory of asteroid motion, in: Proc. of IAU Coll., No. 172, J. Henrard, S. Ferraz-Mello (Eds.) (Kluwer Acad. Publ., 1999). [4] G.E.O. Giacaglia, A note on Hansen’s coefficients in satellite theory, Celest. Mech. 14 (1976) 515. [5] S. Hughes, The computation of tables of Hansen coefficients, Celest. Mech. 25 (1981) 101. [6] C.C. Goad, An efficient algorithm for the evaluation of inclination and eccentricity functions, Manuscripta Geodaetica 12 (1987) 11. [7] S.A. Klioner, A.A. Vakhidov, N.N. Vasiliev, Numerical computation of Hansen-like expansions, Celest. Mech. and Dynam. Astron. 68 (1998) 257.
48
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
[8] B.W. Char, K.O. Geddes, G.H. Gonnet, B.L. Leong, M.B. Monagan, S.M. Watt, MAPLE V Library Reference Manual (Springer, New York, 1993).
[9] V.A. Brumberg, Analytical Techniques of Celestial Mechanics (Springer, Berlin-Heidelberg, 1995).
Construction of polynomial approximation for Hansen coefficients Akmal A.Vakhidov
1
Lohrmann Observatory, Technical University Dresden, Mommsenstrasse 13, Dresden 01062, Germany Received 18 March 1999; received in revised form 17 May 1999
Abstract The problem of constructing efficient approximations of Hansen coefficients using polynomials in terms of the eccentricity is considered. The properties of the following approximating schemes are studied in detail: approximation by fragments of Taylor series, interpolation by Lagrange polynomials, Chebyshev approximation. The accuracy of all these approximating schemes is investigated for different values of eccentricities and for different indices of the Hansen coefficients. 2000 Elsevier Science B.V. All rights reserved. PACS: 95.10 Ce – Celestial Mechanics; 02.30 Mv – Approximations and Expansions Keywords: Hansen coefficients; Polynomial approximation; Computer algebra
1. Introduction If we construct analytical or semianalytical theories of motion of celestial bodies, then in many cases it is very convenient to use the expansions of coordinates into trigonometric series in multiples of the mean anomaly depending on the time linearly. Such an approach gives us the possibility of solving analytically the equations of motion using the method of consequent approximations (see, for example, [1,2]) or to realize the averaging procedure based on von Zeipel or Lie-series methods [3]. Quite generally, the expansion of coordinates in multiples of the mean anomaly in the case of elliptic motion is given by the formula: n ∞ X r exp(j mv) = Xsn,m (e) exp(j sM), a s=−∞ 1 E-mail: [email protected].
where (r, v) are the polar coordinates of the celestial body, a is the semi-major axis of the elliptic orbit, M is√the mean anomaly, n, m, s are some integers, j = −1. The coefficients Xsn,m (called as Hansen coefficients) depend on the eccentricity e of the orbit. The problem of efficient computation of Hansen coefficients was investigated in detail in many papers (see also [4–7]), but this problem is really very complicated and to date we have no standard algorithms for computing these coefficients efficiently for all possible values of indices n, m, s and eccentricity e. A great problem in using the expansions with Hansen coefficients in theories of motion of celestial bodies is connected with the necessity to recompute these coefficients many times for different (very often a little different) values of eccentricity. Indeed, when using a semianalytical theory of motion or theory presented via orbit elements, it is necessary to recompute all Hansen coefficients at each step of the numerical integration of the equations of motion. It
0010-4655/00/$ – see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 9 9 ) 0 0 4 3 1 - 2
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
leads to many inconveniences from the point of view of economy of computing time and decreases the efficiency of the theory itself. Therefore sometimes it seems to be much more efficient to construct some approximations, which give the possibility of computing the Hansen coefficients on the determined interval of values of eccentricity. Such approximations could be obtained before the usage of motion theory and allow one to evaluate the coefficients using compact approximating formulae instead of more complicated direct methods. The accuracy reached in approximations should satisfy to the accuracy of the corresponding analytical or semianalytical theory. For many practical problems of celestial mechanics in theories of the first order the accuracy of computing the Hansen coefficients should be equal to 10−4 . The errors then will be smaller than the terms of the second order, which are not taken into consideration. Such an accuracy is enough for theories of the first order for such problems as the motion of minor planets or the motion of artificial Earth satellites, where the small parameter of the problem is of order 10−3 . If we want to take into consideration the perturbations of the second order then the accuracy of computing the disturbing function (and Hansen coefficients, respectively) should be 10−7 . But of course, the desired accuracy of computing the Hansen coefficients depend on the concrete problem and is different for different problems of celestial mechanics. Therefore it is very important to construct the approximations in a way that could incorporate the simplicity of computing the Hansen coefficients and the desired accuracy simultaneously. In the present paper we investigate three approximating schemes for computing the Hansen coefficients: (1) approximation by fragments of Taylor series; (2) interpolation with the help of Lagrange polynomials; (3) Chebyshev approximation. We present the results of computer experiments to estimate the accuracy of these approximating schemes and advantages and inconveniences of usage of each of them. For numerical experiments we consider the Hansen coefficients both with n < 0 (the case of perturbations from internal body) and with n > 0 (perturbations from external body), different values of indices m and s, and different intervals of eccentricity for approximations. All computer experiments described below were realized on Sun Sparc workstations in the computer alge-
41
bra system MAPLE V (Release 4) [8] with different numbers of digits. The software for constructing the described approximations is available from the author upon request.
2. Approximation of Hansen coefficients by fragments of Taylor series For estimating the accuracy of polynomial approximations of Hansen coefficients we will consider the behaviour of the logarithm of the relative error defined as Xe − Xa , D = log Xe where Xe is the “exact value” of Hansen coefficient and Xa is the value obtained with the help of the approximating polynomial. By the words “exact value” we understand here not the true value computed with absolute accuracy, but the value obtained using direct methods for computing the Hansen coefficients. These “exact values” Xe can be computed, in particular, on the basis of the formulae of Hill’s method [9]. This method is based on the representation of Hansen coefficients via hypergeometric functions and Bessel functions which can be evaluated with arbitrary precision as the standard built-in functions of the computer algebra system MAPLE V [8]. Hill’s method allows also the computation of the Hansen coefficients in symbolic form for the forthcoming expansion of the corresponding analytical expressions of Hansen coefficients into Taylor series. This procedure can also be realized by means of computer algebra. After expanding the Hansen coefficients into Taylor series we take into consideration only terms up to determined degree N with respect to (e − e0 ) and obtain the approximating polynomial in the form: Xsn,m (e) =
N X
Ak (e − e0 )k ,
k=0
where the coefficients Ak are computed numerically for all values of n, m, s and e0 with the help of standard computer-algebraic functions for Taylor expansions. Our computer experiments, realized for different indices n, m, s and in the vicinity of different points e0 , show that the Taylor expansions always converge
42
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
and increasing the degree N of approximating polynomial we always improve also the accuracy of approximation. On the other hand, the construction of Taylor expansions in the neighbourhood of the point e0 6= 0 is very complicated from the point of view of necessary computer resources (large computer memory and computing time) and can be easily realized only for polynomials of low degrees (N 6 15). Let us present now some estimations of the accuracy of Taylor approximations of Hansen coefficients. We would like to point out that all these estimations were obtained only on the basis of numerical experiments, and the aim of the present paper is to help to choose the most reasonable way to construct the efficient polynomial approximations for Hansen coefficients, but not to give the strong analytical estimations on the numerical efficiency of the corresponding approximations. If e0 = 0, n < 0, m 6= | − n − 1|, then the Hansen coefficients Xsn,m (e) are monotonic functions of eccentricity, have a constant sign and increase with increasing the eccentricity e. The logarithm of the relative error of approximation D < 0, but beginning
from a determined boundary value of eccentricity it increases very rapidly up to zero with increasing the eccentricity. This boundary value of eccentricity depends on degree N of approximating polynomial. In particular, for N = 10 the fast increase of parameter D begins already in a small vicinity of the point e0 = 0, and the corresponding approximations are valid only for very small eccentricities. For N = 20 and number of digits 16 the parameter D does not increase and always D < −12 in the domain of eccentricities e < 0.1. In the case N = 50 the logarithm of error D < −12 by digits 16 already for eccentricities e < 0.4, for some values of indices n, m, s even for e < 0.5. It means that the similar approximations (at N = 50) represent the Hansen coefficients with a good precision already on a relative large interval of eccentricities. The behaviour of parameter D as a function of eccentricity for n = −5, m = 2, s = 1 is presented in Fig. 1 for N = 20 (curve a) and N = 50 (curve b). In the case e0 = 0, n < 0, m = |−n−1| the Hansen coefficients Xsn,m (e) are not monotonic functions of eccentricity and show the change of sign at determined
Fig. 1. Accuracy of Taylor approximation of Hansen coefficients as a function of eccentricity for the following cases: (a) e0 = 0, N = 20, n = −5, m = 2, s = 1; (b) e0 = 0, N = 50, n = −5, m = 2, s = 1; (c) e0 = 0, N = 50, n = −5, m = 4, s = 3; (d) e0 = 0.3, N = 15, n = −5, m = 4, s = 1.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
values of e. In the small vicinity of such points the approximation errors becomes very large and sometimes the parameter D becomes even positive. It means that Xe − Xa > Xe , and the approximations do not give reliable results. On the other hand, such large errors appear only in the domain where the parameter D begins to increase with increasing the eccentricity and these errors can be reduced if the degree N of the approximating polynomial is enough large. In particular, for N = 20 the parameter D < −12 in the domain e < 0.1; for N = 50 we have D < −12 in the domain e < 0.4 (see, as an example, Fig. 1, curve c). It means that increasing N we extend the domain of validity of Taylor expansions also in the case n < 0, m = | − n − 1|. The case e0 = 0, n > 0 is similar to the case e0 = 0, n < 0, m = | − n − 1|. The Hansen coefficients here change the sign at determined values of eccentricity and these domains can be approximated by Taylor expansions with low relative accuracy, but increasing N we extend the interval of validity of these approximations: for N = 20 the parameter D < −10 for all eccentricities e < 0.1, and for N = 50 we have D < −10 for all values of eccentricity e < 0.4. We also notice that the accuracy of approximations decreases with increasing the modulus of difference |m − s|, but these effects are very small in comparison with effects obtained by variation of degree of approximating polynomial N . Finally, some remarks concerning the case e0 6= 0. In this case the domain of validity of Taylor expansions becomes very small; it is only a very small vicinity of the point e0 . The behaviour of parameter D is similar to the case e0 = 0 in the domain e > e0 . But if e < e0 , then we observe a very fast increase in the parameter D up to positive values, that makes the Taylor approximations not valid in this region at all (see, for example, Fig. 1, curve d). Taking into consideration that the construction of Taylor expansions in the neighbourhood of the point e0 6= 0 up to high degrees N is very complicated we can recommend: (1) to use the approximation of Hansen coefficients by fragments of Taylor series in the neighbourhood of the point e0 = 0 (up to respective degree N ) in the case of small eccentricities;
43
(2) not to use the Taylor expansions and to choose other types of polynomial approximations of Hansen coefficients in the case of large eccentricities.
3. Approximation of Hansen coefficients with the help of Lagrange polynomials Quite generally, the interpolation of a function X(e) on the determined interval [a, b] by Lagrange polynomial can be made with the help of the following formula: X(e) =
N X i=0
X(ei )
N Y e − ek , ei − ek
k=0 k6=i
where N + 1 is the number of basic points for interpolation (corresponds to the degree of polynomial), ei = a + ih, ek = a + kh and h = (b − a)/N is the step of interpolation. X(ei ) is the value of the function X(e) in the point e = ei . It is very important to notice that by increasing the degree of interpolating polynomial N for Taylor expansions we always increase the accuracy of approximations. In the case of interpolation by Lagrange polynomials it is not so. Indeed, increasing N we increase the number of basic points for interpolation, but decrease the value of divisor ei − ek which leads to substantial loss in the accuracy. The problem of small divisors can always be solved by increasing the number of digits in the process of computation. On the other hand, for each number of digits there is a boundary degree Nmax of interpolating polynomial so that increasing N in the domain N > Nmax we will lose the accuracy of approximation. Therefore it is very efficient to use a large number of digits for construction of interpolating polynomial, although for computing the Hansen coefficients with the help of these interpolating polynomials it is enough to use the smaller number of digits. In all computer experiments described below we have used Digits = 32 for construction of the interpolating polynomial and Digits = 16 for computation of Hansen coefficients on the basis of this polynomial. All experiments were performed using the computer algebra system MAPLE V. Interpolation on large intervals of eccentricity (1e > 0.5) is not very efficient, if we would like to
44
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Fig. 2. Accuracy of interpolation of Hansen coefficients by Lagrange polynomials as a function of eccentricity for the following cases: (a) 0.0 < e < 0.9, N = 30, n = −5, m = 2, s = 1; (b) 0.1 < e < 0.3, N = 20, n = −5, m = 2, s = 1; (c) 0.0 < e < 0.1, N = 15, n = −5, m = 3, s = 1; (d) 0.0 < e < 0.1, N = 15, n = −5, m = 3, s = 3; (f) 0.4 < e < 0.5, N = 10, n = 5, m = 3, s = 1.
approximate Hansen coefficients with a high accuracy. For the interval [0.0; 0.9] the optimal degree of interpolating polynomial Nmax ≈ 30 for above indicated number of digits. On the Fig. 2 (curve a) we can see the behaviour of parameter D (introduced in the previous section) as a function of eccentricity for coefficients with n = −5, m = 2, s = 1. We can see that the maximal accuracy is reached in the center of approximated interval, and the maximal errors are on the borders, especially in the vicinity of the point e = 0, where the absolute values of the Hansen coefficients themselves are also very small. The similar behaviour of parameter D is observed also for other indices n, m, s on this interval of approximation. Interpolation on middle intervals of eccentricity (1e ≈ 0.2–0.4) is quite efficient. Here we have a possibility of reaching the maximal accuracy of interpolation. The optimal degree of approximating polynomial Nmax depends on the interval of interpolation, and, as an example, for the interval of eccentricity [0.1; 0.3] and Digits = 16 we have Nmax = 20. Curve b on the Fig. 2 shows the behaviour of parameter D for the indices: n = −5, m = 2, s = 1. For other indices the
behaviour of function D is similar. For the interval of eccentricity [0.5; 0.7] the optimal degree of polynomial Nmax ≈ 15 and the accuracy of interpolation decreases with increasing N > Nmax . The singularity appears in the small vicinity of the point e = 0, where the Hansen coefficients have small absolute values. In this domain all approximations represent the Hansen coefficients with very large relative errors, although the absolute errors here can be small because of small values of the Hansen coefficients themselves. Interpolation on small intervals of eccentricity (1e ≈ 0.1) gives also quite acceptable accuracy (Fig. 2, curves c and f). The exception is the case of Hansen coefficients Xsn,m with m = s, if the left border of interval of approximation e = 0. The accuracy of approximation of such coefficients is always very low (see, for example, Fig. 2, curve d). This effect takes place also for other indices (also if n < 0), but it does not appear if the left border of interpolation interval e 6= 0. The optimal degree of approximating polynomial in this case is Nmax ≈ 15. With increasing the absolute values of eccentricity Nmax decreases, but not very rapidly.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
On the basis of computer experiments we can conclude: (1) the maximal efficiency is reached by interpolations on middle intervals of eccentricity, the optimal degree of interpolating polynomial in this case is placed between Nmax ≈ 15 and Nmax ≈ 20 for digits = 32; (2) the interpolation in the neighbourhood of the point e = 0 has always large relative errors; (3) the construction of the interpolating polynomial should be made with a large number of digits in order to avoid losing accuracy because of small divisors; the forthcoming computation of Hansen coefficients can be performed with a standard (in particular, double) precision. 4. Chebyshev approximation of Hansen coefficients For constructing the Chebyshev approximations of the Hansen coefficients we can use the standard function chebyshev of the computer algebra system MAPLE V [8]. The accuracy of approximation is a
45
built-in parameter eps of this function. The optimal degree of the approximating polynomial is therefore determined here automatically. We notice that the parameter eps is defined here as the absolute accuracy eps = |Xe − Xa | and not equivalent to the parameter D introduced in Section 2 of the present paper. Here Xe is the “exact value” of the approximated function and Xa is the value obtained from approximating polynomial. Numerical experiments show that the Chebyshev approximation is able to provide in many cases a higher accuracy of approximations than the Taylor and Lagrange polynomials. The relative error D can be very large here only in the domain of very small values of the Hansen coefficients. On Fig. 3 (curves a, b, c) we can see the typical behaviour of parameter D for the Chebyshev approximation on the interval 0.0 < e < 0.9 with the number of digits 16 and parameter eps = 10−16 . The errors of approximation are very large for small eccentricities, where the absolute values of the Hansen coefficients are very small. The behaviour of the Hansen coefficients themselves with the same indices is presented on Fig. 4. The de-
Fig. 3. Accuracy of Chebyshev approximation of Hansen coefficients as a function of eccentricity for the following cases: (a) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = −2; (b) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = 1; (c) 0.0 < e < 0.9, eps = 10−16 , n = −5, m = 4, s = 3; (d) 0.0 < e < 0.9, eps = 10−16 , n = 2; m = 2; s = 1.
46
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Fig. 4. Logarithm of absolute values of Hansen coefficients on the interval 0.0 < e < 0.9 for the following cases: (a) n = −5, m = 4, s = −2; (b) n = −5, m = 4, s = 1; (c) n = −5, m = 4, s = 3.
gree of the approximating polynomial for the interval 0.0 < e < 0.9 and eps = 10−16 evaluated by means of MAPLE-function chebyshev is between N = 40 and N = 50 for all cases which were under consideration. If we decrease the interval of approximation then the degree of the approximating polynomial is also decreased but slowly. For the interval of approximation 1e = 0.1 this degree is between N = 10 and N = 20, if digits = 16 and eps = 10−16 . It means that it is not reasonable to use approximations on many short intervals of eccentricity instead of one approximation on a large interval. The jumps (as in the case of curve c on the Fig. 3 in the vicinity of the point e = 0.66) appear in the neighbourhood of the points where the Hansen coefficients change sign. In the vicinity of such points the absolute values of the Hansen coefficients become small and the parameter eps is also small, but the relative error D can be very large. In Fig. 3 we also show the relative errors of the Chebyshev approximation of the Hansen coefficients in the case of positive index n (curve d). It is also very important to notice that construction of Chebyshev approximations for large positive values of index n is very complicated and requires a large computer memory and computing time.
5. Conclusion Numerical experiments show that the usage of approximating polynomials gives a possibility of computing the Hansen coefficients some hundreds of times faster than using direct methods. Moreover, such an approach allows us also to utilize computer memory very economically. It is very important for increasing the efficiency of theories of motion of celestial bodies containing the expansions with Hansen coefficients. Concerning the efficiency of different approximation schemes themselves, we can make the following conclusions: Lagrange interpolations are very convenient, if we would like to construct the approximating polynomial quickly, but interpolate the Hansen coefficients with not a very high accuracy. These interpolations are especially efficient for middle and short intervals of eccentricity (1e ≈ 0.1–0.3) and in the domain of moderately large eccentricities (0.2 < e < 0.5). It is also very important to evaluate the optimal degree of interpolating polynomial beforehand in order to avoid the loss of accuracy because of small divisors for polynomials of high degrees.
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
Chebyshev approximations are always able to provide very small absolute errors of approximation, but the relative errors can be very large in the domain of small eccentricities. In many practical problems such large relative errors do not play any important role at all, because they are connected with effects, the influence of which is not very important (for small eccentricities there are also always the coefficients with not very small absolute values and these coefficients are approximated with a high relative accuracy). On the other hand, in order to study the determined types of perturbations in the motion of a celestial body it can sometimes be necessary to evaluate the Hansen coefficients with the determined indices. In this case we should be sure that the approximations give for these coefficients small relative and absolute errors simultaneously. Although Chebyshev approximations allow us to reach a higher accuracy than the Lagrange interpolations, we should notice that construction of Chebyshev approximations can be very complicated (especially for large positive indices n) and require a large computer memory and computing time. Taylor approximations are efficient only if they are constructed in the neighbourhood of the point e0 = 0. It is the only type of approximations which can provide the small relative errors in the domain of small eccentricities. Increasing the degree of approximating polynomial always improves the accuracy of approximations. The optimal interval for Taylor approximations: e < 0.5. The optimal degree of approximating polynomial (for the case 0.0 < e < 0.5) is about N = 50. Let us make now some remarks concerning the behaviour of approximations in the points where the Hansen coefficients change sign. We can see that in the neighbourhood of such points the relative errors D become large sometimes (curve c on the Figs. 1 and 3). On the other hand, the absolute values of the Hansen coefficients in these points become small (curve c on the Fig. 4) and in our numerical experiments we can see that the absolute errors eps = |Xe − Xa | also remain small. Our estimations show that ∀n, m if we P compute the full sum s Xsn,m (e) exp(j sM), then the effect of changing the sign of Hansen coefficients with the determined index s does not introduce any additional errors because the corresponding coefficients themselves are very small and can be eliminated from the consideration at all. But if we would like to study
47
the influence of the determined harmonical term of the disturbing function on the motion of a celestial body and the corresponding term contains the Hansen coefficient passing through zero, then the use of approximating polynomials can be not reliable, and the errors of approximations should be evaluated beforehand in order to avoid the loss of accuracy. The same remark concerns also the case of very small eccentricities, where only the Taylor approximations in the vicinity of e0 = 0 give both small relative and small absolute errors. Finally, we would like to point out again that all conclusions made above are based on many computer experiments, and it is not the aim of this paper to give the strong analytical estimations of efficiency of different approximating schemes. The results presented in the paper should help one choose the most reasonable way for constructing the polynomial approximations of Hansen coefficients for concrete problems of celestial mechanics and help to increase the efficiency of corresponding theories of motion of celestial bodies. Acknowledgements The author would like to express the great gratitude to Dr. Irina Tupikova for her help in this research and many useful discussions. The author is very grateful also to Prof. M. Soffel for his kind attention to this work. Some discussions with Dr. S. Klioner were also very useful. References [1] W.M. Kaula, Theory of Satellite Geodesy (Blaisdell, 1966). [2] E. Gaposhkin, Satellite dynamics, SAO Special Report, 1973, p. 353. [3] I.V. Tupikova, A.A. Vakhidov, M. Soffel, Explicit semianalytical theory of asteroid motion, in: Proc. of IAU Coll., No. 172, J. Henrard, S. Ferraz-Mello (Eds.) (Kluwer Acad. Publ., 1999). [4] G.E.O. Giacaglia, A note on Hansen’s coefficients in satellite theory, Celest. Mech. 14 (1976) 515. [5] S. Hughes, The computation of tables of Hansen coefficients, Celest. Mech. 25 (1981) 101. [6] C.C. Goad, An efficient algorithm for the evaluation of inclination and eccentricity functions, Manuscripta Geodaetica 12 (1987) 11. [7] S.A. Klioner, A.A. Vakhidov, N.N. Vasiliev, Numerical computation of Hansen-like expansions, Celest. Mech. and Dynam. Astron. 68 (1998) 257.
48
A.A. Vakhidov / Computer Physics Communications 124 (2000) 40–48
[8] B.W. Char, K.O. Geddes, G.H. Gonnet, B.L. Leong, M.B. Monagan, S.M. Watt, MAPLE V Library Reference Manual (Springer, New York, 1993).
[9] V.A. Brumberg, Analytical Techniques of Celestial Mechanics (Springer, Berlin-Heidelberg, 1995).