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

: The second-order system of non-stiff Initial Value Problems (IVP) is considered and, in particular, the case where the ﬁrst derivatives are absent. This kind of problem is interesting since since it arises in many signiﬁcant 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-deﬁned set of equations with respect to the coefﬁcients. Here, we provide a simpliﬁed form of these equations in a robust algorithm. When creating such pairings for use in double precision arithmetic, numerous conditions are often satisﬁed. First and foremost, we strive to keep the coefﬁcients’ magnitudes small to prevent accuracy loss. We may, however, allow greater coefﬁcients when working with quadruple precision. Then, we may build pairs of orders eight and six with signiﬁcantly 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.


Introduction
Second order initial value problems (IVP) of the specific form with ψ : R × R m → R m sufficiently continuous differentiable and (ζ 0 , ζ 0 ) ∈ R 2m , is conside • red here.We evaluate an approximation to the solution of problem (1) at a set of distinct points (x n , ζ n , ζ n ) 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): where τ 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 wo • rks 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 ψ i , The higher order approximations ζ n , ζ n are employed in all circumstances to propagate the solutions in time.
As a result, we have an estimator of error 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 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 τ n+1 as the new shorter version rather than τ n .The Butcher tableau may be used to represent the coefficients [13].As a result, the method appears in the form c D w, ŵ w , ŵ with D ∈ R s×s , w T , ŵT , w T , ŵT , c ∈ 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 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, • • • , 8.This technique is known as FSAL (First Stage As Last).
RKN p• airs of orders eight and six that use effectively eight st• ages 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].

Runge-Kutta-Nyström Methods of Eighth Order
We apply a RKN method to (1) and deploy the Taylor series expansions for ζ(x n + τ) − ζ n+1 and ζ (x n + τ) − ζ n+1 .When matching expressions up to h 8 for an eighth order method, the following results are obtained: (3) whe • re e ij are exp • ressions depending on w, D, c while ẽij are exp • ressions depending on w , D, c.An algorithm for their symbolic derivation of is given in [15].Expressions Q ij are elementary differentials with respect to ζ , ψ 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 ẽij in expressions ( 3) and ( 4) to the requested order.In Table 2, we list the number of order conditions (i.e., of e ij and ẽij ) for each order.Thus, e.g., for a third algebraic order method we have to nullify 0 + 1 + 1 = 2 equations for ζ and another 1 + 1 + 2 = 4 order conditions for ζ .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 with I s ∈ R s×s the identity matrix and C = diag(c).Using this assumption we automatically satisfy the order conditions for ζ after eliminating the equations of the same order for ζ .Then, we are interested in eliminating only ẽij with respect to w , 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 with the componentwise multiplication among matrices (i.e., Hadamard multiplication), while This multiplication has lower priority than dot product.
We also consider the row simplifying condition for RKN methods and the subsidiary simplifying assumptions 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 , D and c) by the following algorithm.

BEGIN ALGORITHM
Select arbitrary values for the coefficients ŵ 9 , w 9 , c 4 , c 5 , c 6 and c 7 .
Then compute successively and explicitly for w 1 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 .The last Vandermonde equation w • c 7 = 1 8 is automatically satisfied by the selection of c 3 .Then vector w comes explicitly from (5).Solve (w Solve the following three integral equations ) and its respective coordinates.
In the end we get and

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.

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, • • • , 36 and ẽ9j , j = 1, 2, • • • , 72, that appear in series expansions ( 3) and ( 4).
When utilizing double precision arithmetic, we traditionally aim to keep the m• agnitude 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 ε = 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 ŵ(g) 9 , w 9 (g) , c , i = 1, 2, • • • , P is randomly created in the first step of the method.We have also set as fitness function the expression 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 threephase 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.The norms presented there correspond to the Euclidean norm of the coefficients ninth order (i.e., of τ 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.
We also get the number ireject of the rejected steps.

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

The Methods
The explicit 8th order methods selected for testing are the following:
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(τ 7 ).The framework presented in the previous section was used for all runs.

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 • • • , 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.

Inhomogeneous Equation
The inhomogeneous equation is the first test problem selected [21]: We integrated this problem in the interval x ∈ [0, 10π].In Figure 1, we draw the corresponding efficiency plots.

Inhomogeneous Linear System
This problem is given as follows:   The rightmost circle in Figure 2 is justified by the following run of the Mathematica module given in the listing of the previous section.

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.

c 3 4 6
, for d 42 and d 43 .Notice here that (D • c) is a vector and, thus, (D • c) 4 is its fourth component.


We integrated that problem in the interval x ∈ [0, 10π] and the efficiency plots are shown in Figure2.

Figure 1 .
Figure 1.Efficiency plots for the inhomogeneous equation

Figure 4 .
Figure 4. Efficiency plots for the Kepler problem.

Figure 5 .
Figure 5. Efficiency plots for the coupled non-linear pendulum.

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

Table 3 .
Basic characteristics of the RKN pairs considered.

Listing 1 .
Mathematica Module implementing the new pair.