Next Article in Journal
Dynamics of a Predator–Prey Model with the Additive Predation in Prey
Next Article in Special Issue
A Simple Affine-Invariant Spline Interpolation over Triangular Meshes
Previous Article in Journal
Prospective Primary Teachers’ Didactic-Mathematical Knowledge in a Service-Learning Project for Inclusion
Previous Article in Special Issue
Parallel Direct and Iterative Methods for Solving the Time-Fractional Diffusion Equation on Multicore Processors
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Runge-Kutta-Nyström Pairs of Orders 8(6) with Coefficients Trained to Perform Best on Classical Orbits

1
Department of Industrial Engineering, College of Engineering, University of Ha’il, Hail 1234, Saudi Arabia
2
Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah 21589, Saudi Arabia
3
Department of Electrical Engineering, College of Engineering, University of Ha’il, Hail 1234, Saudi Arabia
4
College of Applied Mathematics, Chengdu University of Information Technology, Chengdu 610225, China
5
Scientific and Educational Center “Digital Industry”, South Ural State University, 76 Lenin Ave., 454 080 Chelyabinsk, Russia
6
Department of Medical Research, China Medical University Hospital, China Medical University, Taichung City 40402, Taiwan
7
Data Recovery Key Laboratory of Sichun Province, Neijing Normal University, Neijiang 641100, China
8
Section of Mathematics, Department of Civil Engineering, Democritus University of Thrace, 67100 Xanthi, Greece
9
General Department, Euripus Campus, National & Kapodistrian University of Athens, GR34-400 Psachna, Greece
*
Author to whom correspondence should be addressed.
Submission received: 6 November 2021 / Revised: 4 February 2022 / Accepted: 13 February 2022 / Published: 19 February 2022
(This article belongs to the Special Issue Numerical Analysis and Scientific Computing II)

Abstract

:
In this study, we consider eight stages per step family of explicit Runge-Kutta-Nyström pairs of orders eight and six. The pairs from this family effectively use eight stages for each step. The coefficients provided by such a method are much less than the number of non linear order conditions required to be solved. Thus, we traditionally apply various simplified assumptions in order to address this drawback. The assumptions taken in the family we consider here deliver a subsystem where all the coefficients are evaluated successively and explicitly with respect to five free parameters. We train (adjust) these free parameters in order to derive a certain pair that outperforms other similar pairs of orders 8 ( 6 ) in Keplerian type orbits, e.g., Kepler, perturbed Kepler, Arenstorf orbit, or Pleiades. Differential evolution technique is used for the training. The pair that we finally present offers about an additional digit of accuracy in a variety of orbits.

1. Introduction

The initial value problem of special form
ψ = φ ( x , ψ ) , ψ ( x 0 ) = ψ 0 , ψ ( x 0 ) = ψ 0
where φ : R × R m R m , and ( ψ 0 , ψ 0 ) R 2 m , is considered here.
We approximate the solution of (1) at a set of discrete points x n , ψ n , ψ n with an explicit Runge-Kutta-Nyström (RKN) pair of orders p ( q ) , p > q . This method has the form
φ i = φ ( x n + γ i h n , ψ n + γ i h n ψ n + h n 2 j = 1 s α i j φ j ) , i = 1 , 2 , , s
ψ n + 1 = ψ n + h n ψ n + h n 2 i = 1 s β i φ i , ψ ^ n + 1 = ψ n + h n ψ n + h n 2 i = 1 s β ^ i φ i , ψ n + 1 = ψ n + h n i = 1 s β i φ i , ψ ^ n + 1 = ψ n + h n i = 1 s β ^ i φ i ,
with h n = x n + 1 x n , while the number of the stages of the method is s. The solution is propagated by the higher order approximations ψ n , ψ n . In consequence, β ^ and β ^ correspond to the lower order method.
Some norm of the vector estimating the error
ϵ n = max ( h n 2 i = 1 s ( β i β ^ i ) φ i , h n i = 1 s ( β i β ^ i ) φ i ) ,
is used in comparison with the requested tolerance T O L , for step-size control algorithm
h n + 1 = 0.9 · h n · T O L h n p q + 1 · ϵ n 1 / p .
This formula is used even if T O L < h p q + 1 · ϵ n , but then h n + 1 is actually a new smaller version of h n . For extended details in the issue, see [1].
All the coefficients can be formulated using the Butcher tableau [2,3]. So the method takes the form
γ A β β ^ β β ^
with A R s × s , β T , β ^ T , β T , β ^ T , γ R s . Here, matrix A is strictly lower triangular since the methods considered are explicit. By this we mean that the function evaluations (i.e., ϕ i ’s) are calculated directly and successively. On the contrary, when the method is implicit we have to solve non linear equations for evaluating the stages, which increases the computation time. Implicit methods are used when the problems are stiff [4]. The problems we are dealing with in the following do not fall into the latter category [5].
Here we are interested in studying nine stages (i.e., s = 9 ), FSAL (First stage As Last) pair of orders eight and six (i.e., p = 8 and q = 6 ). This method spends only eight stages per step since the last stage is used again as the first stage in the next step. Thus, the coefficients in the ninth stage coincide with β , i.e., a 9 j = β j for j = 1 , 2 , , 8 .

2. The Dormand-El Mikkawy-Prince Family of Runge-Kutta-Nyström Pairs of Orders 8(6)

This family is used for the derivation of the famous DEP8(6) pair presented in [6]. Further investigation is given in [1]. Thus, we may deploy here the algorithm for the derivation of the parameters of the pair with respect to the free terms γ 4 , γ 5 , γ 6 , γ 7 and β ^ 9 .
Initially, set γ 8 = γ 9 = 1 . All the coefficients not presented below are zero, e.g., β 2 = 0 .
γ 3 = 15 20 γ 4 20 γ 5 + 28 γ 4 γ 5 20 γ 6 + 28 γ 4 γ 6 + 28 γ 5 γ 6 42 γ 4 γ 5 γ 6 20 γ 7 + 28 γ 4 γ 7 + 28 γ 5 γ 7 42 γ 4 γ 5 γ 7 + 28 γ 6 γ 7 42 γ 4 γ 6 γ 7 42 γ 5 γ 6 γ 7 + 70 γ 4 γ 5 γ 6 γ 7 2 ( 10 14 γ 4 14 γ 5 + 21 γ 4 γ 5 14 γ 6 + 21 γ 4 γ 6 + 21 γ 5 γ 6 35 γ 4 γ 5 γ 6 14 γ 7 + 21 γ 4 γ 7 + 21 γ 5 γ 7 35 γ 4 γ 5 γ 7 + 21 γ 6 γ 7 35 γ 4 γ 6 γ 7 35 γ 5 γ 6 γ 7 + 70 γ 4 γ 5 γ 6 γ 7 ) ,
γ 2 = 1 2 γ 3 .
Then solve
b · e = 1 , b · γ = 1 2 , b · γ 2 = 1 3 , b · γ 3 = 1 4 , b · γ 4 = 1 5 , b · γ 5 = 1 6 , b · γ 6 = 1 7
for β 1 , β 3 , β 4 , β 5 , β 6 , β 7 and β 8 .
In these equations, “∗” is to be understood as a component-wise multiplication among vectors and has the lowest priority after all other operations. Also,
γ 2 = γ γ , γ 3 = γ 2 γ , etc .
This latter operation (“raising” a vector to a power) has highest priority and is evaluated before the dot products and “∗”.
Set β i = β i 1 γ i , for i = 1 , 2 , , 8 .
In consequence from FSAL property follows that a 9 i = β i , i = 1 , 2 , , 8 .
Continue evaluating successively
a 32 = γ 3 2 6 γ 2 ,
a 42 = γ 4 3 ( 2 γ 3 + γ 4 ) 12 γ 2 ( γ 2 γ 3 ) ,
a 43 = γ 4 3 ( 2 γ 2 + γ 4 ) 12 γ 3 ( γ 2 + γ 3 ) ,
a 65 = 9 20 γ 3 20 γ 4 + 56 γ 3 γ 4 12 γ 7 + 28 γ 3 γ 7 + 28 γ 4 γ 7 84 γ 3 γ 4 γ 7 10 , 080 ( γ 3 γ 5 ) γ 5 ( γ 4 + γ 5 ) ( 1 + γ 6 ) ( γ 6 γ 7 ) d b 6 ,
a 75 = 9 γ 5 + 20 γ 3 γ 5 + 20 γ 4 γ 5 56 γ 3 γ 4 γ 5 + 15 γ 6 32 γ 3 γ 6 32 γ 4 γ 6 + 84 γ 3 γ 4 γ 6 12 γ 6 2 + 28 γ 3 γ 6 2 + 28 γ 4 γ 6 2 84 γ 3 γ 4 γ 6 2 6 γ 7 + 12 γ 3 γ 7 + 12 γ 4 γ 7 28 γ 3 γ 4 γ 7 + 12 γ 5 γ 7 28 γ 3 γ 5 γ 7 28 γ 4 γ 5 γ 7 + 84 γ 3 γ 4 γ 5 γ 7 10 , 080 ( γ 3 γ 5 ) γ 5 ( γ 4 + γ 5 ) ( γ 5 γ 6 ) ( γ 6 γ 7 ) ( 1 + γ 7 ) β 7 ,
a 76 = 3 6 γ 3 6 γ 4 + 14 γ 3 γ 4 6 γ 5 + 14 γ 3 γ 5 + 14 γ 4 γ 5 42 γ 3 γ 4 γ 5 5040 ( γ 3 γ 6 ) γ 6 ( γ 4 + γ 6 ) ( γ 5 + γ 6 ) ( 1 + γ 7 ) β 7 ) ,
a 85 = β 5 + 2 γ 5 β 5 γ 5 2 β 5 + 2 a 65 β 6 + 2 a 75 β 7 2 β 8 ,
a 86 = β 6 + 2 γ 6 β 6 γ 6 2 β 6 + 2 a 76 β 7 2 β 8 ,
a 87 = β 7 2 γ 7 β 7 + γ 7 2 β 7 2 β 8 .
Then solve simultaneously
β ^ · e = 1 , β ^ · γ = 1 2 , β ^ · γ 2 = 1 3 , β ^ · γ 3 = 1 4 , β ^ · γ 4 = 1 5 , β ^ · γ 5 = 1 6
and
β ^ · A · Γ γ 3 I · Γ γ 4 I · γ γ 2 γ 2 γ 3 γ 2 γ 4 β ^ · A 2 = = x = 0 x = 1 x = 0 x = x x = 0 x = x x γ 3 x γ 4 x d x d x d x
for β ^ 1 , β ^ 3 , β ^ 4 , β ^ 5 , β ^ 6 , b ^ 7 and β ^ 8 , assuming that Γ = diag ( γ ) .
Also set β ^ i = β ^ i 1 γ i , for i = 1 , 2 , , 8 .
Proceed solving
b · A 2 = b · γ A 2 = b · γ 2 A 2 = β ^ · A 2 = 0
for a 52 , a 62 , a 72 and a 82 .
After this, evaluate successively
a 53 = 12 a 52 γ 2 2 + 12 a 52 γ 2 γ 4 2 γ 4 γ 5 3 + γ 5 4 12 γ 3 ( γ 3 γ 4 ) , a 54 = 12 a 52 γ 2 2 + 12 a 52 γ 2 γ 3 2 γ 3 γ 5 3 + γ 5 4 12 γ 4 ( γ 3 + γ 4 ) ) , a 63 = 12 a 62 γ 2 2 + 12 a 62 γ 2 γ 4 + 12 a 65 γ 4 γ 5 12 a 65 γ 5 2 2 γ 4 γ 6 3 + γ 6 4 12 γ 3 ( γ 3 γ 4 ) , a 64 = 12 a 62 γ 2 2 + 12 a 62 γ 2 γ 3 + 12 a 65 γ 3 γ 5 12 a 65 γ 5 2 2 γ 3 γ 6 3 + γ 6 4 12 γ 4 ( γ 3 + γ 4 ) , a 73 = 12 a 72 γ 2 2 + 12 a 72 γ 2 γ 4 + 12 a 75 γ 4 γ 5 12 a 75 γ 5 2 + 12 a 76 γ 4 γ 6 12 a 76 γ 6 2 2 γ 4 γ 7 3 + γ 7 4 12 γ 3 ( γ 3 γ 4 ) , a 74 = 12 a 72 γ 2 2 + 12 a 72 γ 2 γ 3 + 12 a 75 γ 3 γ 5 12 a 75 γ 5 2 + 12 a 76 γ 3 γ 6 12 a 76 γ 6 2 2 γ 3 γ 7 3 + γ 7 4 12 γ 4 ( γ 3 + γ 4 ) , a 83 = β 3 + 2 γ 3 β 3 γ 3 2 β 3 + 2 a 43 β 4 + 2 a 53 β 5 + 2 a 63 β 6 + 2 a 73 β 7 2 β 8 , a 84 = β 4 + 2 γ 4 β 4 γ 4 2 β 4 + 2 a 54 β 5 + 2 a 64 β 6 + 2 a 74 β 7 2 β 8 .
Conclude with
a i 1 = 1 2 γ i 2 j = 2 i 1 a i , j .
Because the ninth stage is used as the first stage of the next step, even though s = 9 , the family wastes only eight stages each 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 coefficients of h 8 in the residual of Taylor error expansions corresponding to the 5th order formula of the underlying RK pair.
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 Outcomes from Various Pairs

Comparison of the results observed by various pairs is an interesting issue. It is of crucial importance in our present work. We thus deal with the Kepler problem. This has a form that follows
1 ψ = 3 ψ , 2 ψ = 4 ψ , 3 ψ = 1 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 3 , 4 ψ = 2 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 3 ,
with x [ 0 , 10 π ] , ψ ( 0 ) = 1 τ , 0 , 0 , 1 + τ 1 τ T and the theoretical solution is [5]
1 ψ ( x ) = cos ( v ) τ , 2 ψ ( x ) = sin ( v ) 1 τ 2 .
In the above, v = τ · sin ( v ) + x , τ is the eccentricity, and the left superscript denotes the components of ψ . We have to avoid confusion with ψ 1 = 1 ψ 1 , 2 ψ 1 , 3 ψ 1 , 4 ψ 1 T , ψ 2 = 1 ψ 2 , 2 ψ 2 , 3 ψ 2 , 4 ψ 2 T , ψ 3 , , that correspond to the vectors approximating the solution at x 1 , x 2 , x 3 , .
When running DEP8(6) for τ = 0.8 and tolerances 10 5 , 10 6 , , 10 11 , we get the results summarized in Table 1. There we recorded the error observed at the end-point and the stages used by each pair.
Then, we try the PT8(6) pair given in [1] and get the results presented in Table 2.
The question we raise now is regarding which pair delivers the better results. It is not clear from what can be seen in the tables. Thus, we use a technique similar to the one described in [5].
We form a linear least squares fit on the logarithms of data (i.e., function evaluations and accuracies achieved). Then we arrive at a slope that corresponds to the order of the method. These parallel slopes for DEP8(6) and TP8(6) can be seen in Figure 1. See [3] (p. 214) for more explanations on this issue.
Now, we are able to compare the results of the pairs in the selected problem.
In Table 3 we present the function evaluations that might have been used for achieving various accuracies by both pairs. These values follow the linear fit, as mentioned in [5], but the constant ratio is explained by parallel slopes amongst techniques of the same order. The latter observation might not be present for every problem. Some special characteristics of various problems may be exploited differently by various pairs of the same orders and thus may provide steeper slopes than anticipated. As a result, the numbers in the second column of Table 3 are interpreted as
function evaluations 10 0.0879 · log 10 ( expected error ) + 2.742 ,
while the stages in the third column can be computed from the relation
function evaluations 10 0.0900 · log 10 ( expected error ) + 2.715 .
Averaging the ratios in the fourth column of Table 3, we conclude that DEP8(6)) is considered to be about 3 % costlier than PT8(6) in the problem under consideration. After this we may ask for free parameters that may derive a new pair, which is cheaper in a wider class of Keplerian type orbits.

4. Pairs’ Performance in a Diverse Range of Problems in Keplerian-like Forms

Our purpose here is to construct a RKN8(6) pair from the family discussed above. This pair has to deliver the best results on problems of Keplerian form. For this reason, we have selected the problems that follow.

4.1. The Kepler Problem

We have already dealt with this. We used the eccentricities τ = 0 , 0.2 , 0.4 , 0.6 , 0.8 , and recorded the cost (i.e., the stages used) and the errors observed at x e n d = 10 π . In the tables with the results, these five problems are referred to as the numbers one through five.

4.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:
1 ψ = 1 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 3 ( 2 + δ ) δ 1 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 5 , 2 ψ = 2 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 3 ( 2 + δ ) δ 2 ψ ( 1 ψ ) 2 + ( 2 ψ ) 2 5 ,
and the analytical solution is
1 ψ = cos ( x + δ x ) , 2 ψ = sin ( x + δ x ) .
We solved the problem through x e n d = [ 0 , 10 π 1 + δ ] . The parameters used are δ = 0.01 , 0.02 , 0.03 , 0.04 , 0.05 . In the tables with the results, these 5 problems are referred to as numbers 6 through to 10, respectively.

4.3. The Arenstorf Orbit

Another fascinating orbit displays a spacecraft’s steady journey around the Earth and Moon [7] (p. 129).
1 ψ = 1 ψ + 2 · 2 ψ ζ · 1 ψ + ζ P 1 ζ · 1 ψ ζ P 2 , 2 ψ = 2 ψ + 2 · 1 ψ ζ · 2 ψ P 1 ζ · 2 ψ P 2 ,
with
P 1 = 1 ψ + ζ 2 + 2 ψ 2 3 , P 2 = 1 ψ ζ 2 + 2 ψ 2 3 , ζ = 0.012277471 , ζ = 0.987722529 ,
and initial values
1 ψ ( 0 ) = 0.994 , 1 ψ ( 0 ) = 0 , 2 ψ ( 0 ) = 0 , 2 ψ ( 0 ) = 2.00158510637908252 ,
and the solution is periodic with period x A = 17.0652165601579625589 . We solved the Arenstorf problem and recorded the cost (i.e., the stages used) and the errors observed to x A and 2 x A . In the tables with the results, these two problems are referred to as the numbers 11 and 12, respectively.

4.4. The Pleiades

Finally, we considered the problem “Pleiades” as given in [7] (p. 245).
i ψ = i j μ j j ψ i ψ ρ i j , i z = i j μ j j z i z ρ i j ,
with
ρ i j = i ψ j ψ 2 + i z j z 2 3 , i , j = 1 , , 7 and μ j = j , j = 1 , , 7 .
The initial values can be retrieved from [7]. We integrated the problem to x e n d = 3 and x e n d = 4 . Similarly, we recorded the cost (i.e., the stages used) and the errors observed at the endpoints. In the tables with the results, these two problems are referred to as the numbers 13 and 14, respectively.
We have formed a set of 14 problems and we tested the pairs DEP8(6) and PT8(6) over 7 tolerances, namely 10 5 , 10 6 , , 10 11 .
In Table 4 we record the corresponding ratios (as those appeared in the rightmost column of Table 3) of DEP8(6) vs. PT8(6). In the first column of Table 4 we see the expected errors. In the last row we record the average performances over all runs (i.e., tolerances used) for each problem. Numbers greater than 1 are in favor of the second pair. The overall average of these means is 0.99, which means that PT8(6) was in average about 1 % more expensive than DEP8(6) in all these 14 × 7 = 98 runs.
Our aim now is to choose the five free parameters in a way to get a new pair that will be more efficient and perform better than the 0.99 achieved above.

5. Training the Free Parameters in a Wide Set of Keplerian-like Problems

The idea we apply here is based on a previous work [8]. After selecting values for the free parameters γ 4 , γ 5 , γ 6 , γ 7 , β ^ 9 we get pairs named NEW8(6) and construct tables similar to Table 4. There we record the efficiency ratios of DEP8(6) vs. NEW8(6) and calculate the average of means. This latter value is a fitness measure and is meant to be maximized. For the maximization process we may choose from a variety of population-dependent techniques.
We have chosen the differential evolution (DE) technique [9]. We already have achieved some very interesting results using DE [10]. The latter work is about Runge-Kutta pairs of orders 5(4) addressing equations of the form ψ = φ ( x , ψ ) , i.e., a clearly different type of methods and problems. DE is an iterative procedure, and in every iteration, named generation g, we work with a “population” of individuals
γ 4 i ( g ) , γ 5 i ( g ) , γ 6 i ( g ) , γ 7 i ( g ) , β ^ 9 i ( g ) , i = 1 , 2 , , N ,
with N the population size. An initial population
γ 4 i ( 0 ) , γ 5 i ( 0 ) , γ 6 i ( 0 ) , γ 7 i ( 0 ) , β ^ 9 i ( 0 ) , i = 1 , 2 , , N ,
is randomly created in the first step of the method. We also used the average of the means from Table 4 as the fitness function. Following that, the fitness function is assessed for each member in the initial population. A three-phase sequential approach updates all of the persons participating in each iteration (generation) g. Differentiation, crossover, and selection are the phases involved. For the latter method, we utilized MATLAB Software DeMat [11].
Among the numerous quintuplets we got by the technique selected, we concluded the following
γ 4 = 0.4556145825203227 , γ 5 = 0.494497106631637 , γ 6 = 0.8105140017857914 ,
γ 7 = 0.898444913211217 and β ^ 9 = 0.02601695275050284 .
In Table 5 we present the derived pair. Coefficients that did not appear there are zero.
For the above selection of the free parameters, we get the efficiency rations tabulated in Table 6.
The average efficiency is 1.29, i.e., DEP8(6) is about 29 % more expensive for these 98 runs over the selected 14 orbital problems. In reverse, this means that about
log 10 1 . 29 8 0.88
digits were gained in average for NEW8(6) at the same costs.

6. Numerical Tests

We tested NEW8(6) presented here in comparison with the DEP8(6) pair, which appeared in [6]. DEP8(6) is the most widely used Runge–Kutta-Nystöm pair of orders 8(6) [12]. Both pairs were run for tolerances of 10 5 , 10 6 , , 10 11 . We set NEW8(6) as the reference pair and thus numbers greater than 1 indicate that NEW8(6) is more efficient. Thus, we can interpret the number 1.30 as DEP8(6) being 0.3 = 30 % more expensive than NEW8(6).
The main difference here is that we recorded the global errors observed over all grid points of the integrations. We proceed by presenting in Table 7 the ratios for global errors for Kepler and perturbed Kepler (i.e., problems 1–10) in the intervals [ 0 , 100 ] changing the parameters. Thus let us name as problems 1 b , 2 b , 3 b , 4 b , 5 b Kepler orbits with various random eccentricities, namely τ = 0.01 , 0.24 , 0.37 , 0.48 and 0.77 , respectively. We also name as problems 6 b , 7 b , 8 b , 9 b , 10 b the perturbed Kepler orbits with parameters δ = 0.015 , 0.025 , 0.035 , 0.045 and 0.055 , respectively.
We also recorded the global errors for Arenstorf for intervals through x e n d = 30 and x e n d = 40 , now naming the problems 11 b and 12 b , respectively. Finally we run Pleiades to x e n d = 1.75 and x e n d = 3 now naming the problems 13 b and 14 b , respectively. The efficiency ratios for these 14 problems are presented in Table 7. The latter four problems come with no analytical solutions. Thus the true values were approximated by a parallel integration with much stringent tolerance.
The overall image is almost the same as the one observed in the previous Table, which is also quite a positive result.
The implementation of Runge-Kutta-Nyström methods is a subject of ongoing interest. This year we presented the Runge-Kutta-Nyström pair RKNST7(5) of orders 7(5) [13]. Another issue of research is producing methods that are symplectic. The latter methods have the property of linear error growth with respect to the length of integration interval when applied to Kepler problems [14]. Their disadvantage is having higher truncation errors and that they need more stages per step. On the contrary, conventional methods (like the one presented here) have quadratic error growth. One of the most known eighth order symplectic methods is the 15-stages Y8, given by Yoshida [15] in Table 1, column D.
We run the methods NEW8(6), DEP8(6), RKNST7(5) and Y8 on problems 1b and 6b. The intervals of integration were lengthened to 10,000 π and 10,000 π / 1.015 , respectively, in favor of symplectic methods. We recorded the stages taken versus the end-point error achieved and plot it in Figure 2 and Figure 3.
RKNTS7(5) pair was less efficient at higher tolerances than eighth order pairs since it attains lesser algebraic accuracy and it was not designed to perform best on Keplerian type orbits. DEP8(6) remains less efficient than NEW8(6) on these intervals and problems. Symplectic methods are not competitive at all. It seems to share interesting properties for geometers but not for numerical analysts. Their poor performance has already been remarked elsewhere [16].

7. Conclusions

In this study, we considered the Runge-Kutta-Nyström pairs of orders 8(6) for addressing the special second order Initial Value Problem. We focused on problems of the Keplerian type. Thus, we proposed a pair with coefficients specially trained in order to address such kinds of orbits. Extensive results justified our effort.

Author Contributions

All authors have contributed equally. All authors have read and agreed to the published version of the manuscript.

Funding

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under grant No. (FP-058-43). The authors therefore gratefully acknowledge DSR for the technical and financial support.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Papakostas, S.N.; Tsitouras, C. High phase-lag order Runge-Kutta and Nyström pairs. SIAM J. Sci. Comput. 1999, 21, 747–763. [Google Scholar] [CrossRef]
  2. Butcher, J.C. On Runge-Kutta processes of high order. J. Aust. Math. Soc. 1964, 4, 179–194. [Google Scholar] [CrossRef] [Green Version]
  3. Butcher, J.C. Numerical Methods for Ordinary Differential Equations, 3rd ed.; John Wiley & Sons: Chichester, UK, 2016. [Google Scholar]
  4. Hairer, E.; Wanner, G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed.; Springer: Berlin, Germany, 1996. [Google Scholar]
  5. Enright, W.H.; Pryce, J.D. Two FORTRAN packages for assessing initial value methods. ACM Trans. Math. Softw. 1987, 13, 1–27. [Google Scholar] [CrossRef]
  6. Dormand, J.R.; El-Mikkawy, M.E.A.; Prince, P.J. High-Order Embedded Runge-Kutta-Nyström Formulae. IMA J. Numer. Anal. 1987, 7, 423–430. [Google Scholar] [CrossRef]
  7. Hairer, E.; Nørsett, S.P.; Wanner, G. Solving Ordinary Differential Equations I: Nonstiff Problems; Springer: Berlin, Germany, 1993. [Google Scholar]
  8. Tsitouras, C. Neural Networks With Multidimensional Transfer Functions. IEEE Trans. Neural Netw. 2002, 13, 222–228. [Google Scholar] [CrossRef] [PubMed]
  9. Storn, R.; Price, K. Differential evolution—A simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Opt. 1997, 11, 341–359. [Google Scholar] [CrossRef]
  10. Kovalnogov, V.N.; Fedorov, R.V.; Karpukhina, T.V.; Simos, T.E.; Tsitouras, C. Runge–Kutta Pairs of Orders 5(4) Trained to Best Address Keplerian Type Orbits. Mathematics 2021, 9, 2400. [Google Scholar] [CrossRef]
  11. DeMat. Available online: https://www.swmath.org/software/24853 (accessed on 23 August 2021).
  12. Brankin, R.W.; Gladwell, I.; Dormand, J.R.; Prince, P.J.; Seward, W.L. Algorithm 670 A Runge-kutta-Nyström Code. ACM Trans. Math. Softw. 1989, 15, 31–40. [Google Scholar] [CrossRef]
  13. Simos, T.E.; Tsitouras, C. On high order Runge–Kutta– Nyström pairs. J. Comput. Appl. Math. 2022, 400, 113753. [Google Scholar] [CrossRef]
  14. Calvo, M.P.; Sanz-Serna, J.M. High order symplectic Runge-Kutta-Nyström methods. SIAM J. Sci. Comput. 1993, 14, 1237–1252. [Google Scholar] [CrossRef]
  15. Yoshida, H. Construction of higher order symplectic integrators. Phys. Lett. A 1990, 150, 262–268. [Google Scholar] [CrossRef]
  16. Tsitouras, C. A tenth order symplectic Runge-Kutta and Nyström method. Celest. Mech. Dyn. Astron. 1999, 74, 223–230. [Google Scholar] [CrossRef]
Figure 1. Results for DEP86 and PT86, as given in Table 1 and Table 2 and their corresponding log-linear fits.
Figure 1. Results for DEP86 and PT86, as given in Table 1 and Table 2 and their corresponding log-linear fits.
Mathematics 10 00654 g001
Figure 2. Efficiency curves for problem 1b in the interval [0, 10,000 π ].
Figure 2. Efficiency curves for problem 1b in the interval [0, 10,000 π ].
Mathematics 10 00654 g002
Figure 3. Efficiency curves for problem 6b in the interval [0, 10,000 π /1.015].
Figure 3. Efficiency curves for problem 6b in the interval [0, 10,000 π /1.015].
Mathematics 10 00654 g003
Table 1. Outcomes of DEP8(6) on Kepler orbit for τ = 0.8 .
Table 1. Outcomes of DEP8(6) on Kepler orbit for τ = 0.8 .
ToleranceStagesEnd-Point Error
10 5 1089 6.4 × 10 4
10 6 1377 2.7 × 10 5
10 7 1769 2.6 × 10 7
10 8 2265 1.3 × 10 8
10 9 2889 6.9 × 10 8
10 10 3497 4.0 × 10 9
10 11 3785 2.5 × 10 10
Table 2. Outcomes of PT8(6) on Kepler orbit for τ = 0.8 .
Table 2. Outcomes of PT8(6) on Kepler orbit for τ = 0.8 .
ToleranceStagesEnd-Point Error
10 5 1161 5.0 × 10 4
10 6 1457 1.6 × 10 6
10 7 1833 4.5 × 10 7
10 8 2361 3.4 × 10 8
10 9 3057 3.7 × 10 9
10 10 3729 1.2 × 10 9
10 11 3769 2.4 × 10 10
Table 3. Stages that might have been used for achieving various accuracies by both pairs for Keplerian orbit with eccentricity τ = 0.8 .
Table 3. Stages that might have been used for achieving various accuracies by both pairs for Keplerian orbit with eccentricity τ = 0.8 .
Expected End-Point ErrorDEPP 8 ( 6 ) PT 8 ( 6 ) Efficiency Ratio
10 3 1013.97 965.53 1.05
10 4 1241.42 1187.91 1.05
10 5 1519.88 1461.51 1.04
10 6 1860.81 1798.13 1.03
10 7 2278.24 2212.28 1.03
10 8 2789.24 2721.81 1.02
10 9 3414.91 3348.71 1.02
10 10 4180.91 4119.99 1.01
Table 4. Efficiency ratios of DEP8(6) vs. PT8(6) for problems 1–14 and the corresponding end-point errors achieved.
Table 4. Efficiency ratios of DEP8(6) vs. PT8(6) for problems 1–14 and the corresponding end-point errors achieved.
Error 1234567891011121314
10 2 *********** 0.95 **
10 3 **** 1.05 ****** 0.98 1.00 1.02
10 4 **** 1.05 ***** 0.94 1.01 1.01 1.04
10 5 0.84 1.07 0.95 0.90 1.04 0.84 0.84 0.84 0.84 0.83 0.99 1.04 1.03 1.06
10 6 0.85 1.08 0.98 0.97 1.03 0.86 0.85 0.85 0.85 0.85 1.03 1.08 1.05 1.07
10 7 0.86 1.09 1.02 1.05 1.03 0.87 0.87 0.87 0.87 0.87 1.08 1.11 1.07 1.09
10 8 0.88 1.09 1.05 1.13 1.02 0.88 0.88 0.88 0.88 0.89 1.14 * 1.09 1.11
10 9 0.89 1.10 1.09 1.23 1.02 0.89 0.90 0.90 0.90 0.91 1.19 * 1.11 1.13
10 10 0.90 1.11 1.12 1.33 1.01 0.90 0.91 0.92 0.91 0.93 ** 1.13 1.15
10 11 0.92 1.11 1.16 ** 0.91 0.93 0.93 0.93 0.95 ****
10 12 0.93 1.12 *** 0.92 0.94 0.95 0.95 0.97 ****
10 13 0.94 **** 0.93 0.96 0.97 0.96 0.99 ****
mean-> 0.89 1.10 1.05 1.10 1.03 0.89 0.90 0.90 0.90 0.91 1.06 1.03 1.06 1.08
Table 5. Coefficients of the proposed here NEW8(6) pair, accurate for double precision computations.
Table 5. Coefficients of the proposed here NEW8(6) pair, accurate for double precision computations.
γ 2 = 0.0854544187688376031 , γ 3 = 0.170908837537675206 , γ 4 = 0.455614582520322714 , γ 5 = 0.494497106631637020 , γ 6 = 0.810514001785791327 , γ 7 = 0.898444913211216931 , γ 8 = γ 9 = 1 , β 1 = 0.0495023778457969496 , β 3 = 0.223315864614348454 , β 4 = 5.864310848696467 · 10 4 , β 5 = 0.176658022702874654 , β 6 = 0.0453762194992222526 , β 7 = 0.00456108425288804292 , β ^ 1 = 0.0493217331530729867 , β ^ 3 = 0.224007190882142852 , β ^ 4 = 0.00580373475137855214 , β ^ 5 = 0.183035611932723099 , β ^ 6 = 0.0443854481831987883 , β ^ 7 = 0.00505375060024082628 , β 1 = β 1 , β ^ 1 = β ^ 1 , β 3 = 0.269350192988574135 , β 4 = 0.00107723510961154486 , β 5 = 0.349469854713854025 , β 6 = 0.239470039616994250 , β 7 = 0.0449124154890862874 , β 8 = 0.0462178842360828093 , β ^ 3 = 0.270184029240960690 , β ^ 4 = 0.0106610768125419417 , β ^ 5 = 0.362086180581648925 , β ^ 6 = 0.234241308600661186 , β ^ 7 = 0.0497636382385428827 , β ^ 8 = 0.0190472342471524293 , β ^ 9 = 0.0260169527505028420 , a 32 = 0.00973661024949315254 , a 42 = 0.122821108259130461 , a 43 = 0.153641587946575897 , a 52 = 0.0264148295270339516 , a 53 = 0.103470702345032179 , a 54 = 0.0103732869329210154 , a 62 = 0.0839513409881428112 , a 63 = 0.142671597223573008 , a 64 = 0.164005790762850565 , a 65 = 0.266751419874429655 , a 72 = 0.273030769247765195 , a 73 = 0.160122716797143754 , a 74 = 1.25849331157904383 , a 75 = 1.02650962278825033 , a 76 = 0.0629905335176362299 , a 82 = 0.0238094759938050803 , a 83 = 0.322215841053004229 , a 84 = 0.448160499830497980 , a 85 = 0.581476734552232745 , a 86 = 0.0318063480094925576 , a 87 = 0.00501106135437686956 , a 9 j = β j , j = 1 , 2 , , 8 , a i 1 = 1 2 γ i 2 j = 2 i 1 a i , j .
Table 6. Efficiency ratios of DEP8(6) vs. NEW8(6) for various problems and end-point errors achieved.
Table 6. Efficiency ratios of DEP8(6) vs. NEW8(6) for various problems and end-point errors achieved.
Error 1234567891011121314
10 1 *********** 0.95 **
10 1 *********** 1.00 * 0.85
10 1 **** 1.03 ***** 0.84 1.05 0.89 0.89
10 1 * 1.15 0.94 * 1.06 ***** 0.91 1.11 0.92 0.93
10 1 1.51 1.14 0.98 0.86 1.09 1.52 1.52 1.51 1.54 1.49 0.98 1.17 0.95 0.97
10 1 1.52 1.14 1.02 0.94 1.11 1.53 1.54 1.52 1.53 1.51 1.06 1.24 0.98 1.02
10 1 1.53 1.13 1.06 1.03 1.14 1.55 1.55 1.54 1.53 1.52 1.15 1.31 1.02 1.07
10 1 1.54 1.13 1.10 1.13 1.17 1.57 1.57 1.55 1.53 1.53 1.24 * 1.06 1.12
10 1 1.55 1.12 1.14 1.24 1.20 1.58 1.58 1.56 1.52 1.54 1.34 * 1.09 1.17
10 1 1.56 1.12 1.19 1.36 * 1.60 1.59 1.58 1.52 1.55 ** 1.13 1.23
10 1 1.57 1.11 1.23 ** 1.62 1.61 1.59 1.52 1.56 ****
10 1 1.58 **** 1.64 1.62 1.60 1.51 1.58 ****
10 1 1.59 **** 1.65 1.64 1.61 1.51 1.59 ****
mean-> 1.55 1.13 1.08 1.10 1.11 1.59 1.58 1.56 1.52 1.54 1.08 1.12 1.01 1.03
Table 7. Efficiency ratios of DEP8(6)(4) vs. NEW8(6) for problems 1 b 14 b and global errors achieved.
Table 7. Efficiency ratios of DEP8(6)(4) vs. NEW8(6) for problems 1 b 14 b and global errors achieved.
Error 1 b 2 b 3 b 4 b 5 b 6 b 7 a 8 b 9 b 10 b 11 b 12 b 13 b 14 b
10 0 *********** 1.13 **
10 1 ********** 0.93 1.14 **
10 2 **** 0.91 ***** 0.98 1.15 **
10 3 * 1.11 ** 0.96 ***** 1.03 1.16 * 0.89
10 4 1.56 1.11 1.03 0.92 1.01 1.54 1.49 1.50 1.54 1.55 1.07 1.17 0.87 0.92
10 5 1.52 1.10 1.07 0.97 1.07 1.55 1.52 1.52 1.55 1.55 1.13 1.18 0.92 0.95
10 6 1.48 1.09 1.11 1.02 1.14 1.56 1.54 1.55 1.56 1.56 1.18 * 0.96 0.98
10 7 1.43 1.08 1.15 1.07 1.20 1.57 1.57 1.57 1.57 1.57 1.23 * 1.02 1.02
10 8 1.39 1.08 1.19 1.12 1.27 1.57 1.60 1.59 1.58 1.58 ** 1.07 1.06
10 9 1.36 1.07 1.24 1.18 * 1.58 1.63 1.61 1.59 1.59 ** 1.13 1.09
10 10 1.32 1.06 1.28 ** 1.59 1.66 1.63 1.60 1.60 ** 1.18 1.13
10 11 1.28 **** 1.60 1.69 1.66 1.61 1.61 ** 1.25 *
10 12 1.25 **** 1.61 1.72 1.68 1.62 1.61 ****
mean-> 1.40 1.09 1.15 1.04 1.08 1.57 1.60 1.59 1.58 1.58 1.08 1.16 1.05 1.01
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jerbi, H.; Omri, M.; Kchaou, M.; Simos, T.E.; Tsitouras, C. Runge-Kutta-Nyström Pairs of Orders 8(6) with Coefficients Trained to Perform Best on Classical Orbits. Mathematics 2022, 10, 654. https://0-doi-org.brum.beds.ac.uk/10.3390/math10040654

AMA Style

Jerbi H, Omri M, Kchaou M, Simos TE, Tsitouras C. Runge-Kutta-Nyström Pairs of Orders 8(6) with Coefficients Trained to Perform Best on Classical Orbits. Mathematics. 2022; 10(4):654. https://0-doi-org.brum.beds.ac.uk/10.3390/math10040654

Chicago/Turabian Style

Jerbi, Houssem, Mohamed Omri, Mourad Kchaou, Theodore E. Simos, and Charalampos Tsitouras. 2022. "Runge-Kutta-Nyström Pairs of Orders 8(6) with Coefficients Trained to Perform Best on Classical Orbits" Mathematics 10, no. 4: 654. https://0-doi-org.brum.beds.ac.uk/10.3390/math10040654

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop