CC BY
21
DOI: 10.17516/1997-1397-2020-13-6-781-791 УДК 519.6
Accuracy of Symmetric Multi-Step Methods for the Numerical Modelling of Satellite Motion
Evgenia D. Karepova*
Institute of computational modeling of SB RAS Krasnoyarsk, Russian Federation
Iliya R. Adaev^
Institute of computational modeling of SB RAS Krasnoyarsk, Russian Federation Siberian Federal University Krasnoyarsk, Russian Federation
Yury V. Shan'ko*
Institute of computational modeling of SB RAS Krasnoyarsk, Russian Federation
Received 10.08.2020, received in revised form 08.09.2020, accepted 07.10.2020 Abstract. Stability of high-order linear multistep Stormer-Cowell and symmetric methods are discussed in detail in this paper. Efficient algorithms for obtaining intervals of absolute stability and periodicity are given for these methods. To demonstrate the accuracy of numerical integration of the orbit over an interval about one year two model problems are considered. First problem is the 3D Kepler problem. Second one is a specially designed 3D model problem that has the exact solution and simulates the Earth-Moon-satellite system.
Keywords: linear multistep method, symmetric method, Stormer-Cowell method, PECE scheme, orbit.
Citation: E.D. Karepova, I.R. Adaev, Y.V. Shan'ko, Accuracy of the Symmetric Multi-Step Methods for the Numerical Modelling of Satellite Motion, J. Sib. Fed. Univ. Math. Phys., 2020, 13(6), 781-791. DOI: 10.17516/1997-1397-2020-13-6-781-791.
Introduction
Accuracy of the numerical integration of a satellite motion still remains one of the top problems associated with Global Navigation Satellite Systems. A review of the approaches used by Analysis centres of International GNSS Service [1] shows that the basic techniques of the numerical integration of a satellite orbit are the Adams-Bashforth/Moulton PECE-algorithms, the nonlinear Everhart's procedure [2] and collocation methods [3,4]. However, a linear multi-step symmetric methods shows considerable promise [5] for near-circular orbits that are typical for navigation satellites.
The theory of multi-step methods, including the Adams family which are traditional for the numerical integration of the motion of celestial objects, are widely discussed in many textbooks
*e.d.karepova<8icm.krasn.ru https://orcid.org/0000-0002-6515-2932 [email protected] https://orcid.org/0000-0002-5670-3747
[email protected] https://orcid.org/0000-0003-2796-4363
© Siberian Federal University. All rights reserved
on numerical methods [6,7,9-12]. The Stormer-Cowell methods were developed and successfully used since the early 20th century. However, in 2016 an interesting result concerning instability for small step size of some Stormer-Cowell methods was presented by N0rsett and Asheim [13]. The general theory of the symmetric multi-step methods was developed by Lambert and Watson [14]. The symmetric methods of high order were discussed in relation to the numerical integration of planetary orbits over a long period of time.
The orbital motion is described by the system of second order ordinary differential equations (ODE). It is generally agreed that is better to solve numerically the second order ODE rather than equivalent system of two first order equations [6,15]. We also confirm this in our numerical experiments.
In this paper, we discuss the accuracy and stability of high-order explicit symmetric multistep methods and their advantage over the Stormer-Cowell methods with/without "predict -evaluate - correct evaluate" (PECE) mode. We propose an efficient way to calculate intervals of absolute stability and periodicity for any linear multi-step method.
To study stability and periodicity we used the general-purpose computer algebra system REDUCE over the complex field with an accuracy of 40 significant digits. Numerical algorithms were implemented in C++ using the library quadmath for quadruple precision calculations.
1. Linear multistep methods
On the discrete point set {tn : tn = t0 + nh, h > 0,n = 0,1,... }, we consider the k-step linear multistep method
k k
E aj xn+j = h2 E j fn+j, k > 2 (1)
j=o j=0
for the numerical solution of the special second-order initial value problem
x" = f (t,x), x(to) = xo, x' (to ) = X. (2)
Here xn is the approximation of the exact solution x(tn) e R and fn = f (tn,xn). Method (1) is characterized by polynomials p(£) and a(£), where
kk p(£) = E aj, a(£) = E j, £ e C.
j=0 j=o
k
We suppose that p and a have no common factors, ak = 1, |a0| + \j0\ = 0, and \Pj I = 0.
j=o
If j3k =0 the method is explicit, otherwise it is implicit. For method (1) to be consistent, it is necessary and sufficient that p(1) = p'(1) = 0 and p''(1) = 2a(1). Method (1) has the order p if for all sufficiently smooth test functions z(t) kk Y^aj z(t + jh) - h2J2 j z ''(t + jh) = Cp+-2hp+2z(p+2'>(t) + O(hp+3). j=0 j=0
We assume that if the Cauchy problem (2) is solved with the use of method (1) the accuracy of first starting values xn, n = 0,... ,k — 1 is at least not less than the order of the method.
All Stormer-Cowell methods have p(£) = £k — 2£k-1 + £k-2. Method (1) is symmetric if a.j = ak-j, jjj = jk-j, j = 0,... ,k. A symmetric method has only even order [7]. We study higher order methods, namely, from 6th to 12th order Stoormer-Cowell methods and even order
symmetric methods. Coefficients a.j and j for these methods are presented in [13] and [5,14], respectively.
These methods are consistent and zero-stable. Hence they are convergent [6,14] and polynomial p has the root of multiplicity two at +1. Let us denote the roots of p by £s, s = 1,... ,k, where = £2 = 1 are the principal roots and the remaining k — 2 roots are spurious. All spurious roots of any Stormer-Cowell method are zero.
We demonstrate main differences between symmetric and Stormer-Cowell methods by the example of the harmonic oscillator equation
x" = —X2x, x(to) = xo, x'(to) = x, X G R. (3)
That has general solution x(t) = A cos Xt + B sin Xt with period T = 2n/X.
Using method (1) to solve (3), we obtain the difference equation
k
a+h 2 j xn+j=0 (4)
j=o
with general solution
k
xn = Dirn + D2T?[ + Y, Dsrn (5)
s=3
Here H = Xh, Ds G C are constant. Let us assume that all the roots rs, s=1,... ,k of the stability polynomial
n(r; H2)= p(r) + HV(t) (6)
are distinct. Since the roots of the polynomial are continuous functions of its coefficients, rs are perturbation of when H2 > 0. Thus, the numerical solution of (3) xn may be represented by the sum of the component (xn)P = D1rrn + D2rn associated with the perturbation of the principal roots and (xn)S that arises from perturbation of spurious roots.
Absolute stability of the Stormer and Cowell methods. Root-locus curves for some Stormer and Cowell methods are shown in Fig. 1 (a-i). They are constructed by the "boundary locus" method [12] which gives a general shape of the boundary |r| = | exp(«p)| = 1 of the open stability region in the complex plane H2. The stability region is always at the left of the curve when we move along the curve as p increases from 0 to 2n. For example, there is no stability region for the 10th order Stormer's method (Fig. 1 c). Moreover, the stability region near the interval of absolute stability is shown in more detail in Figs. 1 (b,f, h) for methods that are used in our numerical experiments. Let us note that in the general case X G C the stability region is determined, while for the harmonic oscillator we obtain the stability interval on the real axis.
In order to determine the stability interval more accurately the Routh-Hurwitz criterion [16] can be used. In this case, a transformation of the region |r| ^ 1 into the region Re(z) ^ 0 is required. There are the Schur-Cohn [12] and the Jury [17-19] criteria that test the strong stability of n(r; H2) directly. The Schur-Cohn and the Jury criteria are convenient for program implementation and they are easily tested for a given H2.
According to the Jury criterion, the problem of determining the set of all values of H2 that all roots of n(r; H2) are inside of the unit circle, is reduced to solving the system of k inequalities, where k is the degree of n(r; H2). The left-hand side of each inequality is the ratio of polynomials in H2 and the right-hand side is zero. The polynomial coefficients are obtained from the coefficients of n(r; H2). Even for small k the system of the inequalities can be analysed analytically only in some cases. For high order methods this task becomes computationally
Table 1. Stability intervals of H2 for Stormer's and Cowell methods and the interval of periodicity of H2 for symmetric methods
Order
7
Störmer Cowell
360; 60N
323' 49 60
°'n
unstable
60
Ml
10
/ 27
(0; 0.3820447...) 0'-
V ' ' V 128
( 87280 189 A /4221504 189<
\308407 52
1824647 71
Order
9
11
12
unstable unstable
(0; 0.3597184 ...) (0; 1.0218233 ...) ^0.
V 595163) unstable
Stoörmer Cowell
51975
1686934
20790
0.1898631'-
28687,
Order
6
8
10
12
Symmetric
(0; 0.8021734 ...) (0; 0.5157665 ...) (0; 0.1724269 ...) (0; 0.0456343 ...)
6
intensive. For example, when the Jury criterion is applied to the 8th order Stormer method the maximum degree of the polynomial equals to 12, and for the 8th order Cowell method it equals to 117!
We propose the following effective method to determine the boundaries of the stability interval of method (1). We have to define all H2 for which the polynomial n(r; H2) has a root that belongs to the unit circle. Consider the roots r* = exp(ip) and r* = exp(—ip) of (6), 0 < p < n. Let us represent n(r; H2) in the form:
n(r; H2) = S(r; H2)(r2 — 2r cos p + 1) + R(r; H2)
where S(r; H2) is a polynomial of the order (k — 2) in r with real coefficients, R(r; H2) = = a0(H2, cos p) + a1(H2, cos p)r, a0, a1 G R. Since r* and r* are the roots of both polynomials n(r; H2) and r2 — 2r cos p + 1, R(r; H2) = 0. Therefore a0(H2, cos p) = 0 and a1(H2, cos p) = 0. Consider solutions (H*,p*) of the last two equations, where —1 < cosp* < 1. In addition, the case p = 0 gives H* = 0 and the case p = n immediately gives H* = —a( — 1)/p( — 1). Choose all H* G R+ only, and they divide R+ into disjoint intervals. We test polynomial (6) using the Jury criterion for strong stability for some value of H2 belonging to each interval. The interval in R+ for which n(r; H2) is strongly stable corresponds to the interval of absolute stability of method (1).
Tab. 1 presents the absolute stability intervals for the Stormer and Cowell methods of orders from 5 to 12 . The results show that not all methods are stable at small H2. For example, the Cowell method of order 8 has a very short stability interval separated from zero. The presented results are the same as those from [13], with the exception of the 7th order Stormer method for which one more root was found. It is close to but it does not agree with that found in [13]. This reduces the stability interval. In addition, rational boundaries of the stability intervals can be found with our approach find if they exist.
Interval of periodicity of symmetric methods. If £ is a root of a symmetric polynomial then 1/£ is also its root. Then for symmetric method (1) there is no such H2 that all roots of the stability polynomial n(r; H2) are in the unit circle. Therefore, any symmetric method is absolutely unstable. On the other hand, symmetric methods can have another useful property, namely, they can have a non-vanishing interval of periodicity [14].
According to [14] method (1) has non-vanishing interval of periodicity (0; H2) if for all H2 £ (0; H2) the roots rs of the stability polynomial n(r; H2) satisfy relations
r1 = exp (iO(H)), r2 = exp (-iO(H)), |rs| = 1, s = 3,...,k, 6(H) £ R
and if the order of (1) is p then 6(H) = H + O(hp+1) £ R.
If method (1) has non-vanishing interval of periodicity then it is symmetric. The opposite is not true, but if polynomial p of symmetric method (1) has all roots in the unit circle and there are no other double roots except the principal ones then the method has a non-vanishing interval of periodicity. In this case, since the roots of the polynomial continuously depend on parameter H2, all roots of n(r; H2) remain in the unit circle when H2 changes from 0 to some H2. Then the principal component (xn)P of the numerical solution is periodic with a period close to 2n/A (the period of the analytical solution of (3)), and (xn)P dominates over (xn)S which is also periodic.
The approach to determine the value of H0 is proposed [14]. Some polynomial is constructed from n(r; H2) by special transformation of variable r [14]. The value of H2 is determined from the condition that all roots of the polynomial are real, distinct and non-negative. This corresponds to the condition that the absolute values of all roots of n(r; H2) are equal to 1 for H2 £ (0; H2).
We propose an alternative method based on determining of H in such a way that multiple root arises for n(r; H2). Let symmetric method (1) has a non-vanishing interval of periodicity (0; H2)> Then for H £ (0; H2) all roots of n(r; H2) are distinct and lie on the unit circle. Moreover, each root that does not lie on the real axis has the conjugate root as the root of a polynomial with real coefficients (Fig. 2 a). If H2 > H2 then there exists £* = r* (cos 6* + i sin 6*) root of n(r; H2), where r* > 1. Therefore 1/£* = (cos 6* — i sin 6*)/r* and its conjugate = = r* (cos 6* — i sin 6*), 1/£* = (cos 6* + i sin 6*)/r* are also the roots of n(r; H2). Because r* is continuously depends on parameter H there exists H = H* for which r* = 1, that is, £* is root of multiplicity 2. Therefore, H* coincides with the right-hand boundary of the interval of periodicity H . Thus, H can be found as the minimum positive real root of the discriminant of the stability polynomial of a symmetric method. In Fig. 2, the behaviour of the roots of the stability polynomial for the 8th order symmetric method is shown when H approaches H0 and when H is greater than H0. Table 1 shows the interval of periodicity for the symmetric methods considered here.
The Stormer methods have a non-vanishing absolute stability interval but do not have an interval of periodicity. Alternatively, symmetric methods are absolutely unstable but they have a non-vanishing interval of periodicity. These differences are shown in Fig. 3 for the following simple numerical example. Let us consider problem (3) with the initial conditions x(0) = 1 and x'(0) = 0. Then the exact solution is x(t) = cos(At). Equation (3) has two the first integrals
x'(t)
E := A2(x(t))2 + (x!(t))2 = const, 6 := At + arctan —= const.
Ax(t)
Although the velocity x'(t) = v(t) is not directly defined by (3), it can be determined by equation (3) through introduction of unknown function v with the initial data v(0) = 0,
v' (0) = x' '(0) = —A2.
The error of the first integral AE = Eh — E for (3) is shown in Fig. 3 (a) and (b) for A =1. equation (3) is integrated with the 8th order symmetric method and Stormer method, respectively. The error of the first integral A6 = 6h — 6 is demonstrated in Fig. 3 (c) for both methods. Here E, 6 are exact values of the first integrals (they equal to 1 and 0, respectively) and Eh, 6h are numerical values of the first integrals. The step-size h = n/128 belongs to the
f
t—1
■ 1
V
0 0.05 0.L 0.15 0.2 0.25
a) 8th order Stormer
b) zoom in 8th order Stormer c) 10th order Stormer (unstable)
-0.5 l 2.5
d) 12th order Stormer
12 3 4 e) 8th order Cowell
f) zoom in 8th order Cowell
0 2 4
g) 9th order Cowell
h) zoom in 9th order Cowell
■3 2 7 12 17 22
i) 10th order Cowell
Fig. 1. The root-locus curves and the stability regions for some Stormer and Cowell methods in the complex plane represented by H2
interval of periodicity of symmetric method and to absolute stability interval of the Stormer method.
The symmetric method gives a periodic solution, therefore Eh is a periodic function with constant amplitude. One can see in Fig. 3 (a) that the energy of the system does not increase with time. Since all roots of n(r; H2) for the Stormer method are less than 1, the energy of the numerical solution decreases (Fig. 3 (b)). Since the period of the numerical solution does not
a)
b)
c)
d)
Fig. 2. The roots of stability polynomial n(r; H2) for the 8th order symmetric method (1) in the complex plane; r is shown for H2 = H2/10 (a), H2 = 9H2/10 (b), H2 = H2 (c), H2 = 11H2/10
(d). The roots r1 and r2 that correspond to the perturbed principal roots with black triangle marker
£2 = 1 are marked
a)
b)
c)
Fig. 3. The errors of the first integrals for the equation of harmonic oscillator for the symmetric (a,c) and Stormer (b,c) methods with step-size h/T = 256, t in radians.
coincide with the theoretical one, the numerical solution is either ahead of or lagging behind the exact solution. For the symmetric method |A6| grows slower than for the Stormer method (Fig. 3 (c)).
2. Numerical experiments
Let us consider two three-dimensional model problems that have exact solutions. By "exact solution" is meant a solution that can be obtained by integrating Kepler's equation. Model problem 1 is the three-dimensional Kepler problem
x''(t) =1X3,
(7)
where n is the standard gravitational parameter, x = (x1,x2,x3) is the radius-vector of the satellite and | x| is the Euclidean norm of x.
Model problem 2 is specially constructed from the restricted three-body problem (Earth-Moon-Earth's satellite of negligible mass). In this problem, the force acting from the Moon on the satellite is compensated by an additional force which depends only on time and it is independent of the position of the satellite on the orbit. This force affects the movement of the
satellite in such a way that the exact solution of the problem describes the movement of the satellite around the Earth in the absence of the Moon.
Let us consider the equation of motion for a satellite of negligible mass in an inertial reference frame centred at the Earth-Moon barycentre
HU\ X - XE X - XM I cn\ to\
X (t) = -¡3 - HM]-¡3 + f (t). (8)
|x - XE |3 |x - XM I3
Here X = (xi,x2,хз), xe = ((xe)i, (xe)2, (xe)з) and xm = ((xm)i, (xm)2, (хм)з) are positions of the satellite, the Earth and the Moon, respectively; цЕ=3.986004419E+14 m3/s2 and fj,M =4.9048696E+12 m3/s2 are the standard gravitational parameters of the Earth and the Moon; f(t) = (fi(t), f2(t), f3(t)) is an additional force. The coordinates of the Earth and the Moon are determined by two-body problem. Let xES be the exact solution of the Kepler problem (7) for the system Earth-satellite. Then x = xES + xE and
f f ,n XES + XE — XM XE — XM
I(t) = hm~,-:-[3 - MM]-л • (9)
IXES + XE - XM I3 IXE - XM |3
Thus, model problem (8)-(9) has the exact solution, and the errors of the numerical solution are calculated directly. Since the Jacobian of the problem coincides with the Jacobian of the restricted three-body problem the stability properties of the numerical methods for these problems coincide.
The following initial orbital parameters are adopted in numerical experiments. For the Moon they are aM = 3.94748E + 08 m, eM = 0.042200, шм = 22°8", Пм = 4o40 ", iM = 18°31", (Mo)m=340° 13 ". For satellite they are asat= 2.5500000004E+07 m, £sat= 0.00068, ^Sat= 135.0000214°, ^Sat=120°, isat= 64.9°, (Mo)sat= 32.6650111°, Tsat= 11h15'44 ".
We compare the accuracy of the orbit integration by the Stormer method and symmetric methods. Additionally, the results for the Bashforth method are shown in the case when problems (7) and (8) are represented in the form of six first-order ODEs. To improve the accuracy of the Adams methods the predictor-corrector scheme is also used in the form P(EC)3E, where the right-hand side (E) and the corrector (C) are evaluated three times at each step. Since the absolute stability intervals of the 8th order Stormer and Cowell methods do not coincide, calculations were carried out for the case when the orders of the predictor and corrector coincide, and they are equal to 8 (P8 (EC8)3E), and for the case of the 9th order corrector (Pg(ECg)3E).
For each Model problem, we are interested in the maximum deviation of the calculated satellite position from the exact one after integration for about a year. Let us denote the numerical and exact solutions at the moment tn by xn and xex(tn), respectively, n = k,... ,K, tK = 779Tsat, Tsat is a period of satellite. The following notations for errors are used
К = xhn - Xex(tn), Ah = maXK | (Д^ |, ph = max^ |ДЩ.
In addition, we consider the decomposition of the error vector ДП in terms of the basis vectors associated with the exact ellipse. They are ro(tn) = Xex(tn )/IXex(tn )|, To (tn) = = Vex(tn)/|Vex(tn)| and no = r0(tn) x T0(tn). Then we have
Sh = mac |ro(tn) • ДЩ, Sh = rnnxK |To(tn) • ДЩ, % = ™axK |no(tn) • ДЩ
n=k,...,K n=k,...,K n=k,...,K
Results of calculations with fixed step-size h = Tsat/512 for the Model problems 1 and 2 In Tabs are presented in 2, 3, respectively. One can see that the direct solving of the second
order ODE is more efficient. Explicit symmetric methods give more accurate results even in comparison with the explicit-implicit PECE algorithms. In Model problem 1 the Stormer-Cowell PECE algorithm offers slight advantage over other algorithms only in the error along the radius. However, symmetric method offers advantage over other algorithms in calculating positions of the satellite. The symmetric algorithm symmetric method offers advantage over other algorithms in the case of Model problem 2.
Table 2. Accuracy of the orbit integration (h = Tsat/512). Model problem 1
Bashforth- Stoormer- Stormer-
Bashforth Moulton Ps(ECg)3E Stoormer Cowell Ps(ECg)3E Cowell PS(ECQ)3E symmetric
Ah, m 7.43E-03 1.77E-04 7.06E-04 2.04E-05 8.45E-06 1.61E-06
Ah, A2 , m 1.07E-02 2.55E-04 1.01E-03 2.93E-05 1.21E-05 2.20E-06
Ah, m 1.08E-02 2.59E-04 1.03E-03 2.98E-05 1.23E-05 2.27E-06
ph, m 1.20E-02 2.86E-04 1.14E-03 3.29E-05 1.36E-05 2.60E-06
sh, m 4.98E-06 1.06E-07 3.82E-07 1.73E-08 1.30E-08 1.42E-07
sh, m 1.20E-02 2.86E-04 1.14E-03 3.29E-05 1.36E-05 2.60E-06
sh, m 1.77E-24 1.74E-24 1.40E-22 7.03E-23 1.17E-22 5.07E-23
Table 3. Accuracy of the orbit integration (h = Tsat/512). Model problem 2
Bashforth- Stoormer- Stoormer-
Bashforth Moulton Ps(ECg)3E Stoormer Cowell Ps(ECg)3E Cowell Ps(ECQ)3E symmetric
Ah, m 1.36E-01 3.23E-03 1.27E-02 3.62E-04 1.56E-04 8.89E-05
Ah, A2 , m 1.95E-01 4.64E-03 1.82E-02 5.20E-04 2.23E-04 1.28E-04
Ah, m 1.97E-01 4.70E-03 1.84E-02 5.27E-04 2.26E-04 1.29E-04
ph, m 2.19E-01 5.22E-03 2.05E-02 5.85E-04 2.51E-04 1.43E-04
sh, m 3.44E-04 8.14E-06 3.20E-05 9.05E-07 4.02E-07 3.70E-07
sh, m 2.19E-01 5.22E-03 2.05E-02 5.85E-04 2.51E-04 1.43E-04
sh, m 7.55E-04 1.80E-05 9.56E-06 2.73E-07 1.17E-07 6.77E-08
Table 4. Accuracy of the orbit integration (h = Tsat/den). Model problem 1
Bashforth- Stormer- Stormer-
Bashforth Moulton Stormer Cowell Cowell symmetric symmetric
Pb(ECb)3 E Pb(ECb )3E Pb(ECq)3E
rhp 486876 1286888 377037 1012680 922316 174497 253176
den 625 413 484 325 296 224 325
h, sec 64.8 98.1 83.7 125 137 181 125
H2 1.01E-04 2.31E-04 1.69E-04 3.74E-04 4.51E-04 7.87E-04 3.74E-04
ph , m 1.99E-03 1.98E-03 1.89E-03 1.98E-03 1.89E-03 1.91E-03 9.80E-05
slh, m 6.56E-07 7.51E-07 6.86E-07 7.64E-07 1.81E-06 1.07E-04 5.41E-06
sh, m 1.99E-03 1.98E-03 1.89E-03 1.98E-03 1.89E-03 1.91E-03 9.80E-05
m 1.22E-24 1.12E-24 1.02E-22 4.49E-23 5.65E-23 2.41E-23 5.57E-23
Another series of calculations were carried out to determine the step at which the maximum deviation of the numerical solution from the exact one does not exceed 2 mm for a year. The
Table 5. Accuracy of the orbit integration (h = Tsat/den). Model problem 2
Bashforth- Stormer- Stöormer-
Bashforth Moulton Stöormer Cowell Cowell symmetric symmetric
Ps (ECs)3E Ps(ECs)3E Ps(ECg)3E
rhp 671499 1779216 519594 1402180 1274424 289789 350551
den 862 571 667 450 409 372 450
h, sec 47.0 71.0 60.8 90.1 99.1 109 90.1
H2 5.31E-05 1.21E-04 8.87E-05 1.95E-04 2.36E-04 2.85E-04 1.95E-04
ph, m 2.00E-03 1.95E-03 1.87E-03 1.88E-03 1.89E-03 1.84E-03 4.02E-04
S1}, m 3.10E-06 3.03E-06 2.91E-06 2.92E-06 3.03E-06 4.75E-06 1.04E-06
Sh, m 2.00E-03 1.95E-03 1.87E-03 1.88E-03 1.89E-03 1.84E-03 4.02E-04
5vh, m 6.89E-06 6.72E-06 8.75E-07 8.78E-07 8.86E-07 8.69E-07 1.90E-07
results of calculations are presented in Tabs. 4, 5. The first row marked "rhp" shows the number of evaluations of the right-hand side that were required to achieve the accuracy. In the last column, the results are presented for the symmetric method with the step it takes the Stormer-Cowell PECE algorithm to achieve the specified accuracy. The advantage of the symmetric method is obvious, especially for Model problem 2. In addition, the symmetric methods have the lowest number of right-hand side evaluations in comparison with other methods considered.
This work was supported by the Krasnoyarsk Mathematical Center and financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Censers for Mathematics Research and Education (Agreement no. 075-02-2020-1631).
References
[1] IGS ftp archives, ftp://ftp.igs.org/pub/center/analysis/. Last accessed 4 Aug 2020
[2] E.Everhart, Implicit Single-Sequence Methods for Integrating Orbits, Celestial Mechanics, 10(1974), 35-55.
[3] G.Beutler, Numerische Integration gewöhnlicher Differentialgleichungssysteme: Prinzipien und Algorithmen. Mitt. Satell, Beobachtungsstn. Zimmerwald, 23(1990).
[4] G.Beutler, Methods of Celestial Mechanics I: Physical, Mathematical, and Numerical Principles, Springer-Verlag, Berlin, 2005.
[5] G.Quinlan, S.Tremaine, Symmetric multistep methods for the numerical integration of planetary orbits, Astron. J., 100(1990), no. 5, 1694-1700.
[6] P.Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley and Sons, New York, 1969.
[7] J.D.Lambert, Computational Methods in Ordinary Differential Equations, John Wiley and Sons, New York, 1973.
[8] T.Bordovitcina, The modern numerical methods in problems of celestial mechanics, Nauka, Moscow, 1984 (in Russian).
[9] E.Yairer, S.Norsett, G.Wanner, Solving Ordinary Differential Equations, Springer-Verlag, Berlin, 1987.
[10] E.Vergbitckii, Basis of Numerical Methods, Vysshaya shkola, Moscow, 2004 (in Russian).
[11] V.Avdushev, Numerical modeling of orbits, Izdat. NTI, Tomsk, 2010 (in Russian).
[12] J.C.Butcher, Numerical methods for ordinary differential equations, John Wiley and Sons, New York, 2016.
[13] S.N0rsett, A.Asheim Regarding the absolute stability of Stormer-Cowell methods, Discrete and Continuous Dynamical Systems, 34(2014), no. 3, 1131-1146.
DOI: 10.3934/dcds.2014.34.1131
[14] J.D.Lambert, Symmetric Multistep Methods for Periodic Initial Value Problems, J. Inst. Maths Applies, 18(1976), 189-202.
[15] P.Chakravarti, P.Worland, A class of self-starting methods for the numerical solution of y'' = f (x,y), BIT Numerical Mathematics, 11(1971), no 4, 368-383.
[16] A.Hurwitz, On the conditions under which an equation has only roots with negative real parts (English translation by H. G. Bergmann), in Selected Papers on Mathematical Trends in Control Theory, R. Bellman and R. Kalaba Eds., Dover, New York, 1964, 70-82.
[17] E.Jury, J.Blanchard, A stability test for linear discrete systems in table form, I.R.E. Proc., 49(1961), 1947-1948.
[18] E.Jury A modified stability table for linear discrete systems, Proc. IEEE, 53(1965), 184-185.
[19] E.Jury, Inners and the Stability of Linear Systems, John Wiley and Sons, New York, 1982.
Точность симметричных многошаговых методов численного моделирования движения спутника
Евгения Д. Карепова
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Илья Р. Адаев
Сибирский федеральный университет Красноярск, Российская Федерация Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Юрий В. Шанько
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Аннотация. В статье мы подробно обсуждаем устойчивость линейных многошаговых симметричных методов высокого порядка в задаче гармонического осциллятора. Приведены эффективные алгоритмы вычисления интервалов абсолютной устойчивости и периодичности. Численные эксперименты демонстрируют точность вычисления орбиты на интервале около одного года для трехмерной задачи Кеплера и для специально разработанной трехмерной тестовой задачи, которая моделирует систему Земля-Луна-спутник и имеет точное решение.
Ключевые слова: линейные многошаговые методы, симметричный метод, методы Адамса-Штермера-Коуэлла, PECE-схема, орбита.
