1. Introduction
A system of ordinary differential equations of the form
with
,
and
, is called Initial Value problem (IVP).
For addressing the above problem (
1), we usually try Runge–Kutta (RK) pairs, which is perhaps the most popular choice among the numerical methods. These pairs are compactly presented by the Butcher tableau [
1,
2] that follows,
with
and
This method spends
s function evaluations per step that are evaluated explicitly when
D is a strictly lower triangular. The estimate of the solution steps from
to
after producing two approximations for
Namely,
and
, given by
and
with
for
These two estimations
and
are, respectively, of algebraic orders
and
. Thus, in every step, we form the quantity
and it is supposed to be an estimate of the local error. Actually, it is used in a step-size changing algorithm
where
is a small enough positive tolerance set by the user. Whenever
we use the above formula for setting the length of the next step forward, i.e., from
to
. On the contrary, the same formula is used, but then
is changed with
since we have to estimate a new smaller step size for advancing from the current point
to a new point
. Details in this subject are given in [
3]. An abbreviation for these methods is the terminology
pairs.
Runge was the first to present a second order RK method by combining a sequence of Euler formulas [
4]. Some years later, Kutta managed to construct a four stages fourth order method [
5]. Since the system of order conditions required to be solved (see [
6]) is too complicated, there was a rather slow progress in the issue. Nyström showed the correct coefficients of a six stages method of fifth order [
7]. Higher order methods appeared in the 50s. Huta gave a sixth order method at a cost of eight stages per step [
8]. RK pairs were introduced in the 60s. It was E. Fehlberg who constructed at first a series of celebrated pairs of orders 5(4), 6(5) and 8(7), [
9,
10]. His pioneering work followed in the early 80s by Dormand and Prince who presented their famous pairs [
11,
12]. Our research group has also derived a series of interesting RK pairs [
13,
14,
15,
16].
Almost every non-stiff problem of the form (
1) may be solved effectively with RK pairs. The accuracy required explains the wide range of pairs. As a result, the lowest RK pairs are more efficient the lower the accuracy on demand. A high order pair, on the other hand, should be chosen for demanding accuracies at quadruple precision [
17].
We’ll concentrate on RK5(4) pairs, which are best for intermediate accuracies. Problems (
1) with structure similar to Keplerian forms are of particular interest to us. As a result, we shall offer a specific RK5(4) pair for dealing with this type of cases.
2. Families of Runge–Kutta Pairs of Orders 5(4)
Runge–Kutta pairs of algebraic orders five and four are perhaps the most famed. Their coefficients have to satisfy 25 order conditions. As a result, over time, families of solutions have been developed. The minimum number of stages required to construct a fifth-order RK method or a 5(4) pair is known to be six [
1].
Usually these pairs are constructed according to various types of simplifications. After these simplified assumptions applied to the original equations (order conditions), we arrive at a simpler system where some coefficients remain as free parameters. The values of these free parameters were chosen by Fehlberg for minimization the truncation error coefficients of the pair’s fourth-order formula. After considerable numerical testing, Shampine in [
18] concluded that propagating the higher-order solution (i.e., here, the fifth order) of a pair is desirable from a numerical standpoint. This technique is called Local Extrapolation.
Dormand and Prince went on to expand a family of fifth-order methods already presented by Butcher in [
1]. They suggested a family that borrows the first function evaluation from the next step to embed a fourth-order method within a fifth-order method at no additional cost [
11]. For a long time, the best fifth-order pair was commonly considered to be DP5(4), an individual pair in their four-parameter family with small truncation error coefficients of its fifth-order formula. Later on, Papakostas and Papageorgiou extended this family by adding a parameter more [
19]. The later authors presented a pair with minimal truncation coefficients and claimed that its performance was in average about
more efficient than DP5(4) in the standard IVP test set named DETEST [
20].
Ten years ago, Tsitouras presented a new family that makes use of the least simplifying assumptions possible [
15]. Only the assumptions
with
,
the identity matrix, were made there.
The common set of simplifications
is not used. In the above, we set as
the component-wise multiplication among the elements of
q. Generally, in the following, by
we mean the component by component multiplication among matrices with the same number of rows.
Here, we will make use of the five parameters family introduced in [
19]. The choice is justified, since this family of solutions is a superset of the solutions given by [
11]. Besides, the derivation of the coefficients is explicit and therefore cheap. Meanwhile, the family given in [
15], although groundbreaking, is implicit. Thus, the later technique of solutions is somewhat expensive for our goals here.
In the following, we describe the derivation of the coefficients according to the guidelines given in [
19].
Notice that, at first, and for .
Set arbitrary values for has to be different from each other and not equal to 1.
Solve for
Set
Solve for
Solve for and By is meant the th component of the underlying vector.
Then set and
Finally, the FSAL (First Stage As Last) property holds, and
Because the seventh stage is used as the initial stage of the following step, even if
, the family only spends six stages every step.
The challenge now is how to choose the free parameters. The norm of the principal coefficients of the local truncation error is traditionally minimized. i.e., the terms of
in the residual of Taylor error expansions corresponding to the fifth order method of the underlying RK pair [
13].
Since our concern here is to deal with Keplerian type orbits, we shall try another approach and train these free parameters in order for the resulting pair to perform best in our problems of interest.
3. Comparison of the Results among Pairs
Comparison of the results observed by various pairs is an interesting issue. It is of crucial importance in our present work. So, let us work on Kepler problem. This has the following form:
with
,
and the theoretical solution
In the above, , is the eccentricity, and the the left superscript denotes the components of x. They shall not be confused with , that correspond to the vectors approximating the solution at .
When running DP5(4) for
and tolerances
,
,
, we get the results summarized in
Table 1. There we recorded the error observed at the end-point as long as the function evaluations are taken by the pair.
Then, we try T5(4) pair given in [
15] and get the results presented in
Table 2.
Which pair is more efficient? How can we derive this? Of course, we may not simply check the errors after each tolerance since the number of stages differ. We want to eliminate the presence of tolerances. They do not offer something in comparing the results. We proceed in a way similar to that proposed in [
20].
Thus, we form a linear least squares fit on the logarithms of data (i.e., stages and errors). Then, we arrive at a slope that justifies the order of the method. These parallel slopes for DP5(4) and PP5(4) can be seen in
Figure 1. See ([
2], p. 214) for details in the issue.
Now, we are able to compare the results of the pairs in the selected problem.
In
Table 3, we present the functions evaluations that might have been spend for achieving various errors by both pairs. These numbers come after the linear fit as explained in [
20], while the constant ratio is explained by the parallel slopes among methods of the same order. The later observation might not be present for every problem. Some special characteristics of various problems may exploited differently by various pairs of the same orders and thus may give steeper slopes than expected. Thus, the numbers in the second column of
Table 3 are evaluated as
while the stages in the third column can be found by the relation
In conclusion, DP5(4) is considered to be 13–14% cheaper than T5(4) in the problem under consideration. After this, we may ask for free parameters that may produce a new pair that is cheaper in a wider class of Keplerian type orbits.
Both methods outperformed their order expectations for this problem. According to theoretical aspects, a fifth order method is expected to perform as
since the global error behaves as
, then, with
h a constant step through the integration, ([
21], p. 162).
4. Performance of Pairs in a Wide Set of Keplerian-like Problems
We indent to construct a specific RK5(4) pair from the above mentioned family. The resulting pair has to perform best on Keplerian type problems. Thus, we have chosen the following problems.
1. The Kepler problem
This problem has already been discussed above. We ran it for five different eccentricities (i.e., ), while we recorded the stages used and the endpoint errors for . These problems are referred to with numbers 1–5, respectively, in the tables with the results below.
2. The perturbed Kepler
The Schwarzschild potential is used to explain the motion of a planet according to Einstein’s general relativity theory. The equations are:
and the analytical solution is
We solved this problem through
, by converting it into a system of four first-order equations. We run this problem for
and referred with numbers 6–10, respectively, in the tables with the results below.
3. The Arenstorf orbit
Another fascinating orbit depicts a spacecraft’s stable travel around the Earth and moon ([
21], p. 129).
with
initial values
and with
,
the solution is periodic.
We also transformed this problem to a system of four first-order equations and solved it to and . Similarly, we recorded the endpoint errors and the costs and used numbers 11 and 12 for these two runs in the tables with the numerical results.
4. The Pleiades
Finally, we considered the problem “Pleiades” as given in ([
21], p. 245).
with
We set
. We solved the problem by transforming it into a system of fourteen first-order equations once more. As end points, we used to
and 4. After estimating the solution with
Mathematica [
22] and quadruple precision, we recorded the endpoint errors made by the RK5(4) pairs. Similarly we recorded the endpoint errors and the costs and used numbers 13 and 14 for these two runs in the tables with the numerical results.
In conclusion, we set 14 problems (i.e., five Keplerian, five perturbed Kepler, two Arenstorf orbits and two Pleiades) run for 7 tolerances each (i.e., ).
At first we record the corresponding ratios (as those recorded in the rightmost column of
Table 3) of DP5(4) vs. T5(4) in
Table 4. In the first column, we list the expected errors. In the last row we see the average performances over all runs (i.e., tolerances used) for each problem. Numbers greater than 1 ar in favor of the second pair. The overall average of these means is
, which means that DP5(4) was in average about
more expensive than T5(4) in all these
runs.
Thus, our question is to find free parameters in the algorithm given in the previous section that furnish a pair doing much better than as it can be.
5. Training the Free Parameters in a Wide Set of Keplerian-like Problems
The original idea of our present work is based in [
23]. After choosing the free parameters
, we get pairs named NEW5(4) and form tables like
Table 4. There, we record the efficiency ratios of DP5(4) vs. NEW5(4) and calculate the average of means. This later value is a fitness measure and meant to be maximized. For the maximization process we applied Differential Evolution technique [
24]. We used MATLAB Software DeMat [
25] for implementing the latter technique.
The differential evolution process has already used by us for producing numerical methods for solving IVP problems and obtaining very interesting results. However, until now, we tried optimization in one or two problems. Notice that here we actually operate on runs for 98 problems! In [
26], we produced Nyström pairs of orders 4(3) and 6(4), while in [
27], we derived Numerov type algorithms of sixth orders for integration of orbits. Lately, we have tried this technique for Runge–Kutta pairs of orders 6(5) [
28,
29].
The optimization through differential evolution furnished various quintuplets for the parameters and we present the selected one in double precision below,
The resulting pair is presented in
Table 5. Coefficients not appearing there are zero.
We fixed
, since this coefficient actually affects only the tolerance. Indeed we may choose
and a new
Then, the tolerance
simply becomes
. A small value of
shifts the results towards stringent tolerances. We have to choose
(or
) with some care. Here we tried to get tables of the form given in
Table 4, as full as possible.
The norm of the principal truncation error coefficients is
, which is a little smaller than the corresponding value
of DP5(4). On the other hand, it is of the same magnitude with
observed for T5(4). The interval of absolute stability is
, and it is not extended. No extra phase-lag order is observed since
[
30].
In conclusion, no extra property seems to hold. The pair given in
Table 5 does not possess something interesting. It is hard to believe a special performance after seeing its traditional characteristics.
For the above selection of free parameters, we obtained the efficiency rations tabulated in
Table 6.
The average efficiency is , i.e., DP5(4) is about more expensive for these 98 runs over the selected 14 orbital problems. In reverse, this means that about digits were gained on average for NEW5(4) at the same costs.
6. Numerical Tests
We tested NEW5(4), presented here along with DP5(4) pair, given in [
11]. The latter pair is by far the most used Runge–Kutta pair. Everything else presented until now is hardly more efficient [
15,
19]. Both pairs were run for tolerances of
We set NEW5(4) as the reference pair and thus numbers greater than 1 indicate that NEW5(4) is more efficient. Thus, we can interpret the number
as DP5(4) being
more expensive than NEW5(4) while an entry of
means that DP5(4) is
more expensive (i.e., has twice the cost for achieving the same accuracy).
In the numerical tests, the errors are measured over the whole mesh on the integration interval. An estimation of the true solution is made every step by a parallel integration at stringent tolerance using an eighth order pair [
31]. Thus, the almost true global error is recorded. In the previous sections and for optimization, we preferred using the end point error since this reduced considerably the evaluation times. Thus, we repeat
Table 6 using global errors and report the results in
Table 7. The average is
justifying the results over the end-point only.
We proceed presenting in
Table 8 ratios for global errors for Kepler and perturbed Kepler (i.e., problems 1–10) in the interval
. We almost observe a slightly improvement in the results in comparison to the previous
Table 7. The overall average observed was
for problems 1–10.
We conclude by running Kepler and perturbed Kepler with different parameters. Thus, let us name as problems
Kepler orbits with eccentricities
and
, respectively. We also name as problems
the perturbed Kepler orbits with parameters
, respectively. The efficiency ratios for these 10 problems are presented in
Table 9.
The overall average observed was for problems –, which is also a quite positive result.
Finally, we mention that we manage to construct various pairs with coefficients near the one presented. By this, we mean that the distance was no more than, say . These pairs also gave good results. Perhaps 5–10% worse than the one presented. This means that no strict property holds for the pair, i.e., no equations like order conditions, conservation laws or some symplectic property holds. Our new proposed pair lies within an area of coefficients where it seems good for the kind of Keplerian-type obits we are interested here.
7. Conclusions
Training the coefficients of a Runge–Kutta pair for addressing a certain kind problems is considered. We concentrated on problems of Keplerian-type orbits and an extensively studied family of Runge–Kutta pairs of orders five and four. After optimizing the free parameters (coefficients) in an extensive set of problems and tolerances we concluded to a certain pair. This pair is found to outperform other representatives from this family in a wide range of relevant problems, intervals and tolerances.