Runge–Kutta–Nyström Pairs of Orders 8(6) for Use in Quadruple Precision Computations

Abstract

The second-order system of non-stiff Initial Value Problems (IVP) is considered and, in particular, the case where the first derivatives are absent. This kind of problem is interesting since since it arises in many significant problems, e.g., in Celestial mechanics. Runge–Kutta–Nyström (RKN) pairs are perhaps the most successful approaches for solving such type of IVPs. To achieve a pair attaining orders eight and six, we have to solve a well-defined set of equations with respect to the coefficients. Here, we provide a simplified form of these equations in a robust algorithm. When creating such pairings for use in double precision arithmetic, numerous conditions are often satisfied. First and foremost, we strive to keep the coefficients’ magnitudes small to prevent accuracy loss. We may, however, allow greater coefficients when working with quadruple precision. Then, we may build pairs of orders eight and six with significantly smaller truncation errors. In this paper, a novel pair is generated that, as predicted, outperforms state-of-the-art pairs of the same orders in a collection of important problems.

1. Introduction

Second order initial value problems (IVP) of the specific form

$${\zeta }^{″}=\psi (x,\zeta ),\zeta \left(x_0\right)={\zeta }_0,{\zeta }^{\prime }\left(x_0\right)={\zeta }_0^{\prime }$$

(1)

with $\psi :\mathbb{R}\times {\mathbb{R}}^m\to {\mathbb{R}}^m$ sufficiently continuous differentiable and $({\zeta }_0,{\zeta }_0^{\prime })\in {\mathbb{R}}^{2m},$ is considered here.

We evaluate an approximation to the solution of problem (1) at a set of distinct points $\left(x_n,{\zeta }_n,{\zeta }_n^{\prime }\right)$ using an explicit Runge–Kutta–Nyström method of algebraic order p. This method’s format is as follows (see [1] and ([2], p. 283) for further explanations on these methods):

$${\psi }_i=\psi (x_n+c_i{\tau }_n,{\zeta }_n+c_i{\tau }_n{\zeta }_n^{\prime }+{\tau }_n^2\sum _{j=1}^{i-1}{d}_{ij}{\psi }_j),i=1,2,\cdots ,s$$

$$\begin{array}{c}{\zeta }_{n+1}={\zeta }_n+{\tau }_n{\zeta }_n^{\prime }+{\tau }_n^2\sum _{i=1}^s{w}_i{\psi }_i,\hfill \\ {\zeta }_{n+1}^{\prime }={\zeta }_n^{\prime }+{\tau }_n\sum _{i=1}^s{w}_i^{\prime }{\psi }_i,\hfill \end{array}$$

where ${\tau }_n=x_{n+1}-x_n,$ is the size of the step. The last 50 years have seen a persistent interest in these methods. Follow the works of E. Fehlberg [3], Dormand et al. [4,5], El-Mikkawy and Rahmo [6], Papakostas et al. [7], Simos et al. [8], and Jerbi et al. [9]. There have also been RKN methods presented with special features. Houven et al. investigated RKN methods with lower phase lags, while Yoshida [10] and Calvo et al. [11] built RKN algorithms with the symplectic property.

In this paper, we set $p=8$ and we combine the aforementioned method with an additional formula of order six. In light of this, we also calculate an estimate of fifth order using the same values of ${\psi }_i$,

$$\begin{array}{c}{\widehat{\zeta }}_{n+1}={\zeta }_n+{\tau }_n{\zeta }_n^{\prime }+{\tau }_n^2\sum _{i=1}^s{\widehat{w}}_i{\psi }_i,\hfill \\ {\widehat{\zeta }}_{n+1}^{\prime }={\zeta }_n^{\prime }+{\tau }_n\sum _{i=1}^s{\widehat{w}}_i^{\prime }{\psi }_i.\hfill \end{array}$$

The higher order approximations ${\zeta }_n,{\zeta }_n^{\prime }$ are employed in all circumstances to propagate the solutions in time.

As a result, we have an estimator of error

$$ϵ=\parallel {\zeta }_{n+1}-{\widehat{\zeta }}_{n+1}\parallel =O\left({\tau }^7\right).$$

Then, we compare $ϵ$ with $TOL$ which is a small positive number defined by the user. Then, using this little value, known as tolerance, we can estimate the length of the following step as

$${\tau }_{n+1}=0.9·{\tau }_n·{\left(\frac{TOL}{ϵ}\right)}^{1/7},$$

(2)

which is in common use for the RKN8(6) pairs [4,12]. When $TOL<ϵ$, we prevent the solution from propagating. We essentially repeat the current step, but this time we use ${\tau }_{n+1}$ as the new shorter version rather than ${\tau }_n$.

The Butcher tableau may be used to represent the coefficients [13]. As a result, the method appears in the form

$$\begin{array}{cc}c\hfill & D\hfill \\ & \begin{array}{c}w,\widehat{w}\hfill \\ w^{\prime },{\widehat{w}}^{\prime }\hfill \end{array}\end{array}$$

with $D\in {\mathbb{R}}^{s\times s},$$w^{\mathrm{T}},{\widehat{w}}^{\mathrm{T}},{w^{\prime }}^{\mathrm{T}},{\widehat{w}}^{\mathrm{T}},c\in {\mathbb{R}}^s$, i.e. the weights are represented as row vectors.

Below we take into account a triplet with nine stages ($s=9$). The Butcher tableau shown in

Table 1. The Butcher tableau of RKN pairs of orders 8(6).

$0$
$c_2$	$d_{21}$
$c_3$	$d_{31}$	$d_{32}$
$c_4$	$d_{41}$	$d_{42}$	$d_{43}$
$c_5$	$d_{51}$	$d_{52}$	$d_{53}$	$d_{54}$
$c_6$	$d_{61}$	$d_{62}$	$d_{63}$	$d_{64}$	$d_{65}$
$c_7$	$d_{71}$	$d_{72}$	$d_{73}$	$d_{74}$	$d_{75}$	$d_{76}$
$1$	$d_{81}$	$d_{82}$	$d_{83}$	$d_{84}$	$d_{85}$	$d_{86}$	$d_{87}$
$1$	$w_1$	0	$w_3$	$w_4$	$w_5$	$w_6$	$w_7$	0
8th-order $w$	$w_1$	0	$w_3$	$w_4$	$w_5$	$w_6$	$w_7$	0	0
6th-order $\widehat{w}$	${\widehat{w}}_1$	0	${\widehat{w}}_3$	${\widehat{w}}_4$	${\widehat{w}}_5$	${\widehat{w}}_6$	${\widehat{w}}_7$	0	0
8th-order $w^{\prime }$	$w_1^{\prime }$	0	$w_3^{\prime }$	$w_4^{\prime }$	$w_5^{\prime }$	$w_6^{\prime }$	$w_7^{\prime }$	$w_8^{\prime }$	$w_9^{\prime }$
6th-order ${\widehat{w}}^{\prime }$	${\widehat{w}}_1^{\prime }$	0	${\widehat{w}}_3^{\prime }$	${\widehat{w}}_4^{\prime }$	${\widehat{w}}_5^{\prime }$	${\widehat{w}}_6^{\prime }$	${\widehat{w}}_7^{\prime }$	${\widehat{w}}_8^{\prime }$	${\widehat{w}}_9^{\prime }$

displays its coefficients.

Using such a method, we spend only eight stages per step since the final stage is reused as first stage in the following step. Thus, the numbers in the ninth stage coincide with w, i.e., $d_{9j}=w_j$ for $j=1,2,\cdots ,8$. This technique is known as FSAL (First Stage As Last).

RKN pairs of orders eight and six that use effectively eight stages per step were studied in [5,7]. Eighth order RKN methods that share seven stages per step have only been constructed for the special case of linear inhomogeneous problems [14].

2. Runge–Kutta–Nyström Methods of Eighth Order

We apply a RKN method to (1) and deploy the Taylor series expansions for $\zeta (x_n+\tau )-{\zeta }_{n+1}$ and ${\zeta }^{\prime }(x_n+\tau )-{\zeta }_{n+1}^{\prime }$. When matching expressions up to $h^8$ for an eighth order method, the following results are obtained:

$$\zeta (x_n+\tau )-{\zeta }_{n+1}={\tau }^2{e}_{2,1}{Q}_{2,1}+{\tau }^3{e}_{3,1}{Q}_{3,1}+\cdots +{\tau }^8\left(e_{8,1}{Q}_{8,1}+\dots +e_{8,20}{Q}_{8,20}\right)+O\left({\tau }^9\right)$$

(3)

$${\zeta }^{\prime }(x_n+\tau )-{\zeta }_{n+1}^{\prime }=\tau {\tilde{e}}_{1,1}{Q}_{1,1}+{\tau }^2{\tilde{e}}_{2,1}{Q}_{2,1}+\cdots +{\tau }^8\left({\tilde{e}}_{8,1}{Q}_{8,1}+\dots +{\tilde{e}}_{8,36}{Q}_{8,36}\right)+O\left({\tau }^9\right)$$

(4)

where $e_{ij}$ are expressions depending on $w,D,c$ while ${\tilde{e}}_{ij}$ are expressions depending on $w^{\prime },D,c$. An algorithm for their symbolic derivation of is given in [15]. Expressions $Q_{ij}$ are elementary differentials with respect to ${\zeta }^{\prime },\psi$ and its partial derivatives. The elementary differentials come from the problem and can not be handled by the method. However, for an eighth order RKN method we have to eliminate the coefficients $e_{ij}$ and ${\tilde{e}}_{ij}$ in expressions (3) and (4) to the requested order. In

Table 2. Number of equations of conditions for RKN methods.

		Order
$method$	$\underbrace{numberof}$ - order→	1	2	3	4	5	6	7	8	9	10
RKN	order conditions for $\zeta$	0	1	1	2	3	6	10	20	36	72
	order conditions for ${\zeta }^{\prime }$	1	1	2	3	6	10	20	36	72	137

, we list the number of order conditions (i.e., of $e_{ij}$ and ${\tilde{e}}_{ij}$) for each order. Thus, e.g., for a third algebraic order method we have to nullify $0+1+1=2$ equations for $\zeta$ and another $1+1+2=4$ order conditions for ${\zeta }^{\prime }$.

Inspecting from Butcher tableau above (i.e., Table 1) the number of coefficients available for a nine stages method and compared it with the equations of condition up to eighth order as reported in Table 2, we see that are far too less. Thus, we proceed making various simplifying assumptions that drastically reduce the number of order conditions. Firstly, we set

$$w=w^{\prime }·(I_s-C),$$

(5)

with $I_s\in {\mathbb{R}}^{s\times s}$ the identity matrix and $C=\mathrm{diag}\left(c\right)$. Using this assumption we automatically satisfy the order conditions for $\zeta$ after eliminating the equations of the same order for ${\zeta }^{\prime }$. Then, we are interested in eliminating only ${\tilde{e}}_{ij}$ with respect to $w^{\prime },D,c$.

Again summing the numbers in the last row of Table 2, we see that are remaining to too many for the coefficients at hand. Thus, we proceed making the following assumptions

$$D·\mathbb{I}=\frac{1}{2}{c}^2,D·c=\frac{1}{6}{c}^3,D·c^2=\frac{1}{12}{c}^4,$$

(6)

with

$$c^i=\underset{i-\mathrm{times}}{\underbrace{c\circ \cdots \circ c\circ c}},$$

the componentwise multiplication among matrices (i.e., Hadamard multiplication), while $c^0=\mathbb{I}={[1,1,\cdots ,1]}^T\in {\mathbb{R}}^s$. This multiplication has lower priority than dot product.

We also consider the row simplifying condition for RKN methods

$$w^{\prime }·(D+C-\frac{1}{2}(C\circ C)-\frac{1}{2}{I}_s)$$

and the subsidiary simplifying assumptions

$${(w·D)}_2=0,{(w^{\prime }·D)}_2=0,{(w^{\prime }·(C\circ C)·D)}_2=0,{(\widehat{w}·D)}_2=0.$$

Then we achieve a severe reduction in the number of order conditions and we may continue deriving the coefficients of an eighth order method (i.e., $w,w^{\prime },D$ and c) by the following algorithm.

• BEGIN ALGORITHM

Select arbitrary values for the coefficients ${\widehat{w}}_9^{\prime },w_9^{\prime },c_4,c_5,c_6$ and $c_7$.

Then compute successively and explicitly

$$c_3=\frac{\left\{\begin{array}{c}15-20c_4-20c_5+28c_4{c}_5-20c_6+28c_4{c}_6+28c_5{c}_6-42c_4{c}_5{c}_6\\ -20c_7+28c_4{c}_7+28c_5{c}_7-42c_4{c}_5{c}_7+28c_6{c}_7\\ -42c_4{c}_6{c}_7-42c_5{c}_6{c}_7+70c_4{c}_5{c}_6{c}_7\end{array}\right\}}{\left\{\begin{array}{c}2(10-14c_4-14c_5+21c_4{c}_5-14c_6+21c_4{c}_6+21c_5{c}_6\\ -35c_4{c}_5{c}_6-14c_7+21c_4{c}_7+21c_5{c}_7-35c_4{c}_5{c}_7\\ +21c_6{c}_7-35c_4{c}_6{c}_7-35c_5{c}_6{c}_7+70c_4{c}_5{c}_6{c}_7)\end{array}\right\}}$$

$$c_2=\frac{1}{2}{c}_3$$

Solve Vandermonde equations

$$w^{\prime }·e=1,w^{\prime }·c=\frac{1}{2},w^{\prime }·c^2=\frac{1}{3},w^{\prime }·c^3=\frac{1}{4},$$

$$w^{\prime }·c^4=\frac{1}{5},w^{\prime }·c^5=\frac{1}{6},w^{\prime }·c^6=\frac{1}{7},$$

for $w_1^{\prime },w_3^{\prime },w_4^{\prime },w_5^{\prime },w_6^{\prime },w_7^{\prime },w_8^{\prime }$. The last Vandermonde equation $w^{\prime }·c^7=\frac{1}{8}$ is automatically satisfied by the selection of $c_3$. Then vector w comes explicitly from (5).

Solve ${(D·c)}_4=\frac{c_4^2}{2},{(D·c^2)}_4=\frac{c_4^3}{6}$, for $d_{42}$ and $d_{43}$. Notice here that $(D·c)$ is a vector and, thus, ${(D·c)}_4$ is its fourth component.

Solve ${(w·D)}_2=0,{(w^{\prime }·D)}_2=0$, ${(w^{\prime }·(C\circ C)·D)}_2=0$ ans ${(\widehat{w}·D)}_2=0$ for $d_{72},d_{62},d_{52},d_{82}$.

Solve the following three integral equations

$$\begin{array}{ccc}& & w^{\prime }·(C-I_s)·(C-c_7{I}_s)·D·(C-c_3{I}_s)·(C-c_4{I}_s)·c\hfill \\ & =& {\int }_0^1(x-1)(x-c_7){\int }_0^x{\int }_0^x(x-x_3)(x-c-4)xdxdxdx,\hfill \end{array}$$

$$\begin{array}{ccc}& & w^{\prime }·(C-I_s)·D·(C-c_3{I}_s)·(C-c_4{I}_s)·(C-c_5{I}_s)·c\hfill \\ & =& {\int }_0^1(x-1){\int }_0^x{\int }_0^x(x-x_3)(x-c-4)(x-c_5)xdxdxdx,\hfill \end{array}$$

$$\begin{array}{ccc}& & w^{\prime }·(C-I_s)·D·(C-c_3{I}_s)·(C-c_4{I}_s)·(C-c_5{I}_s)·c\hfill \\ & =& {\int }_0^1(x-1){\int }_0^x{\int }_0^x(x-x_3)(x-c-4)(x-c_5)xdxdxdx,\hfill \end{array}$$

for $d_{65},d_{76},d_{75}$.

Evaluate $a_{53},a_{54},a_{63},a_{64},a_{73},a_{74}$ from ${(D·c)}_j=\frac{c_j^2}{2},{(D·c^2)}_j=\frac{c_j^3}{6}$, for $j=5,6,7$.

Evaluate $a_{83},a_{84},a_{85},a_{86},a_{87}$ from $w^{\prime }·(D+C-\frac{1}{2}(C\circ C)-\frac{1}{2}{I}_s)$ and its respective coordinates.

In the end we get

$$a_{9i}=w_i,\mathrm{for}i=1,2,\cdots ,8$$

and

$$a_{j1}=\frac{c_j^2}{2}-\sum _{k=2}^{j-1}{a}_{jk},j=2,3,\cdots ,8.$$

• END OF ALGORITHM

It is of interest to remark that no such simplified algorithm has ever appeared until now. It helped us very much in deriving our triplet.

3. Producing a RKN Pair of Orders 8 and 6

By the algorithm given in the previous section we may construct an eighth order RKN method at an actual cost of eight stages per step. This algorithm is that it offers six free parameters. Thus, we may exploit them in order to improve the efficiency of our new method. We choose to minimize the terms of the principal error terms, i.e., the Euclidean norm of the ninth order coefficients $e_{9j},j=1,2,\cdots ,36$ and ${\tilde{e}}_{9j},j=1,2,\cdots ,72$, that appear in series expansions (3) and (4).

When utilizing double precision arithmetic, we traditionally aim to keep the magnitude of the coefficients as small as possible. The margins of the available digits would therefore be severely tested by a coefficient of size ${10}^3$, a function value of size ${10}^2$ and a tolerance of $\epsilon ={10}^{-11}$. However, with quadruple precision, we are still able to accept these large coefficients with tolerances as low as approximately ${10}^{-25}$. By allowing the coefficients to increase, we can now move on to a new minimization procedure.

In order to address this, we choose to use Differential Evolution Algorithm [16,17]. Differential Evolution is an iterative procedure and in every iteration, named generation g, we work with a “population” of individuals ${\left({\widehat{w}}_9^{\left(g\right)},w_9^{\prime }{}^{\left(g\right)},c_4^{\left(g\right)},c_5^{\left(g\right)},c_6^{\left(g\right)},c_7^{\left(g\right)}\right)}_i$, $i=1,2,\cdots ,P$ with P the population size. An initial population ${\left({\widehat{w}}_9^{\left(0\right)},w_9^{\prime }{}^{\left(0\right)},\cdots ,c_7^{\left(0\right)}\right)}_i$, $i=1,2,\cdots ,P$ is randomly created in the first step of the method. We have also set as fitness function the expression

$$\mathsf{\Phi }=\sqrt{e_{9,1}^2+e_{9,2}^2+\cdots +e_{9,36}^2}+\sqrt{{\tilde{e}}_{9,1}^2+{\tilde{e}}_{9,2}^2+\cdots +{\tilde{e}}_{9,72}^2}={∥T^{\left(9\right)}∥}_2+{∥{T^{\prime }}^{\left(9\right)}∥}_2$$

which measures the loss from a ninth order method. The fitness function is then evaluated for each individual in the initial population and it is meant to be minimized. A three-phase sequential procedure updates every participant in each generation (iteration) g. The stages are Differentiation, Crossover, and Selection. We used the software DeMat [18] in MATLAB [19] where the latter technique is implemented. Success is not guaranteed with one optimization. Thus, we run the procedure several hundred times. Then we manage to get a solution. The result is further refined in order to get more digits of accuracy. This was done working on multi-precision arithmetic and using the function NMinimize of Mathematica [20].

The main characteristics of the major RKN pairs of eighth order studied here are given in

Table 3. Basic characteristics of the RKN pairs considered.

Pair	Stages	FSAL	${∥{\mathit{T}}^{\left(9\right)}∥}_2$	${∥{\mathit{T}}^{\prime }{}^{\left(9\right)}∥}_2$
PT$8\left(6\right)$ [7]	9	YES	$1.7·{10}^{-7}$	$1.6·{10}^{-7}$
DEP$8\left(6\right)$ [4]	9	YES	$8.3·{10}^{-7}$	$8.2·{10}^{-7}$
RKNT$8\left(6\right)$	9	YES	$1.7·{10}^{-8}$	$1.6·{10}^{-8}$

. The norms presented there correspond to the Euclidean norm of the coefficients ninth order (i.e., of ${\tau }^9$) in the expressions (3) and (4). We expect the new method to perform better since it produces significantly reduced local truncation errors.

The coefficients of the method constructed can be found in the following Mathematica module where we also implemented the integration algorithm we used in numerical tests in Listing 1.

Listing 1.Mathematica Module implementing the new pair.

RKNT86[f_List, vars_List, initialvalues_List, dinitialvalues_List, finalx_, errorTolerance_] := Block[{x = SetAccuracy[First[initialvalues], 33], y = SetAccuracy[Rest[initialvalues], 33], dy = SetAccuracy[dinitialvalues, 33], xend = SetAccuracy[finalx, 33], h, err, hnew, y7, y8, dy7, dy8, solution, hmax, ireject, k = Table[i + j, {i, 9}, {j, Length[vars] - 1}], b=SetAccuracy[{46704396222138759/1124501888012545693,0, 84069894477030747/424535379079037893,60269691739898297/328032958547368465, 2009963068113133/27794099874007722,162341471393132/140140455957185117, 6086576956589044/1882413506280312633,0,0},33], bb=SetAccuracy[{10769958754260247/261191895425614637,0, 104933541030533329/527807735255158343,8187542127950603/44863180380403502, 50493885750265423/674323734860213804,-396215365808089/252398506959352750, 5468871271464350/1319483122963052413,0,0},33], db=SetAccuracy[{46704396222138759/1124501888012545693,0, 90371972523959954/390135632629351589,118990880894033457/367654647557162744, 180830119624415039/624884373647391279,16628088200566168/1643600751401035359, 1524820183138666476/417332398303375801,-942444174868320016/265473221553563103, 0},33], dbb=SetAccuracy[{10769958754260247/261191895425614637,0, 58861559987617091/253105545276009947,142913350550568712/444546485690175277, 8398007711885933/28026591338889651,-8440103966850896/615634893567208211, 1592393294195924241/339999309740023022, -6699802037196600096/1421037300124099357,3/20},33], c=SetAccuracy[{0,8065253268/111157879849,16130506536/111157879849,99/229, 1855/2473,116/131,1129/1130,1,1},33], a=SetAccuracy[{{0},{502615833312847/190946037812928939}, {1601030787675953/456179150746555700,1601030787675953/228089575373277850}, {47478115875661981/518814108724307373,-64883723802385428/357040639400014459, 25666007926449694/139746227660637731}, {-328112826298039228/251912779790891183,969895830706346953/297412056373654755, -958305119264262743/492487831928632961,151603443293999467/564549369158251216}, {44079989458325648760/345626831710945999,-267609305840442666747/859338149021870938, 130442442641184422881/655209191357439877,-7381158156698807543/475346800759815547, 594932629852457670/835908452635682287}, {-10802627635977292643/544607328597417370,22047268993379696720/454307750813938153, -9705881798108421635/315306127829247354,1078781161885226048/413453123878982063, -8616008188673363/388077019471353686,365346507915481/466435620062528214}, {-13306779498890004275/660225117657805349,22208114914951831801/450387553598953907, -6398475501845852180/204556450443208783,1412284034546646006/533270054097053815, -19179472816466775/820785347597843378,14435103384615/18331075303513484, -364401779978/904202609357507829}, {46704396222138759/1124501888012545693,0,84069894477030747/424535379079037893, 60269691739898297/328032958547368465,2009963068113133/27794099874007722, 162341471393132/140140455957185117,6086576956589044/1882413506280312633,0}},33], half = SetAccuracy[1/2, 33], two = SetAccuracy[2, 33], ninetenth = SetAccuracy[9/10, 33], ode = Function[Release[vars], Release[f]]}, hmax = SetAccuracy[xend - x, 33]; h = SetAccuracy[errorTolerance^(1/8), 33]; ireject = 0; solution = SetAccuracy[{Flatten[{x, y, dy}]}, 33]; k[[9]] = SetAccuracy[Apply[ode, Flatten[{x, y}]], 33]; While[x < xend, If[x + h == x, Break[]]; If[x + h - xend > 0, h = SetAccuracy[xend - x, 33]]; k[[1]] = k[[9]]; k[[2]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[2]] h, y + c[[2]] h*dy + h^2*a[[2, 1]] k[[1]]}]], 33]; k[[3]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[3]] h, y + c[[3]] h*dy + h^2*a[[3]].Take[k, 2]}]], 33]; k[[4]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[4]] h, y + c[[4]] h*dy + h^2*a[[4]].Take[k, 3]}]], 33]; k[[5]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[5]] h, y + c[[5]] h*dy + h^2*a[[5]].Take[k, 4]}]], 33]; k[[6]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[6]] h, y + c[[6]] h*dy + h^2*a[[6]].Take[k, 5]}]], 33]; k[[7]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[7]] h, y + c[[7]] h*dy + h^2*a[[7]].Take[k, 6]}]], 33]; k[[8]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[8]] h, y + c[[8]] h*dy + h^2*a[[8]].Take[k, 7]}]], 33]; k[[9]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[9]] h, y + c[[9]] h*dy + h^2*a[[9]].Take[k, 8]}]], 33]; y8 = SetAccuracy[y + h*dy + h^2 b.k, 33]; y7 = SetAccuracy[y + h*dy + h^2 bb.k, 33]; dy8 = SetAccuracy[dy + h db.k, 33]; dy7 = SetAccuracy[dy + h dbb.k, 33]; err = SetAccuracy[Max[{Max[Abs[y8 - y7]/10], Max[Abs[dy8 - dy7]/10]}], 33]; hnew = SetAccuracy[ Min[SetAccuracy[hmax, 33],SetAccuracy[ h/Max[half, Min[two, SetAccuracy[(Rationalize[err/errorTolerance, 10^-33]^(1/7))/ninetenth, 33]]], 33]], 33]; If[err <= errorTolerance, x += h; y = y8; dy = dy8; AppendTo[solution, SetAccuracy[Flatten[{x, y8, dy8}], 33]], hnew = SetAccuracy[Min[hnew, h], 33]; ireject = ireject + 1; k[[9]] = k[[1]]];h = hnew]; Return[{solution,ireject}]]

In the listing above we give in the input

• f: A list with function $\psi$, i.e., in the form $\{{\psi }_1(x,{}^1\zeta ,{}^2\zeta ,\cdots ),{\psi }_2(x,{}^1\zeta ,{}^2\zeta ,\cdots ),\cdots \}$
• vars: A list with variables, i.e., in the form $\{x,{}^1\zeta ,{}^2\zeta ,\cdots \}$
• initial values: A list with initial values, i.e., in the form $\{x_0,{}^1\zeta \left(x_0\right),{}^2\zeta \left(x_0\right),\cdots \}$
• dinitialvalues: A list with initial values for ${\zeta }^{\prime }$, i.e., in the form ${\{}^1{\zeta }^{\prime }\left(x_0\right),{}^2{\zeta }^{\prime }\left(x_0\right),\cdots \}$
• finalx: The final value of x
• errorTolerance: The tolerance $TOL$

Notice here that by ${}^1\zeta ,{}^2\zeta ,\cdots$ we denote the coordinates of $\zeta$.

In the output we get the list solution where we collected all the steps taken, i.e., the values

$$x_0,{}^1{\zeta }_0,{}^2{\zeta }_0,\cdots ,{}^m{\zeta }_0,{}^1{\zeta }_0^{\prime },{}^2{\zeta }_0^{\prime },\cdots ,{}^m{\zeta }_0^{\prime }$$

$$x_1,{}^1{\zeta }_1,{}^2{\zeta }_1,\cdots ,{}^m{\zeta }_1{,}^1{\zeta }_1^{\prime },{}^2{\zeta }_1^{\prime },\cdots ,{}^m{\zeta }_1^{\prime }$$

$$x_2,{}^1{\zeta }_2,{}^2{\zeta }_2,\cdots ,{}^m{\zeta }_2,{}^1{\zeta }_2^{\prime },{}^2{\zeta }_2^{\prime },\cdots ,{}^m{\zeta }_2^{\prime }$$

$$\cdots$$

etc. We also get the number ireject of the rejected steps.

4. Numerical Results

In the following we may present some numerical tests to support the value of our new proposal.

4.1. The Methods

The explicit 8th order methods selected for testing are the following:

• The RKN pair of orders 8(6) given in [5], named DEP8(6).
• The RKN pair of orders 8(6) given in [7], named PT8(6).
• The RKN pair of orders 8(6) presented here, named RKNT8(6).

These pairs were run as normal, with an error estimation $ϵ$ being evaluated at each step. The Formula (2) is then used to create the new step since the error’s asymptotical form for them is $O\left({\tau }^7\right)$. The framework presented in the previous section was used for all runs.

4.2. The Problems

For our tests, a number of well-known problems were selected out of the literature. These problems were run for tolerances ${10}^{-14},{10}^{-15}\cdots ,{10}^{-22}$. For all these runs we recorded the steps taken (accepted and rejected) and the maximum global error observed at the end-point. The observed numbers (i.e., steps vs errors) are shown in various efficiency plots. All computations were done in Mathematica.

4.2.1. Inhomogeneous Equation

The inhomogeneous equation is the first test problem selected [21]:

$${\zeta }^{″}=-100\zeta \left(x\right)+99·sin\left(x\right),\zeta \left(0\right)=1,{\zeta }^{\prime }\left(0\right)=11,$$

with theoretical solution

$$\zeta \left(t\right)=cos\left(10x\right)+sin\left(10x\right)+sin\left(x\right).$$

We integrated this problem in the interval $x\in \left[0,10\pi \right]$. In Figure 1, we draw the corresponding efficiency plots.

4.2.2. Inhomogeneous Linear System

This problem is given as follows:

$${\zeta }^{″}=\left[\begin{array}{cc}\frac{1}{100}& -\frac{1}{10}\\ -\frac{1}{10}& \frac{1}{100}\end{array}\right]·\zeta +\left[\begin{array}{c}0\\ sinx\end{array}\right],$$

with theoretical solution

$$\zeta =\left[\begin{array}{c}cos\frac{3}{10}x-\frac{1000}{10,101}sinx\\ cos\frac{3}{10}x-\frac{10,100}{10,101}sinx\end{array}\right]$$

We integrated that problem in the interval $x\in \left[0,10\pi \right]$ and the efficiency plots are shown in Figure 2.

The rightmost circle in Figure 2 is justified by the following run of the Mathematica module given in the listing of the previous section.

In[1]:={solution, ireject} = RKNT86[{1/100*z1-1/10*z2,-1/10*z1+1/100*z2+Sin[x]}, {x,z1,z2}, {0,1,1},{-1000/10101,-(10100/10101)},10*Pi,10^-22]; In[2]={Length[solution] + ireject - 1, Max[Abs[Last[solution][[2 ;; 5]] - {-1, -1, -1000/10101, -10100/10101}]]} Out[2]={6957,2.419274*10^-26}

4.2.3. Problem F

Next we consider the following problem:

$${}^1{\zeta }^{″}=-4x^2·{}^1\zeta -2\frac{{}^2\zeta }{\sqrt{{}^1{\zeta }^2+{}^2{\zeta }^2}},{}^2{\zeta }^{″}=-4x^2·{}^2\zeta +2\frac{{}^1\zeta }{\sqrt{{}^1{\zeta }^2+{}^2{\zeta }^2}},x\in \left[\sqrt{\frac{1}{2}\pi },10\right]$$

with initial values

$${}^1\zeta \left(\sqrt{\frac{1}{2}\pi }\right)=0,{}^2\zeta \left(\sqrt{\frac{1}{2}\pi }\right)=1,{}^1{\zeta }^{\prime }\left(\sqrt{\frac{1}{2}\pi }\right)=-\left(\sqrt{2\pi }\right),{}^2{\zeta }^{\prime }\left(\sqrt{\frac{1}{2}\pi }\right)=0$$

and theoretical solution

$${}^1\zeta \left(x\right)=cos{x}^2,{}^2\zeta =sin{x}^2.$$

Here, ${}^1\zeta$ and ${}^2\zeta$ may understood as components and not as time steps.

We integrated that problem in the interval $x\in \left[0,10\pi \right]$. The theoretical solution and the efficiency plots are shown in Figure 3.

4.2.4. Kepler Problem

Next, we choose the celebrated two body problem with eccentricity $e=0.5$,

$$\begin{array}{c}\hfill {}^1{\zeta }^{″}=-\frac{{\zeta }_1}{{{(}^1{\zeta }^2+{}^2{\zeta }^2)}^{3/2}},{}^2{\zeta }^{″}=-\frac{{}^2\zeta }{{{(}^1{\zeta }^2+{}^2{\zeta }^2)}^{3/2}},\\ \hfill {}^1\zeta \left(0\right)=\frac{1}{2},{}^1{\zeta }^{\prime }\left(0\right)=0,{}^2\zeta \left(0\right)=0,{}^2{\zeta }^{\prime }\left(0\right)=\sqrt{3}.\end{array}$$

We solved the above equations in the interval $\left[0,10\pi \right]$ as $\zeta \left(10\pi \right)=\zeta \left(0\right)$ and the theoretical solution is given in [22]. The efficiency plots are shown in Figure 4.

4.2.5. Coupled Non-Linear Pendulum

Finally, we considered a smooth variant of the non-linear problem from [2] pg. 297. The equations of motion are

$$\begin{array}{ccc}\hfill {}^1{\zeta }^{″}& =& -sin{(}^1\zeta )-0.2(sin{(}^1\zeta )-sin{(}^2\zeta ))cos{(}^1\zeta )+e^{-10x},\hfill \\ \hfill {}^2{\zeta }^{″}& =& -sin{(}^2\zeta )-0.1(sin{(}^2\zeta )-sin{(}^1\zeta ))cos{(}^2\zeta ).\hfill \end{array}$$

We integrated the problem in the interval $x\in \left[0,496\right]$ starting with no move or speed at all, i.e., ${}^1\zeta \left(0\right)={}^2\zeta \left(0\right)={}^1{\zeta }^{\prime }\left(0\right)={}^2{\zeta }^{\prime }\left(0\right)=0$.

No theoretical solution is available. Thus, the almost true solution is derived by running an integration using $TOL={10}^{-26}$. The efficiency plots are shown in Figure 5.

4.3. Discussion on the Results

The results show that the new pair outperforms by far other RKN8(6) pairs on the problems tested. Almost 1–2 digits of accuracy were gained in most cases. The findings demonstrate that when high accuracies are required in the solution of special second order IVPs, the new approach significantly beats earlier ones.

5. Conclusions

In this paper, we looked at Runge–Kutta–Nyström pairs that are specifically optimized for solving initial value problems of second order when the first derivative is missing. We exploited the large coefficients that can be handled when working in quadruple precision arithmetic. The main contribution of our try is that the proposed method possesses very smaller truncation error terms in comparison with the other eighth order pairs appeared in the literature. Numerical tests in relevant problems justify our effort.

Reconsideration of other high order Runge–Kutta and Runge–Kutta–Nyström pairs may considered in the future in the direction shown here. Minimization of truncation errors for use in quadruple precision, regardless the increase in the magnitude of the coefficients.