- Email: [email protected]
Increasing the accuracy in the numerical integration of perturbed oscillators with new methods
Applied Numerical Mathematics 41 (2002) 295–304 www.elsevier.com/locate/apnum
Increasing the accuracy in the numerical integration of perturbed oscillators with new methods Ana B. González ∗ , Pablo Martín, Amelia García, José M. Farto Departamento de Matemática Aplicada a la Ingeniería, E.T.S. de Ingenieros Industriales, Universidad de Valladolid, Paseo del Cauce s/n, 47011 Valladolid, Spain Received 10 October 2001
Abstract An important problem in celestial mechanics is the precise long-term prediction of satellite orbits. This is a difficult problem due to different factors. When the problem is formulated in BF or KS variables its structure is of the form of perturbed oscillators. The effect of the choice of the frequency in the integration of perturbed oscillators leads us to construct new methods which exhibit a very good behaviour. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction The precise long-term prediction of Earth satellite orbits is an important problem in celestial mechanics. Numerical solutions can move or spiral away from the exact solution due to Lyapunov or orbital instability, and the integration remains valid only for a span of time that could be too short. It is very difficult to increase the accuracy, since the problem is affected by different factors, ranging from deficiencies in the force model, to uncertainties in the initial values, to the simple accumulation of errors during the integration. We are concerned here only with integration errors, and in order to minimize these, we transform the equations of motion to others suitable for the numerical integration by means of changes of variables such as the linearizing transformations as KS (Stiefel and Scheifele [13]) or BF variables (Ferrándiz [5]). In these variables the equations of motion are reduced to a system of perturbed oscillators. In 1971 Scheifele [12] obtained a refinement of the Taylor method based on his G-functions for the integration of perturbed oscillators. But, the Scheifele G-function method has the disadvantage that it can only be applied in very particular cases, due to the complexity of the preliminary calculations that it requires. Recently, new methods based on the Scheifele G-functions (SMF codes of Martín and * Corresponding author.
E-mail address: [email protected] (A.B. González). 0168-9274/01/$22.00 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 1 2 1 - 0
296
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
Ferrándiz [7] and RKGM schemes of González et al. [6]) have been developed to avoid such problems. But most of the special methods require a previous knowledge of a good approximation to the frequency. The effect of this frequency in the integration of perturbed oscillators is very important as we can see in Kirchgraber [9], Novo and Rojo [10] and Martín and Farto [8]. This fact leads us to formulate new purely Scheifele-based methods which inherit its good behaviour in order to improve by selecting a suitable frequency. In particular, we construct a fourth order scheme which shows very good results in the numerical integration of perturbed oscillators. The choice of a suitable approximation to the frequency for this method, leads to an important reduction in error. Finally, the new algorithm is applied to the propagation of artificial Earth satellite orbits once the corresponding equations of motion have been reduced to oscillatory form by using linearizing variables.
2. The solution of a perturbed oscillator Let us consider the model problem ˙ y¨ + ω2 y = ε f (t, y, y),
y(0) = y0 ,
y(0) ˙ = y˙0 ,
(1)
where ε is a small parameter. If we suppose that the function f (t, y(t), y(t)) ˙ admits an expansion of the form
˙ = f t, y(t), y(t)
∞
(j ) t
f0
j =0
∞
,
the j th derivative of f at t = 0), then the initial value problem (1) can be
(j ) t
j
f0
j =0
j!
(j )
(we have denoted by f0 written as y¨ + ω2 y = ε
j
j!
,
y(0) = y0 ,
y(0) ˙ = y˙0 .
(2)
The solution of (2) can be obtained in the usual way, building the solution of the homogeneous equation with the given initial conditions and adding to it the solution of the non-homogeneous equation in which the solution and its derivative at t = 0 are cancelled out. The second solution can be calculated by solving the initial value problems y¨j + ω2 yj =
tj , j!
yj (0) = y˙j (0) = 0,
j 0,
(3) (j )
and linearly combining their solutions with the coefficients εf0 . Scheifele defines his G-functions as the solutions to the problems (3) with the notation Gj (t) = yj −2 (t) for j 2. The solution to the homogeneous problem will be obtained by simply combining the solutions of the problems y¨ + ω2 y = 0, y¨ + ω2 y = 0,
y(0) = 1, y(0) = 0,
y(0) ˙ = 0, y(0) ˙ = 1,
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
297
which Scheifele denotes as G0 (t) and G1 (t), respectively, and suggests a systematic way to obtain the solution of the problem (1) which is given by: Y (t) = y0 G0 (t) + y˙0 G1 (t) + ε
∞
(j )
f0 Gj +2 (t).
(4)
j =0
The Scheifele G-functions have the following expressions in terms of the elementary functions
j −1
k t 2k (−1)j −ω2 , G2j (t) = 2j cos ωt − ω (2k)! k=0
(−1)j G2j +1 (t) = ω2j
j −1
k t 2k+1 1 sin ωt − −ω2 , ω (2k + 1)! k=0
j 0,
(5)
or in terms of powers of t as follows Gj (t) =
∞
−ω2
k
k=0
t 2k+j , (2k + j )!
j 0.
(6)
See Scheifele [12]. Let h a real number. If we truncate expression (4) up to j = p − 2, we obtain an approximation of order p to the solution at the point t = h. In spite of its good behaviour, the Scheifele G-function method has the great disadvantage that it can only be applied in very particular cases, due to the complexity of the preliminary calculations that it requires. This problem was solved by Martín and Ferrándiz [7] in a drastic way, through the conversion of the G-function method into a multistep scheme, or more recently, by González et al. [6] who construct a new one-step approximation to the solution. But most of the special methods need a good approximation to the frequency of the problem. We shall analyze in the following section the effect of such frequency on the numerical solution provided by some special methods.
3. The effect of the frequency on the solution of a perturbed oscillator It is well known that the choice of a good approximation to the actual frequency of the solution, improves notably the accuracy of some kinds of numerical methods when integrating weakly perturbed oscillators (see Kirchgraber [9], Novo and Rojo [10] and Martín and Farto [8]). We shall discuss this phenomenon. To this end, we consider the simple and widely used test problem y¨ + y = εy 3 ,
y(0) = 1,
y(0) ˙ = 0.
(7)
The true frequency of the solution can be easily obtained by means of perturbation techniques (Lindstedt method, Bush [1]) as a power series of ε. A first order approximation to this frequency is 3 ω¯ = 1 − ε. 8 In this way, we can rewrite problem (7) as
y¨ + ω¯ 2 y = εy 3 − 1 − ω¯ 2 y,
y(0) = 1,
y(0) ˙ = 0.
(8)
298
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
Fig. 1. y¨ + y = εy 3 , y(0) = 1, y(0) ˙ = 0, ε = 10−3 , h = 0.1.
While for the integration with a Taylor method, the use of one or other form is irrelevant, when integrating with a Scheifele method, the selection of the expression (7) or (8) is very important. This last method depends explicitly on the frequency of the unperturbed oscillator considered. Fig. 1 shows the integration of the test problem with a classical embedded pair RKN4(3) (Dormand et al. [2]), and the schemes depending on the G-functions (and therefore on the frequency), the multistep code SMF of Martín and Ferrándiz [7] and the one-step scheme RKGM of González et al. [6]. We can observe that the results obtained with the SMF method when considering as frequency ω = 1 − 38 ε are better than those corresponding to ω = 1. The right choice leads to a remarkable stabilization of the error. These results have been improved thanks to the elimination of the resonance that appeared. But surprisingly, the RKGM methods do not seem to inherit the good behaviour to the change of the frequency of the Scheifele scheme. The RKGM scheme do not show very much reduction in error with the correction to the frequency. In fact, the improvement achieved is irrelevant. A further analysis is necessary to understand this behaviour. This will be done in the next section.
4. A new approximation for the perturbed oscillator The RKGM method
y1 = y0 G0 (h) + y˙0 G1 (h) + ε
p−2
s
j =0
i=1
y˙1 = y˙0 G0 (h) − ω y0 G1 (h) + ε 2
ki = f ci h, y˜i ,
i = 1, . . . , s,
bj i ki Gj +2 (h),
p−2
s
j =0
i=1
bj i ki Gj +1 (h),
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
y˜1 = y0 ,
y˜i = y0 + ci hy˙0 + h2
i−1
299
aij εkj − ω2 y˜j ,
j =1
appears as a kind of mixed method. While the final approximation to the solution is a Scheifele-like formula, the internal stages ki are obtained using arguments based on Taylor series to approximate y(ci h). But y(t) is the solution of a perturbed oscillator, and thus it seems more natural to approximate the exact solution at the points ci h using the solution (4) provided by Scheifele. This characteristic may be the cause of the unexpected response of the RKGM method to the change of frequency. We should not want to relinquish to the excellent properties of the RKGM method when integrating satellite orbits. If they could be arranged so as to avoid the deficiency noted, then the performance could be improved. Unfortunately, the problem of RKGM methods is structural. For this reason, the development of a new approximation preserving the good properties of the RKGM method, and purely Scheifele-based, is worth the effort. In fact, this can be done as in the following example. We begin by considering the Scheifele method of order 4, Y1 = y0 G0 (h) + y˙0 G1 (h) + ε
2
(j )
f0 Gj +2 .
j =0
We shall try to approximate the first two derivatives of f as follows (j )
f0 ≈ βj 1 g1 + βj 2 g2 + βj 3 g3 ,
j = 0, 1, 2,
with βj i real numbers and gi evaluations of f at the points t¯ = 0, t¯ = h/2 and t¯ = h. Obviously, y(0) = y0 , then we take g1 = f (0, y(0)) = f (0, y0 ). Now we need to approximate y(h/2) and y(h). From (4) we have y(h/2) = y0 G0 (h/2) + y˙0 G1 (h/2) + ε
∞
(j )
f0 Gj +2 (h/2),
(9)
j =0
then y(h/2) ≈ y¯2 = y0 G0 (h/2) + y˙0 G1 (h/2) + εg1 G2 (h/2), thus we take g2 = f (t¯2 , y¯2 ). A similar idea for y(h) allows us to consider y¯3 = y0 G0 (h) + y˙0 G1 (h) + εg2 G2 (h), and g3 = f (t¯3 , y¯3 ). In order to obtain the parameters βj i , it is obvious that β01 = 1 and β02 = β03 = 0. For the values ˙ y )|t =0 and the following β1i , i = 1, 2, 3, we estimate the agreement between terms of f0(1) = (ft + yf Taylor series about h = 0 β11 g1 + β12 g2 + β13 g3 = (β11 + β12 + β13 )f + ft (β12 /2 + β13 ) + fy y˙0 (β12 /2 + β13 ) h
+ ft t (β12 /8 + β13 /2) + fty y˙0 (β12 /4 + β13 ) + fy εf − ω2 y (β12 /8 + β13 /2)
+ fyy y˙ 2 (β12 /8 + β13 /2) h2 + O h3
300
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
obtaining the system: β11 + β12 + β13 = 0, β + 2β13 = 2/ h, (10) 12 β12 + 4β13 = 0, whose solution is β11 = −3/ h, β12 = 4/ h, β13 = −1/ h. For the values β2i , i = 1, 2, 3, we estimate the agreement between terms of f0(2) = (ft t + 2fty y˙ + y˙ 2 fyy + fy (εf − ω2 y))|t =0 and the Taylor series about h = 0 β21 g1 + β22 g2 + β23 g3 = (β21 + β22 + β23 )f + ft (β22 /2 + β23 ) + fy y˙0 (β22 /2 + β23 ) h
+ ft t (β22 /8 + β23 /2) + fty y˙0 (β22 /4 + β23 ) + fy εf − ω2 y (β22 /8 + β23 /2)
+ fyy y˙ 2 (β22 /8 + β23 /2) h2 + O h3
obtaining the system: β21 + β22 + β23 = 0, β + 2β23 = 0, 22 β22 + 4β23 = 8/ h2 ,
(11)
whose solution is β21 = 4/ h2 , β22 = −8/ h2 , β23 = 4/ h2 . The final approximation to the solution takes the form y1 = y0 G0 (h) + y˙0 G1 (h) + ε
2 3
βj i gi Gj +2 (h),
j =0 i=1
y˙1 = y˙0 G0 (h) − ω2 y0 G1 (h) + ε
2 3
βj i gi Gj +1 (h).
j =0 i=1
Fig. 2. y¨ + y = εy 3 , y(0) = 1, y(0) ˙ = 0, h = 0.1, ε = 10−3 . Frequency for RKGM and SMF ω = 1 − (3/8)ε.
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
301
From the construction and property (6) we have (j )
(βj 1 g1 + βj 2 g2 + βj 3 g3 )Gj +2 = f0 Gj +2 + O h5 , (βj 1 g1 + βj 2 g2 + βj 3 g3 )Gj +1 =
(j ) f0 Gj +1
+O h , 4
j = 1, 2, j = 1, 2.
Therefore, at least we achieve order four in the approximation to the solution and order three in the derivative. But a simple Taylor series of the approximation to the derivative with these particular coefficients shows that we also achieve order four in this case. As a result of removing all the Taylor inheritance from the new method, we can see in Fig. 2 its good response to the frequency modification technique. Thus, we have reached our target, and also we have constructed a new method with better behaviour than the RKGM scheme. In a variety of satellite problems, the true frequency is very difficult to approximate, or such approximation simply can not be obtained. In these cases we have obtained an improved numerical method. This fact justifies the construction of the new formula. 5. Reducing the error growth in long-term computation of satellite orbits This section is devoted to a study of the behaviour of the new approximation constructed in the previous section when integrating the satellite problem. To this end, we consider the problem expressed in linearizing, using KS (Stiefel and Scheifele [13]) or BF (Ferrándiz [5] and Ferrándiz et al. [11]) variables. In KS variables the equations of motion become 1 |Y |2 , ω˙ = R, (12) t˙ = y¨i + yi = Ri , i = 1, 2, 3, 4, 4 ω √ where Y = (y1 , y2 , y3 , y4 ) is the vector of the KS coordinates, ω = 2p0 , p0 stands for the opposite of the total energy and Ri and R denote the perturbation terms. The dot denotes the derivative with respect to an eccentric anomaly-like variable s. The problem in focal variables (BF variables) is reduced to four perturbed harmonic oscillators with unit frequency. The coordinates of the basic set of focal variables are three components (y1 , y2 , y3 ) of the direction vector of the particle and the inverse u of the radial distance. The system is completed with the addition of an equation which provides the variation of the angular momentum c (which has a stabilizing effect on the solutions) and with the equation corresponding to the change of time. Then the equations of motion can be described as y˙i y¨i + yi = Qi − c˙ , i = 1, 2, 3, c u µ u¨ + u = 2 + Q − c˙ , (13) c c 1 F t˙ = 2 , c˙ = , c cu where µ is the reduced mass, Qi , Q and F stand for the corresponding disturbing terms, and the dot denotes the derivative with respect to a generalized true anomaly variable τ . We have considered some orbits with different initial values and different sets of perturbations and we have investigated the error produced during the numerical propagation of the orbits over several tens of revolutions.
302
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
We consider an almost periodic equatorial orbit with the coefficient J2 due to the oblateness as the perturbation parameter. We have neglected terms of the order of J22 . The eccentricities considered are e = 0 and e = 0.99. Thus, for BF variables the problem can be written as µ J2 2 + 12 u, (14) c2 c2 where µ/c2 = 20/21 and J2 /c2 = (10/21) × 10−3 for e = 0 and µ/c2 = 100/20895 and J2 /c2 = (50/20895) × 10−3 for e = 0.99. We present efficiency curves for each problem, where the common logarithm of the maximum global error over all steps and variables is plotted against the number of evaluations of the function required. These graphics give a good indication of the efficiency of the considered methods when they are carefully implemented. To this end, we rewrite the final approximation of the RKGM method as y¨i + yi = 0,
i = 1, 2, 3,
y1 = y0 G0 (h) + y˙0 G1 (h) + ε
u¨ + u =
p−2 s i=1
y˙1 = y˙0 G0 (h) − ω y0 G1 (h) + ε 2
bj i Gj +2 (h) ki ,
j =0
p−2 s i=1
bj i Gj +1 (h) ki ,
j =0
and we compute the coefficients of the internal stages ki at the beginning of the integration, taking into account the expression (6) for the G-functions. We use the same idea for the new codes. We compute at the beginning of the integration the coefficients of gi and the G-functions evaluated at ci h. Thus, the overhead per function evaluation is approximately the same for these two methods and a fixed step size RKN scheme. On the other hand, when considering problems with much more costly right-hand side functions, the overhead costs become much less significant. Fig. 3 shows the results obtained when integrating 1000 revolutions of the problem in BF variables and eccentricities e = 0 and e = 0.99. In this case we have a reference solution (see Farto et al. [4])
Fig. 3. BF variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 1000 revolutions. Left: e = 0. Right: e = 0.99.
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
303
and the error is computed at each point of the numerical integration. The suitable frequency for e = 0 is ω = 0.9945578231292517 and for eccentricity e = 0.99 is ω = 0.9986257476147703. We have denoted by New1 and New2 the results of the new method with frequencies ω = 1 and the corrected frequency, respectively, and similarly for SMF1 and SMF2. We have used as step sizes h = 1/2i , i = 1, . . . , 7, for the RKGM and new methods and h = 1/2i , i = 1, . . . , 8, for the SMF scheme. The tolerances used are 10−3 , . . . , 10−9 for the embedded pair. In both cases, low and high eccentricity, the best results correspond to the new approximation. Figs. 4 and 5 correspond to the case of KS variables. In this case, an expression for the solution is not available. Thus we have obtained at perigees and apogees a reference numerical solution, in which we have also computed the numerical methods. For this reference solution, we have considered the RKN
Fig. 4. KS variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 20 revolutions. Left e = 0. Right: e = 0.5.
Fig. 5. KS variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 20 revolutions. e = 0.99.
304
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
method of order eight constructed by Dormand et al. [3] and implemented with fixed step size. We have used arithmetic precision of 32 significant figures and step size h = π/29 . For e = 0 and e = 0.5, there are six runs corresponding to h = π/23, . . . , π/28 for the RKGM method and h = π/22 , . . . , π/27 for the new scheme, and tolerances 10−3 , . . . , 10−6 for the embedded pair. For e = 0.99, there are runs with h = π/23 , . . . , π/27 for the new scheme and the same tolerances as before for the RKN. For the first two eccentricities the best results correspond to the new method, while for e = 0.99, new and RKGM schemes present a similar behaviour. This fact is not surprising, since the equations in KS variables are singular at e = 1 (see Ferrándiz [5]). Acknowledgements The authors acknowledge partial financial support received from Junta de Castilla y León under project VA11/99 and by MCyT under Project AYA2000-1787. The authors are deeply gratefully to Dr. Alan C. Hindmarsh for his valuable comments in the elaboration of this paper.
References [1] A.W. Bush, Perturbation Methods for Engineers and Scientists, CRC Press, Boca Raton, FL, 1992. [2] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, Families of Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 235–250. [3] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, High-order embedded Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 423–430. [4] J.M. Farto, A.B. González, P. Martín, An algorithm for the systematic construction of solutions to perturbed problems, Comput. Phys. Comm. 111 (1998) 110–132. [5] J.M. Ferrándiz, A general canonical transformation increasing the number of variables with application to the two-body problem, Celestial Mech. 41 (1988) 345–357. [6] A.B. González, P. Martín, J.M. Farto, A new family of Runge–Kutta type methods for the numerical integration of perturbed oscillators, Numer. Math. 82 (1999) 635–646. [7] P. Martín, J.M. Ferrándiz, Multistep numerical methods based on the Scheifele G-functions with application to satellite dynamics, SIAM J. Numer. Anal. 34 (1997) 359–375. [8] P. Martín, J.M. Farto, Improved numerical integration of perturbed oscillators via average, Appl. Math. Comput. 99 (1999) 129–139. [9] U. Kirchgraber, An ODE-solver based on the method of averaging, Numer. Math. 53 (1989) 621–652. [10] S. Novo, J. Rojo, Some remarks on an ODE-solver of Kirchgraber, Numer. Math. 61 (1992) 261–264. [11] J.M. Ferrándiz, M.E. Sansaturio, J.R. Pojman, Increased accuracy of computations in the main satellite problem through linearization methods, Celestial Mech. Dynam. 53 (1992) 347–363. [12] G. Scheifele, On numerical integration of perturbed linear oscillating system, Z. Angew. Math. Phys. 22 (1971) 186–210. [13] E.L. Stiefel, G. Scheifele, Linear and Regular Celestial Mechanics, Springer, New York, 1971.
Increasing the accuracy in the numerical integration of perturbed oscillators with new methods Ana B. González ∗ , Pablo Martín, Amelia García, José M. Farto Departamento de Matemática Aplicada a la Ingeniería, E.T.S. de Ingenieros Industriales, Universidad de Valladolid, Paseo del Cauce s/n, 47011 Valladolid, Spain Received 10 October 2001
Abstract An important problem in celestial mechanics is the precise long-term prediction of satellite orbits. This is a difficult problem due to different factors. When the problem is formulated in BF or KS variables its structure is of the form of perturbed oscillators. The effect of the choice of the frequency in the integration of perturbed oscillators leads us to construct new methods which exhibit a very good behaviour. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction The precise long-term prediction of Earth satellite orbits is an important problem in celestial mechanics. Numerical solutions can move or spiral away from the exact solution due to Lyapunov or orbital instability, and the integration remains valid only for a span of time that could be too short. It is very difficult to increase the accuracy, since the problem is affected by different factors, ranging from deficiencies in the force model, to uncertainties in the initial values, to the simple accumulation of errors during the integration. We are concerned here only with integration errors, and in order to minimize these, we transform the equations of motion to others suitable for the numerical integration by means of changes of variables such as the linearizing transformations as KS (Stiefel and Scheifele [13]) or BF variables (Ferrándiz [5]). In these variables the equations of motion are reduced to a system of perturbed oscillators. In 1971 Scheifele [12] obtained a refinement of the Taylor method based on his G-functions for the integration of perturbed oscillators. But, the Scheifele G-function method has the disadvantage that it can only be applied in very particular cases, due to the complexity of the preliminary calculations that it requires. Recently, new methods based on the Scheifele G-functions (SMF codes of Martín and * Corresponding author.
E-mail address: [email protected] (A.B. González). 0168-9274/01/$22.00 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 1 2 1 - 0
296
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
Ferrándiz [7] and RKGM schemes of González et al. [6]) have been developed to avoid such problems. But most of the special methods require a previous knowledge of a good approximation to the frequency. The effect of this frequency in the integration of perturbed oscillators is very important as we can see in Kirchgraber [9], Novo and Rojo [10] and Martín and Farto [8]. This fact leads us to formulate new purely Scheifele-based methods which inherit its good behaviour in order to improve by selecting a suitable frequency. In particular, we construct a fourth order scheme which shows very good results in the numerical integration of perturbed oscillators. The choice of a suitable approximation to the frequency for this method, leads to an important reduction in error. Finally, the new algorithm is applied to the propagation of artificial Earth satellite orbits once the corresponding equations of motion have been reduced to oscillatory form by using linearizing variables.
2. The solution of a perturbed oscillator Let us consider the model problem ˙ y¨ + ω2 y = ε f (t, y, y),
y(0) = y0 ,
y(0) ˙ = y˙0 ,
(1)
where ε is a small parameter. If we suppose that the function f (t, y(t), y(t)) ˙ admits an expansion of the form
˙ = f t, y(t), y(t)
∞
(j ) t
f0
j =0
∞
,
the j th derivative of f at t = 0), then the initial value problem (1) can be
(j ) t
j
f0
j =0
j!
(j )
(we have denoted by f0 written as y¨ + ω2 y = ε
j
j!
,
y(0) = y0 ,
y(0) ˙ = y˙0 .
(2)
The solution of (2) can be obtained in the usual way, building the solution of the homogeneous equation with the given initial conditions and adding to it the solution of the non-homogeneous equation in which the solution and its derivative at t = 0 are cancelled out. The second solution can be calculated by solving the initial value problems y¨j + ω2 yj =
tj , j!
yj (0) = y˙j (0) = 0,
j 0,
(3) (j )
and linearly combining their solutions with the coefficients εf0 . Scheifele defines his G-functions as the solutions to the problems (3) with the notation Gj (t) = yj −2 (t) for j 2. The solution to the homogeneous problem will be obtained by simply combining the solutions of the problems y¨ + ω2 y = 0, y¨ + ω2 y = 0,
y(0) = 1, y(0) = 0,
y(0) ˙ = 0, y(0) ˙ = 1,
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
297
which Scheifele denotes as G0 (t) and G1 (t), respectively, and suggests a systematic way to obtain the solution of the problem (1) which is given by: Y (t) = y0 G0 (t) + y˙0 G1 (t) + ε
∞
(j )
f0 Gj +2 (t).
(4)
j =0
The Scheifele G-functions have the following expressions in terms of the elementary functions
j −1
k t 2k (−1)j −ω2 , G2j (t) = 2j cos ωt − ω (2k)! k=0
(−1)j G2j +1 (t) = ω2j
j −1
k t 2k+1 1 sin ωt − −ω2 , ω (2k + 1)! k=0
j 0,
(5)
or in terms of powers of t as follows Gj (t) =
∞
−ω2
k
k=0
t 2k+j , (2k + j )!
j 0.
(6)
See Scheifele [12]. Let h a real number. If we truncate expression (4) up to j = p − 2, we obtain an approximation of order p to the solution at the point t = h. In spite of its good behaviour, the Scheifele G-function method has the great disadvantage that it can only be applied in very particular cases, due to the complexity of the preliminary calculations that it requires. This problem was solved by Martín and Ferrándiz [7] in a drastic way, through the conversion of the G-function method into a multistep scheme, or more recently, by González et al. [6] who construct a new one-step approximation to the solution. But most of the special methods need a good approximation to the frequency of the problem. We shall analyze in the following section the effect of such frequency on the numerical solution provided by some special methods.
3. The effect of the frequency on the solution of a perturbed oscillator It is well known that the choice of a good approximation to the actual frequency of the solution, improves notably the accuracy of some kinds of numerical methods when integrating weakly perturbed oscillators (see Kirchgraber [9], Novo and Rojo [10] and Martín and Farto [8]). We shall discuss this phenomenon. To this end, we consider the simple and widely used test problem y¨ + y = εy 3 ,
y(0) = 1,
y(0) ˙ = 0.
(7)
The true frequency of the solution can be easily obtained by means of perturbation techniques (Lindstedt method, Bush [1]) as a power series of ε. A first order approximation to this frequency is 3 ω¯ = 1 − ε. 8 In this way, we can rewrite problem (7) as
y¨ + ω¯ 2 y = εy 3 − 1 − ω¯ 2 y,
y(0) = 1,
y(0) ˙ = 0.
(8)
298
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
Fig. 1. y¨ + y = εy 3 , y(0) = 1, y(0) ˙ = 0, ε = 10−3 , h = 0.1.
While for the integration with a Taylor method, the use of one or other form is irrelevant, when integrating with a Scheifele method, the selection of the expression (7) or (8) is very important. This last method depends explicitly on the frequency of the unperturbed oscillator considered. Fig. 1 shows the integration of the test problem with a classical embedded pair RKN4(3) (Dormand et al. [2]), and the schemes depending on the G-functions (and therefore on the frequency), the multistep code SMF of Martín and Ferrándiz [7] and the one-step scheme RKGM of González et al. [6]. We can observe that the results obtained with the SMF method when considering as frequency ω = 1 − 38 ε are better than those corresponding to ω = 1. The right choice leads to a remarkable stabilization of the error. These results have been improved thanks to the elimination of the resonance that appeared. But surprisingly, the RKGM methods do not seem to inherit the good behaviour to the change of the frequency of the Scheifele scheme. The RKGM scheme do not show very much reduction in error with the correction to the frequency. In fact, the improvement achieved is irrelevant. A further analysis is necessary to understand this behaviour. This will be done in the next section.
4. A new approximation for the perturbed oscillator The RKGM method
y1 = y0 G0 (h) + y˙0 G1 (h) + ε
p−2
s
j =0
i=1
y˙1 = y˙0 G0 (h) − ω y0 G1 (h) + ε 2
ki = f ci h, y˜i ,
i = 1, . . . , s,
bj i ki Gj +2 (h),
p−2
s
j =0
i=1
bj i ki Gj +1 (h),
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
y˜1 = y0 ,
y˜i = y0 + ci hy˙0 + h2
i−1
299
aij εkj − ω2 y˜j ,
j =1
appears as a kind of mixed method. While the final approximation to the solution is a Scheifele-like formula, the internal stages ki are obtained using arguments based on Taylor series to approximate y(ci h). But y(t) is the solution of a perturbed oscillator, and thus it seems more natural to approximate the exact solution at the points ci h using the solution (4) provided by Scheifele. This characteristic may be the cause of the unexpected response of the RKGM method to the change of frequency. We should not want to relinquish to the excellent properties of the RKGM method when integrating satellite orbits. If they could be arranged so as to avoid the deficiency noted, then the performance could be improved. Unfortunately, the problem of RKGM methods is structural. For this reason, the development of a new approximation preserving the good properties of the RKGM method, and purely Scheifele-based, is worth the effort. In fact, this can be done as in the following example. We begin by considering the Scheifele method of order 4, Y1 = y0 G0 (h) + y˙0 G1 (h) + ε
2
(j )
f0 Gj +2 .
j =0
We shall try to approximate the first two derivatives of f as follows (j )
f0 ≈ βj 1 g1 + βj 2 g2 + βj 3 g3 ,
j = 0, 1, 2,
with βj i real numbers and gi evaluations of f at the points t¯ = 0, t¯ = h/2 and t¯ = h. Obviously, y(0) = y0 , then we take g1 = f (0, y(0)) = f (0, y0 ). Now we need to approximate y(h/2) and y(h). From (4) we have y(h/2) = y0 G0 (h/2) + y˙0 G1 (h/2) + ε
∞
(j )
f0 Gj +2 (h/2),
(9)
j =0
then y(h/2) ≈ y¯2 = y0 G0 (h/2) + y˙0 G1 (h/2) + εg1 G2 (h/2), thus we take g2 = f (t¯2 , y¯2 ). A similar idea for y(h) allows us to consider y¯3 = y0 G0 (h) + y˙0 G1 (h) + εg2 G2 (h), and g3 = f (t¯3 , y¯3 ). In order to obtain the parameters βj i , it is obvious that β01 = 1 and β02 = β03 = 0. For the values ˙ y )|t =0 and the following β1i , i = 1, 2, 3, we estimate the agreement between terms of f0(1) = (ft + yf Taylor series about h = 0 β11 g1 + β12 g2 + β13 g3 = (β11 + β12 + β13 )f + ft (β12 /2 + β13 ) + fy y˙0 (β12 /2 + β13 ) h
+ ft t (β12 /8 + β13 /2) + fty y˙0 (β12 /4 + β13 ) + fy εf − ω2 y (β12 /8 + β13 /2)
+ fyy y˙ 2 (β12 /8 + β13 /2) h2 + O h3
300
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
obtaining the system: β11 + β12 + β13 = 0, β + 2β13 = 2/ h, (10) 12 β12 + 4β13 = 0, whose solution is β11 = −3/ h, β12 = 4/ h, β13 = −1/ h. For the values β2i , i = 1, 2, 3, we estimate the agreement between terms of f0(2) = (ft t + 2fty y˙ + y˙ 2 fyy + fy (εf − ω2 y))|t =0 and the Taylor series about h = 0 β21 g1 + β22 g2 + β23 g3 = (β21 + β22 + β23 )f + ft (β22 /2 + β23 ) + fy y˙0 (β22 /2 + β23 ) h
+ ft t (β22 /8 + β23 /2) + fty y˙0 (β22 /4 + β23 ) + fy εf − ω2 y (β22 /8 + β23 /2)
+ fyy y˙ 2 (β22 /8 + β23 /2) h2 + O h3
obtaining the system: β21 + β22 + β23 = 0, β + 2β23 = 0, 22 β22 + 4β23 = 8/ h2 ,
(11)
whose solution is β21 = 4/ h2 , β22 = −8/ h2 , β23 = 4/ h2 . The final approximation to the solution takes the form y1 = y0 G0 (h) + y˙0 G1 (h) + ε
2 3
βj i gi Gj +2 (h),
j =0 i=1
y˙1 = y˙0 G0 (h) − ω2 y0 G1 (h) + ε
2 3
βj i gi Gj +1 (h).
j =0 i=1
Fig. 2. y¨ + y = εy 3 , y(0) = 1, y(0) ˙ = 0, h = 0.1, ε = 10−3 . Frequency for RKGM and SMF ω = 1 − (3/8)ε.
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
301
From the construction and property (6) we have (j )
(βj 1 g1 + βj 2 g2 + βj 3 g3 )Gj +2 = f0 Gj +2 + O h5 , (βj 1 g1 + βj 2 g2 + βj 3 g3 )Gj +1 =
(j ) f0 Gj +1
+O h , 4
j = 1, 2, j = 1, 2.
Therefore, at least we achieve order four in the approximation to the solution and order three in the derivative. But a simple Taylor series of the approximation to the derivative with these particular coefficients shows that we also achieve order four in this case. As a result of removing all the Taylor inheritance from the new method, we can see in Fig. 2 its good response to the frequency modification technique. Thus, we have reached our target, and also we have constructed a new method with better behaviour than the RKGM scheme. In a variety of satellite problems, the true frequency is very difficult to approximate, or such approximation simply can not be obtained. In these cases we have obtained an improved numerical method. This fact justifies the construction of the new formula. 5. Reducing the error growth in long-term computation of satellite orbits This section is devoted to a study of the behaviour of the new approximation constructed in the previous section when integrating the satellite problem. To this end, we consider the problem expressed in linearizing, using KS (Stiefel and Scheifele [13]) or BF (Ferrándiz [5] and Ferrándiz et al. [11]) variables. In KS variables the equations of motion become 1 |Y |2 , ω˙ = R, (12) t˙ = y¨i + yi = Ri , i = 1, 2, 3, 4, 4 ω √ where Y = (y1 , y2 , y3 , y4 ) is the vector of the KS coordinates, ω = 2p0 , p0 stands for the opposite of the total energy and Ri and R denote the perturbation terms. The dot denotes the derivative with respect to an eccentric anomaly-like variable s. The problem in focal variables (BF variables) is reduced to four perturbed harmonic oscillators with unit frequency. The coordinates of the basic set of focal variables are three components (y1 , y2 , y3 ) of the direction vector of the particle and the inverse u of the radial distance. The system is completed with the addition of an equation which provides the variation of the angular momentum c (which has a stabilizing effect on the solutions) and with the equation corresponding to the change of time. Then the equations of motion can be described as y˙i y¨i + yi = Qi − c˙ , i = 1, 2, 3, c u µ u¨ + u = 2 + Q − c˙ , (13) c c 1 F t˙ = 2 , c˙ = , c cu where µ is the reduced mass, Qi , Q and F stand for the corresponding disturbing terms, and the dot denotes the derivative with respect to a generalized true anomaly variable τ . We have considered some orbits with different initial values and different sets of perturbations and we have investigated the error produced during the numerical propagation of the orbits over several tens of revolutions.
302
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
We consider an almost periodic equatorial orbit with the coefficient J2 due to the oblateness as the perturbation parameter. We have neglected terms of the order of J22 . The eccentricities considered are e = 0 and e = 0.99. Thus, for BF variables the problem can be written as µ J2 2 + 12 u, (14) c2 c2 where µ/c2 = 20/21 and J2 /c2 = (10/21) × 10−3 for e = 0 and µ/c2 = 100/20895 and J2 /c2 = (50/20895) × 10−3 for e = 0.99. We present efficiency curves for each problem, where the common logarithm of the maximum global error over all steps and variables is plotted against the number of evaluations of the function required. These graphics give a good indication of the efficiency of the considered methods when they are carefully implemented. To this end, we rewrite the final approximation of the RKGM method as y¨i + yi = 0,
i = 1, 2, 3,
y1 = y0 G0 (h) + y˙0 G1 (h) + ε
u¨ + u =
p−2 s i=1
y˙1 = y˙0 G0 (h) − ω y0 G1 (h) + ε 2
bj i Gj +2 (h) ki ,
j =0
p−2 s i=1
bj i Gj +1 (h) ki ,
j =0
and we compute the coefficients of the internal stages ki at the beginning of the integration, taking into account the expression (6) for the G-functions. We use the same idea for the new codes. We compute at the beginning of the integration the coefficients of gi and the G-functions evaluated at ci h. Thus, the overhead per function evaluation is approximately the same for these two methods and a fixed step size RKN scheme. On the other hand, when considering problems with much more costly right-hand side functions, the overhead costs become much less significant. Fig. 3 shows the results obtained when integrating 1000 revolutions of the problem in BF variables and eccentricities e = 0 and e = 0.99. In this case we have a reference solution (see Farto et al. [4])
Fig. 3. BF variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 1000 revolutions. Left: e = 0. Right: e = 0.99.
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
303
and the error is computed at each point of the numerical integration. The suitable frequency for e = 0 is ω = 0.9945578231292517 and for eccentricity e = 0.99 is ω = 0.9986257476147703. We have denoted by New1 and New2 the results of the new method with frequencies ω = 1 and the corrected frequency, respectively, and similarly for SMF1 and SMF2. We have used as step sizes h = 1/2i , i = 1, . . . , 7, for the RKGM and new methods and h = 1/2i , i = 1, . . . , 8, for the SMF scheme. The tolerances used are 10−3 , . . . , 10−9 for the embedded pair. In both cases, low and high eccentricity, the best results correspond to the new approximation. Figs. 4 and 5 correspond to the case of KS variables. In this case, an expression for the solution is not available. Thus we have obtained at perigees and apogees a reference numerical solution, in which we have also computed the numerical methods. For this reference solution, we have considered the RKN
Fig. 4. KS variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 20 revolutions. Left e = 0. Right: e = 0.5.
Fig. 5. KS variables. J2 perturbation. Inclination 0. Perigee distance 1.05RE . Efficiency curves when integrating 20 revolutions. e = 0.99.
304
A.B. González et al. / Applied Numerical Mathematics 41 (2002) 295–304
method of order eight constructed by Dormand et al. [3] and implemented with fixed step size. We have used arithmetic precision of 32 significant figures and step size h = π/29 . For e = 0 and e = 0.5, there are six runs corresponding to h = π/23, . . . , π/28 for the RKGM method and h = π/22 , . . . , π/27 for the new scheme, and tolerances 10−3 , . . . , 10−6 for the embedded pair. For e = 0.99, there are runs with h = π/23 , . . . , π/27 for the new scheme and the same tolerances as before for the RKN. For the first two eccentricities the best results correspond to the new method, while for e = 0.99, new and RKGM schemes present a similar behaviour. This fact is not surprising, since the equations in KS variables are singular at e = 1 (see Ferrándiz [5]). Acknowledgements The authors acknowledge partial financial support received from Junta de Castilla y León under project VA11/99 and by MCyT under Project AYA2000-1787. The authors are deeply gratefully to Dr. Alan C. Hindmarsh for his valuable comments in the elaboration of this paper.
References [1] A.W. Bush, Perturbation Methods for Engineers and Scientists, CRC Press, Boca Raton, FL, 1992. [2] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, Families of Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 235–250. [3] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, High-order embedded Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 423–430. [4] J.M. Farto, A.B. González, P. Martín, An algorithm for the systematic construction of solutions to perturbed problems, Comput. Phys. Comm. 111 (1998) 110–132. [5] J.M. Ferrándiz, A general canonical transformation increasing the number of variables with application to the two-body problem, Celestial Mech. 41 (1988) 345–357. [6] A.B. González, P. Martín, J.M. Farto, A new family of Runge–Kutta type methods for the numerical integration of perturbed oscillators, Numer. Math. 82 (1999) 635–646. [7] P. Martín, J.M. Ferrándiz, Multistep numerical methods based on the Scheifele G-functions with application to satellite dynamics, SIAM J. Numer. Anal. 34 (1997) 359–375. [8] P. Martín, J.M. Farto, Improved numerical integration of perturbed oscillators via average, Appl. Math. Comput. 99 (1999) 129–139. [9] U. Kirchgraber, An ODE-solver based on the method of averaging, Numer. Math. 53 (1989) 621–652. [10] S. Novo, J. Rojo, Some remarks on an ODE-solver of Kirchgraber, Numer. Math. 61 (1992) 261–264. [11] J.M. Ferrándiz, M.E. Sansaturio, J.R. Pojman, Increased accuracy of computations in the main satellite problem through linearization methods, Celestial Mech. Dynam. 53 (1992) 347–363. [12] G. Scheifele, On numerical integration of perturbed linear oscillating system, Z. Angew. Math. Phys. 22 (1971) 186–210. [13] E.L. Stiefel, G. Scheifele, Linear and Regular Celestial Mechanics, Springer, New York, 1971.