Runge–Kutta Pairs of Orders 6(5) with Coefﬁcients Trained to Perform Best on Classical Orbits

: We consider a family of explicit Runge–Kutta pairs of orders six and ﬁve without any additional property (reduced truncation errors, Hamiltonian preservation, symplecticness, etc.). This family offers ﬁve parameters that someone chooses freely. Then, we train them in order for the presented method to furnish the best results on a couple of Kepler orbits, a certain interval and tolerance. Consequently, we observe an efﬁcient performance on a wide range of orbital problems (i.e., Kepler for a variety of eccentricities, perturbed Kepler with various disturbances, Arenstorf and Pleiades). About 1.8 digits of accuracy is gained on average over conventional pairs, which is truly remarkable for methods coming from the same family and order.


Introduction
The initial value problem (IVP) is with x 0 ∈ R, y, y ∈ R m and f : R × R m → R m . Runge-Kutta (RK) pairs are amongst the most popular numerical methods for addressing (1). They are characterized by the following Butcher tableau [1,2]: with b T ,b T , c ∈ R s and A ∈ R s×s . Then, the method shares s stages, and when c 1 = 0 and A is strictly lower triangular, it is evaluated explicitly. The approximated solution steps from (x n , y n ) to x n+1 = x n + h n by producing two estimations for y(x n+1 ). Namely, y n+1 andŷ n+1 , given by for i = 1, 2, · · · , s. These two approximations y n+1 and y n+1 are of algebraic orders p and q < p respectively. Thus, a local error estimation is formed in every step and is combined in an algorithm for changing the step size: where t is a tolerance given by the user. When n < t, the above formula is used for the new step forward. In reverse, we also use it, but the solution is not advanced and h n+1 is a new version of h n . Details can be found in [3]. As an abbreviation, these methods are named RKp(q) pairs. Runge-Kutta methods were introduced back in the late 19th century [4,5]. After 1960, RK pairs appeared. Fehlberg gave the first celebrated such pairs of orders 5(4), 6(5) and 8 (7) [6,7]. Dormand and Prince followed in the early 1980s [8,9]. Our group has also presented a series of successful RK pairs [10][11][12][13].
Runge-Kutta pairs are suited for efficiently solving almost every non-stiff problem of the form (1). The variety of pairs is explained by the accuracy required. Thus, when less accuracy is required, the lowest RK pairs are more efficient. In contract, for stringent accuracies at quadruple precision, a high-order pair should be chosen [14].
Here we focus on RK6(5) pairs, which are preferred for moderate to higher accuracies. We are especially interested in problems (1) that resemble Kepler-like orbits. Thus, we will propose a particular RK6(5) pair for addressing this type of problem.

Producing Runge-Kutta Pairs of Orders 6(5) and Training Their Coefficients
Runge-Kutta pairs of orders six and five are amongst the most frequently used. The coefficients have to satisfy 54 order conditions. Thus, families of solutions have been discovered over the years. Here we chose the Verner-DLMP [15,16] family, which has the advantage of being solved linearly. Then, we freely chose the coefficients c 2 , c 4 , c 5 , c 6 , c 7 andb 9 . Pairs from this family have been proven to perform most efficiently in various classes of problems [17].
We proceed by explicitly evaluating the remaining coefficients. The algorithm is discussed in [10] and is given as a Mathematica [18] package in the Appendix A.
Although s = 9, the family spends only eight stages per step since the ninth stage is used as the first stage of the next step. This property is called FSAL (first stage as last).
The next question to be answered is regarding how to select the free parameters. Traditionally we try to minimize the norm of the principal term of the local truncation error. That is, the coefficients of h 7 in the residual of Taylor error expansions corresponding to the sixth-order method of the underlying RK pair.
We intend to derive a particular RK6(5) pair belonging to the family of interest here. The resulting pair has to perform best on Kepler orbits and other problems of this nature. Thus, we concentrate on the particular orbit 1 y = 3 y, 2 y = 4 y, 1−ecc T and the theoretical solution In the above, v = ecc · sin(u) + x, ecc is the eccentricity, and the components of y are denoted by the left superscript. They should not be confused with y 1 = 1 y 1 , 2 y 1 , 3 y 1 , 4 y 1 T , y 2 = 1 y 2 , 2 y 2 , 3 y 2 , 4 y 2 T , y 3 , · · · , which represent the vectors approximating the solution at x 1 , x 2 , x 3 , · · · . This problem can be solved with an RK6(5) pair from the family we are interested in here. After a certain run, we recorded the number f ev of function evaluations (stages) needed and the global error ge observed over the mesh (grid) in the interval of integration. Then, we formed the efficiency measure u = f ev · ge 1/6 . (2) Running DLMP6(5) twice for Kepler we obtained the efficiency measuresû 1 andû 2 as reported in Table 1. For example, we needed 1121 stages in order to achieve a global error of 2.14 · 10 −6 when running Kepler for ecc = 0, x end = 10π and tol = 10 −7 . Thus, we obtainedû 1 = 1123 · 2.14 · 10 −6 ≈ 127.22, as reported in Table 1. Analogously, for a second run with ecc = 0.6, x end = 20π and tol = 10 −11 we observedû 2 = 833.27.
Let us suppose that any new pair NEW6(5) furnishes corresponding efficiency measuresũ 1 andũ 2 for the same runs. Then, as a fitness function we form the sum and try to maximize it. That is, the fitness function is actually two whole runs of initial value problems. The value u changes according to the selection of the free parameters c 2 , c 4 , c 5 , c 6 , c 7 andb 9 .
The original idea is based on [19]. For the minimization of u we used the differential evolution technique [20]. We have already tried this approach and obtained some interesting results in producing Numerov-type methods for integrating orbits [21]. In this latter work we trained the coefficients of a Numerov-type method on a Kepler orbit. We observed very pleasant results over a set of Kepler orbits as well as other known orbital problems. The optimization furnished six values for the free parameters. They are included with the rest of the coefficients in the resulting pair NEW6(5) presented in Table 2.
In conclusion, no extra property seemed to hold. The pair given in Table 2 does not possess any interesting properties. It is difficult to believe a special performance could be obtained after seeing its traditional characteristics.

Numerical Tests
We tested the following pairs chosen from the family studied above.
DLMP6 (5) is the best representative of conventional RK pairs. Everything else presented until now is hardly more efficient [10]. Both pairs were run for tolerances of 10 −5 , 10 −6 , · · · , 10 −11 , and the efficiency measures (2) were recorded for each one. We set NEW6(5) as the reference pair. Then we divided each efficiency measure of DLMP6 (5) with the corresponding efficiency measure of NEW6 (5). Numbers greater than 1 indicate that NEW6(5) is more efficient. Thus, we can interpret the number 1.1 as DLMP6(5) being 0.1 = 10% more expensive than NEW6(5) while an entry of 2 means that DLMP6(5) is 100% more expensive (i.e., has twice the cost for achieving the same accuracy).
The problems we tested were as follows.

The Kepler problem
This problem is explained above. We ran it for five different eccentricities (i.e., ecc = 0, 0.2, 0.4, 0.6, 0.8), while we recorded the efficiency measures using the endpoint errors for x end = 10π and x end = 20π.
We transformed this problem into a system of four first-order equations and solved for x end = 10π and x end = 20π. After recording the endpoint errors and the costs, we present the efficiency measures ratios of DLMP6(5) vs. NEW6(5) for perturbed Kepler in Tables 5 and 6.

The Arenstorf orbit
Another interesting orbit describes the stable movement of a spacecraft around Earth and Moon ( [22], pg. 129). We also transformed this problem to a system of four first-order equations and solved it to x A and 2x A . After recording the endpoint errors and the costs we present the efficiency measures ratios of DLMP6(5) vs. NEW6(5) for Arenstorf in Table 7.

Conclusions
This paper is concerned with training the coefficients of a Runge-Kutta pair for addressing a certain kind of problem. We concentrated on problems with Kepler-type orbits and an extensively studied family of Runge-Kutta pairs of orders six and five. After optimizing the free parameters (coefficients) in a couple of runs on Kepler orbits, we concluded to a certain pair. This pair was found to outperform other representatives from this family in a wide range of relevant problems.

Conflicts of Interest:
The authors declare no conflicts of interest.