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

Abstract

We consider a family of explicit Runge–Kutta pairs of orders six and five without any additional property (reduced truncation errors, Hamiltonian preservation, symplecticness, etc.). This family offers five 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 efficient 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.

1. Introduction

The initial value problem (IVP) is

$$y^{\prime }=f(x,y),y\left(x_0\right)=y_0$$

(1)

with $x_0\in \mathbb{R}$, $y,y^{\prime }\in {\mathbb{R}}^m$ and $f:\mathbb{R}\times {\mathbb{R}}^m\to {\mathbb{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,{\widehat{b}}^T,c\in R^s$ and $A\in R^{s\times 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\left(x_{n+1}\right).$ Namely, $y_{n+1}$ and ${\widehat{y}}_{n+1}$, given by

$$y_{n+1}=y_n+h_n\sum _{i=1}^s{b}_i{f}_i$$

and

$$y_{n+1}=y_n+h_n\sum _{i=1}^s{\widehat{b}}_i{f}_i$$

with

$$f_i=f(x_n+c_i{h}_n,y_n+h_n\sum _{j=1}^{i-1}{a}_{ij}{f}_j),$$

for $i=1,2,\cdots ,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

$${ϵ}_n=h_n^{p-q-1}·\parallel y_{n+1}-{\widehat{y}}_{n+1}\parallel$$

is formed in every step and is combined in an algorithm for changing the step size:

$$h_{n+1}=0.9·h_n^{}·{\left(\frac{t}{{ϵ}_n}\right)}^{1/p},$$

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.

2. 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$ and ${\widehat{b}}_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

$$\begin{array}{c}{}^1{y}^{\prime }={}^{3}y,\hfill \\ {}^2{y}^{\prime }={}^{4}y,\hfill \\ {}^3{y}^{\prime }=-{\displaystyle \frac{{}^{1}y}{{\left(\sqrt{{}^1{y}^2+{}^2{y}^2}\right)}^3}},\hfill \\ {}^4{y}^{\prime }=-{\displaystyle \frac{{}^{2}y}{{\left(\sqrt{{}^1{y}^2+{}^2{y}^2}\right)}^3}},\hfill \end{array}$$

with $x\in [0,10\pi ]$, $y\left(0\right)={\left[1-ecc,0,0,\sqrt{\frac{1+ecc}{1-ecc}}\right]}^T$ and the theoretical solution

$${}^{1}y\left(x\right)=cos\left(v\right)-ecc,{}^{2}y\left(x\right)=sin\left(v\right)\sqrt{1-ec{c}^2}.$$

In the above, $v=ecc·sin\left(u\right)+x$, $ecc$ is the eccentricity, and the components of y are denoted by the left superscript. They should not be confused with $y_1={\left[{}^1{y}_1{,}^2{y}_1{,}^3{y}_1{,}^4{y}_1\right]}^T,y_2={\left[{}^1{y}_2{,}^2{y}_2{,}^3{y}_2{,}^4{y}_2\right]}^T,y_3,\cdots$, which represent the vectors approximating the solution at $x_1,x_2,x_3,\cdots$.

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 $fev$ 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=fev·g{e}^{1/6}.$$

(2)

Running DLMP6(5) twice for Kepler we obtained the efficiency measures ${\widehat{u}}_1$ and ${\widehat{u}}_2$ as reported in

Table 1. Efficiency measures for both runs and pairs.

	$\mathit{ecc}$	${\mathit{x}}_{\mathit{end}}$	$\mathit{tol}$	DLMP6(5)	NEW6(5)	Ratio $\widehat{\mathit{u}}/\tilde{\mathit{u}}$
First run	0	$10\pi$	${10}^{-7}$	${\widehat{u}}_1=127.22$	${\tilde{u}}_1=50.64$	$2.51$
Second run	$0.6$	$20\pi$	${10}^{-11}$	${\widehat{u}}_2=833.27$	${\tilde{u}}_2=386.64$	$2.16$

. 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\pi$ and $tol={10}^{-7}$. Thus, we obtained ${\widehat{u}}_1=1123·2.14·{10}^{-6}\approx 127.22$, as reported in Table 1. Analogously, for a second run with $ecc=0.6$, $x_{end}=20\pi$ and $tol={10}^{-11}$ we observed ${\widehat{u}}_2=833.27$.

Let us suppose that any new pair NEW6(5) furnishes corresponding efficiency measures ${\tilde{u}}_1$ and ${\tilde{u}}_2$ for the same runs. Then, as a fitness function we form the sum

$$u=\frac{{\widehat{u}}_1}{{\tilde{u}}_1}+\frac{{\widehat{u}}_2}{{\tilde{u}}_2},$$

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$ and ${\widehat{b}}_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. Coefficients of the proposed NEW6(5) pair, accurate for double precision computations.

$c_2=0.173146279530013,$	$c_3=0.163620769891761,$	$c_4=0.245431154837642,$
$c_5=0.452502877641229,$	$c_6=0.902924768667267,$	$c_7=0.8101151362080617,$
$c_1=0,c_8=c_9=1,$	$b_1=0.0794169052387116,$	$b_2=b_3=0,$
$b_4=0.320063598496390,$	$b_5=0.179217292937057,$	$b_6=-0.2872484367615202,$
$b_7=0.573172758378662,$	$b_8=0.135377881710699,$	$b_9=0$
${\widehat{b}}_1=0.148854176113754,$	${\widehat{b}}_2={\widehat{b}}_3=0,$	${\widehat{b}}_4=0.291009331941132,$
${\widehat{b}}_5=0.229278395578701,$	${\widehat{b}}_6=-0.1155397766857130,$	${\widehat{b}}_7=0.429687174664803,$
${\widehat{b}}_8=0.0167106983873234,$	${\widehat{b}}_9=0.064345053530889,$
$a_{21}=0.173146279530013,$	$a_{31}=0.0863111204651556,$	$a_{32}=0.077309649426606,$
$a_{41}=0.061357788709411,$	$a_{42}=0,$	$a_{43}=0.184073366128232,$
$a_{51}=0.178735636864969,$	$a_{52}=0,$	$a_{53}=-0.430121641642955,$
$a_{54}=0.703888882419215,$	$a_{61}=-0.3492563988707026,$	$a_{62}=0,$
$a_{63}=4.2286674995349015,$	$a_{64}=-5.131590895887595,$	$a_{65}=2.155104563890663,$
$a_{71}=-0.004184382566843,$	$a_{72}=0,$	$a_{73}=1.062724280290705,$
$a_{74}=-1.188530484293243,$	$a_{75}=0.8944565948851806,$	$a_{76}=0.045649127892262,$
$a_{81}=-0.518393300452978,$	$a_{82}=0,$	$a_{83}=4.607278279969559,$
$a_{84}=-5.004120306973807,$	$a_{85}=1.510536380616834,$	$a_{86}=-0.399249451366671,$
$a_{87}=0.803948398207063,$	$a_{9j}=b_j,j=1,2,\cdots ,8.$

.

Interpreting Table 1, we observe that DLMP6(5) was $151\%$ and $116\%$ more expensive than NEW6(5) for the first and second run, respectively. The norm of the principal truncation error coefficients was $\parallel T^{\left(7\right)}{\parallel }_2\approx 2.64·{10}^{-4}$, which is much larger than the corresponding value $\parallel T^{\left(7\right)}{\parallel }_2\approx 4.37·{10}^{-5}$ for DLMP6(5). The interval of absolute stability for the new pair was $(-4.24,0]$ while for for DLMP6(5) it was $(-4.21,0]$.

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.

3. Numerical Tests

We tested the following pairs chosen from the family studied above.

• The DLMP6(5) pair, given in [15].
• NEW6(5), presented here.

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},\cdots ,{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.

1. 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\pi$ and $x_{end}=20\pi$.

The efficiency measure ratios of DLMP6(5) vs. NEW6(5) for Kepler are presented in

Table 3. Efficiency measure ratios of DLMP6(5) vs. NEW6(5) for Kepler in $[0,10\pi ]$.

$\frac{\mathbf{tol}->}{\underbrace{\mathbf{ecc}}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
0	$1.90$	$2.36$	$3.66$	$3.06$	$2.46$	$3.18$	$3.59$	$2.89$
$0.2$	$1.15$	$1.11$	$1.14$	$1.01$	$1.07$	$1.46$	$2.07$	$1.29$
$0.4$	$1.09$	$1.10$	$1.11$	$1.51$	$0.80$	$1.27$	$1.71$	$1.23$
$0.6$	$1.49$	$1.49$	$1.86$	$1.27$	$1.18$	$1.43$	$1.73$	$1.49$
$0.8$	$1.13$	$1.17$	$1.32$	$1.21$	$1.36$	$1.08$	$2.73$	$1.43$

and

Table 4. Efficiency measure ratios of DLMP6(5) vs. NEW6(5) for Kepler in $[0,20\pi ]$.

$\frac{\mathbf{tol}->}{\underbrace{\mathbf{ecc}}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
0	$1.87$	$2.34$	$2.59$	$2.35$	$2.67$	$3.36$	$3.36$	$2.65$
$0.2$	$1.22$	$1.27$	$1.33$	$1.24$	$1.29$	$1.67$	$2.39$	$1.49$
$0.4$	$1.31$	$1.55$	$1.21$	$0.96$	$1.08$	$2.36$	$2.27$	$1.54$
$0.6$	$1.38$	$1.32$	$1.17$	$1.10$	$0.98$	$1.14$	$1.99$	$1.30$
$0.8$	$1.16$	$1.24$	$1.85$	$1.25$	$1.08$	$0.94$	$1.61$	$1.30$

.

2. The perturbed Kepler   

This problem describes the motion of a planet according to Einstein’s general relativity theory and the Schwarzschild potential is applied. The equations are:

$$\begin{array}{c}{}^1{y}^{″}=-{\displaystyle \frac{{}^{1}y}{{\sqrt{{}^1{y}^2+{}^2{y}^2}}^3}}-(2+\delta )\delta {\displaystyle \frac{{}^{1}y}{{\sqrt{{}^1{y}^2+{}^2{y}^2}}^5}},\hfill \\ {}^2{y}^{″}=-{\displaystyle \frac{{}^{2}y}{{\sqrt{{}^1{y}^2+{}^2{y}^2}}^3}}-(2+\delta )\delta {\displaystyle \frac{{}^{2}y}{{\sqrt{{}^1{y}^2+{}^2{y}^2}}^5}},\hfill \end{array}$$

and the analytical solution is

$$\begin{array}{c}{}^{1}y=cos(x+\delta x),{}^{2}y=sin(x+\delta x).\hfill \end{array}$$

We transformed this problem into a system of four first-order equations and solved for $x_{end}=10\pi$ and $x_{end}=20\pi$. After recording the endpoint errors and the costs, we present the efficiency measures ratios of DLMP6(5) vs. NEW6(5) for perturbed Kepler in

Table 5. Efficiency measures ratios of DLMP6(5) vs. NEW6(5) for perturbed Kepler in $[0,10\pi ]$.

$\frac{\mathbf{tol}->}{\underbrace{\mathit{\delta }}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
$0.01$	$1.90$	$2.48$	$2.39$	$2.33$	$3.04$	$3.64$	$3.49$	$2.75$
$0.02$	$1.93$	$2.69$	$2.19$	$2.11$	$2.58$	$3.65$	$3.45$	$2.66$
$0.03$	$1.91$	$2.76$	$2.14$	$2.07$	$2.50$	$3.45$	$3.47$	$2.61$
$0.04$	$1.87$	$2.63$	$2.18$	$2.12$	$2.53$	$3.28$	$3.57$	$2.60$
$0.05$	$1.87$	$2.40$	$2.32$	$2.30$	$2.47$	$3.24$	$3.72$	$2.62$

and

Table 6. Efficiency measures ratios of DLMP6(5) vs. NEW6(5) for perturbed Kepler in $[0,20\pi ]$.

$\frac{\mathbf{tol}->}{\underbrace{\mathit{\delta }}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
$0.01$	$1.88$	$2.44$	$2.35$	$2.18$	$2.43$	$3.11$	$3.34$	$2.53$
$0.02$	$1.86$	$2.49$	$2.29$	$2.14$	$2.38$	$3.09$	$3.36$	$2.52$
$0.03$	$1.85$	$2.44$	$2.33$	$2.17$	$2.43$	$3.17$	$3.39$	$2.54$
$0.04$	$1.83$	$2.38$	$2.36$	$2.17$	$2.44$	$3.21$	$3.43$	$2.55$
$0.05$	$1.82$	$2.32$	$2.37$	$2.16$	$2.43$	$3.19$	$3.44$	$2.53$

.

3. The Arenstorf orbit  

Another interesting orbit describes the stable movement of a spacecraft around Earth and Moon ([22], pg. 129).

$$\begin{array}{ccc}\hfill {}^1{y}^{″}& =& {}^{1}y+2·{{}^{2}y}^{\prime }-{\zeta }^{\prime }·\frac{{}^{1}y+\zeta }{P_1}-\zeta ·\frac{{}^{1}y-{\zeta }^{\prime }}{P_2},\hfill \\ \hfill {}^2{y}^{″}& =& {}^{2}y+2·{{}^{1}y}^{\prime }-{\zeta }^{\prime }·\frac{{}^{2}y}{P_1}-\zeta ·\frac{{}^{2}y}{P_2},\hfill \end{array}$$

with

$$\begin{array}{ccc}\hfill P_1& =& {\sqrt{{\left({}^{1}y+\zeta \right)}^2+{{}^{2}y}^2}}^3,P_2={\sqrt{{\left({}^{1}y-{\zeta }^{\prime }\right)}^2+{{}^{2}y}^2}}^3,\hfill \\ \hfill \zeta & =& 0.012277471,{\zeta }^{\prime }=0.987722529,\hfill \end{array}$$

initial values

$${}^{1}y\left(0\right)=0.994,{{}^{1}y}^{\prime }\left(0\right)=0,{}^{2}y\left(0\right)=0,{{}^{2}y}^{\prime }\left(0\right)=-2.00158510637908252,$$

and with $x_A=17.0652165601579625589$ the solution is periodic.

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. Efficiency measure ratios of DLMP6(5) vs. NEW6(5) for Arenstorf.

$\frac{\mathbf{tol}->}{\underbrace{{\mathit{x}}_{\mathbf{end}}}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
$x_A$	$1.13$	$1.42$	$1.89$	$1.34$	$0.90$	$1.51$	$1.71$	$1.41$
$2x_A$	$0.96$	$1.82$	$1.64$	$1.07$	$1.09$	$1.62$	$1.58$	$1.40$

.

4. The Pleiades  

Finally, we considered the problem “Pleiades” as given in ([22], pg. 245).

$${}^i{y}^{″}=\sum _{i\ne j}\frac{{\mu }_j\left({}^{j}y-{}^{i}y\right)}{{\rho }_{ij}},{{}^{i}z}^{″}=\sum _{i\ne j}\frac{{\mu }_j\left({}^{j}z{-}^{i}z\right)}{{\rho }_{ij}},$$

with

$${\rho }_{ij}={\sqrt{{\left({}^{i}y-{}^{j}y\right)}^2+{\left({}^{i}z-{}^{j}z\right)}^2}}^3,i,j=1,\cdots ,7.$$

The initial values are

$${}^{1}y\left(0\right)=3,{}^{2}y\left(0\right)=3,{}^{3}y\left(0\right)=-1,{}^{4}y\left(0\right)=-3,{}^{5}y\left(0\right)=2,{}^{6}y\left(0\right)=-2,{}^{7}y\left(0\right)=2,$$

$${}^{1}z\left(0\right)=3,{}^{2}z\left(0\right)=-3,{}^{3}z\left(0\right)=2,{}^{4}z\left(0\right)=0,{}^{5}z\left(0\right)=0,{}^{6}z\left(0\right)=-4,{}^{7}z\left(0\right)=4,$$

$${}^1{y}^{\prime }\left(0\right)=0,{}^2{y}^{\prime }\left(0\right)=0,{}^3{y}^{\prime }\left(0\right)=0,{}^4{y}^{\prime }\left(0\right)=0,{}^5{y}^{\prime }\left(0\right)=0,{}^6{y}^{\prime }\left(0\right)=1.75,{}^7{y}^{\prime }\left(0\right)=-1.5,$$

$${}^1{z}^{\prime }\left(0\right)=0,{}^2{z}^{\prime }\left(0\right)=0,{}^3{z}^{\prime }\left(0\right)=0,{}^4{z}^{\prime }\left(0\right)=-1.25,{}^5{z}^{\prime }\left(0\right)=1,{}^6{z}^{\prime }\left(0\right)=0,{}^7{z}^{\prime }\left(0\right)=0,$$

We again transformed this problem to a system of fourteen first-order equations and solved it to $x_{end}=3$ and 4. We recorded the endpoint errors after we estimated the solution there by a very accurate integration using Mathematica and quadruple precision. The efficiency measure ratios of DLMP6(5) vs. NEW6(5) for Pleiades can be found in

Table 8. Efficiency measures ratios of DLMP6(5) vs. NEW6(5) for Pleiades.

$\frac{\mathbf{tol}->}{\underbrace{{\mathit{x}}_{\mathbf{end}}}}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$	${10}^{-11}$	Mean
3	$1.12$	$1.02$	$0.97$	$1.03$	$0.89$	$1.10$	$1.63$	$1.11$
4	$1.13$	$1.03$	$0.97$	$0.96$	$0.94$	$1.15$	$1.64$	$1.12$

.

We estimated 168 (i.e., 12 problems times 7 tolerances times two end points) efficiency measures for each pair. In average we observed a ratio of $1.98$, meaning that DLMP6(5) is about $98\%$ more expensive! This is quite remarkable since a great deal of effort has been put over the years towards achieving 10–20% efficiency [23,24]. In reverse, this means that about ${log}_{10}1.{98}^6\approx 1.8$ digits were gained on average at the same costs.

4. 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.