Evolutionary Derivation of Runge–Kutta Pairs of Orders 5(4) Specially Tuned for Problems with Periodic Solutions

The purpose of the present work is to construct a new Runge–Kutta pair of orders five and four to outperform the state-of-the-art in these kind of methods when addressing problems with periodic solutions. We consider the family of such pairs that the celebrated Dormand–Prince pair also belongs. The chosen family comes with coefficients that all depend on five free parameters. These latter parameters are tuned in a way to furnish a new method that performs best on a couple of oscillators. Then, we observe that this trained pair outperforms other well known methods in the relevant literature in a standard set of problems with periodic solutions. This is remarkable since no special property holds such as high phase-lag order or an extended interval of periodicity.


Introduction
The Initial Value Problem (IVP) is with x 0 ∈ R and the vectors y, y ∈ R m . The function f is defined as f : R × R m → R m . Runge-Kutta (RK) pairs are perhaps the most used numerical methods for addressing (1). They usually presented in a so-called Butcher tableau [1,2] as given below. In this type of tableau, we have b T ,b T , c ∈ R s while A ∈ R s×s . Then, the method shares s stages and in case that A is strictly lower triangular, it is evaluated explicitly. The numerical approximations of the solution step from (x n , y n ) ∈ R 1+m to x n+1 = x n + h n Mathematics 2021, 9,2306 2 of 11 by producing two numerical estimations for y(x n+1 ). Namely, y n+1 ∈ R m andŷ n+1 ∈ R m , given by for i = 1, 2, · · · , s. These two approximations y n+1 ∈ R m andŷ n+1 ∈ R m are of algebraic orders p and q < p, respectively. This means that when expanding them in Taylor series, they attain orders O(h p ) and O(h q ), respectively, with h being the proper step-length. Then, a local error estimation is formed in every step and is combined in an algorithm for changing the steplength.
with t a small positive number which is set by the user and is named tolerance. The safety factor 0.8 is in common use in such formulas and offers increased reliability to the results. Whenever n < t, we use the above formula for defining the length of the next step forward. In reverse, when n ≥ t we again use it without advancing the solution in this case and using instead the value h n+1 as a new trial step h n . Information in the issue can be found with details in [3]. As an abbreviation these methods are named RKp(q) airs. Carl David Tolmé Runge [4] and Martin Wilhelm Kutta [5] introduced the methods bearing their names almost in the turning of the 19th century. For almost 60 years the these methods were implemented with constant step sizes. Richardson extrapolation appeared in the meantime [6] and was used in a kind of step control through doubling and halving [7]. Runge-Kutta pairs appeared 60 years ago. The first series of well-known Runge-Kutta pairs of orders 5(4), 6(5), and 8 (7) were presented by Fehlberg [8,9]. In the early 1980s, Dormand and Prince gave their celebrated pairs [10,11].
Non-stiff problems having the form (1) are well suited for being efficiently solved by Runge-Kutta pairs. There is a number of different pairs sharing various orders. We may explain this by the accuracy on demand. Thus, the lesser the accuracy required, the lesser order pairs are chosen. Otherwise, for stringent accuracies at quadruple precision, a high order pair has to be preferred.
Here, we concentrate on RK5(4) pairs which are the first choice when middle tolerances are used. Our special interest is in problems of the form (1) with periodic/oscillatory solutions. In the following, we focus on producing a RK5(4) pair that best address the latter type of problems.
The paper is organized in sections as follows.
Training the coefficients 4.

Theory of Runge-Kutta Pairs for Orders 5(4)
Runge-Kutta pairs of orders five and four are perhaps the most used ones. The coefficients of these pairs have to satisfy 25 order conditions when the almost obligatory holds. Namely, 17 order conditions for the higher order formula and another 8 conditions for the fourth order formula. For a seven stages (i.e., s = 7) pair with an FSAL (First Stage As Last) property there are 28 coefficients after considering (2). Over the years, various techniques for solving this system have been presented. The solutions form different families. Dormand and Prince presented such a family and perhaps the most famous pair of all in [10]. Papakostas and Papageorgiou studied this family extensively [12]. Then, we may choose arbitrarily the coefficients c 2 , c 3 , c 4 , c 5 ,b 7 and produce all the rest coefficients explicitly. The particular pair DP5(4) appeared in [10] shares small principal truncation error coefficients and it is implemented in the builtin function ode45 of MATLAB [13]. We now proceed presenting explicit formulas for the remaining coefficients that only depend on the five free parameters. , , , , , , and finally the FSAL property holds This means that although s = 7, the method wastes only six stages per step and the seventh stage is reused as first stage of the next step.
The question raising now is how to choose the free parameters. Traditionally, the method developers try to minimize some norm for the principal term of the local truncation error, i.e., the terms of h 6 in the residual of Taylor error expansions corresponding to the fifth order method of the underlying RK pair. Another choice is to increase the phase-lag order. This means that we try to reduce the gap in the angle among the numerical and the theoretical solution in a free oscillator [14]. The latter approach is well suited for usage in problems with periodic solutions.

Training the Coefficients
We intent to derive a particular RK5(4) pair that belongs to the family discussed above. The resulting pair has to perform best on harmonic oscillators and other problems with periodic solutions. For achieving this, we will first try to find a pair that performs best on a couple of harmonic oscillators. Then, we will check if this performance expands to other problems with periodic solutions.
Thus, we concentrate on the harmonic oscillator with theoretical solution y(x) = cos µx. This problem can be transformed to a first-order system, and then solved it numerically by a RK5(4) pair picked from the family of solutions we are concerned here. We use tolerance t = 10 −11 and µ = 3 and µ = 7. The choice of µ is tailored by the numerical tests we will present below. The values selected are best placed when µ ∈ [0, 10]. Notice that the above selection of µ is for the training phase. It is hoped that the resulted method will furnish better results for every µ.
We record the number k 1 and k 2 of stages (i.e., function evaluations) needed and the global errors g 1 and g 2 observed over the grid (mesh) in the interval of integration, respectively (i.e for both selections of µ). Then, we form two efficiency measures in the sense that higher values mean lower efficiency. These measures were introduced in [15] for comparing pairs of the same order. Now we may set as fitness function the sum u 1 + u 2 which is meant to be minimized. Thus, the fitness function consists of two runs of an Initial Value Problem. The value u 1 + u 2 changes according to the selection of the free parameters c 2 , c 3 , c 4 , c 5 , andb 7 . We actually do not care at this stage forb 7 , as this coefficient actually affects only the tolerance. Indeed we may chooseb 7 = λ ·b 7 , λ = 0, and set a newb = λb + (1 − λ)b, as the new fourth order formula. Then, the tolerance simply becomes λt.
This idea was originally appeared in [16]. Here, for the minimization of u 1 + u 2 we tried Differential Evolution [17]. We have already got positive results using this approach for methods in integrating of orbits [18,19]. In these latter works, we trained the coefficients of the methods on a Kepler orbit. Then we observed very pleasant results over a set of Kepler orbits as long on other known orbital problems.
DE is an iterative procedure and in every iteration, named generation g, we work with a "population" of individuals c , i = 1, 2, · · · , P with P the population size. An initial population c , i = 1, 2, · · · , P is randomly created in the first step of the method. We have also set as fitness function the measure u 1 + u 2 obtained after two runs of harmonic oscilator. The fitness function is then evaluated for each individual in the initial population. In each generation (iteration) g, a three-phase sequential scheme updates all of the individuals involved. These phases are Differentiation, Crossover, and Selection. For further details in the issue see in [20]. We used MATLAB Software DeMat [21] for implementing the latter technique.
The optimization furnished five values for the parameters. The result is rather robust, i.e., we get almost the same optimal value for u even for neighboring parameters. Thus, we present the selected parameters in 6 significant decimal digits below, The resulting pair is presented in Table 1.  (4) is 216% (interpreting number 3.16) and 180% (interpreting number 2.80), respectively, more expensive for delivering the same accuracy in the two oscillators chosen above. The Euclidean norm of the principal truncation error coefficients for the new pair is T (6) 2 ≈ 2.82 × 10 −4 which is a little smaller than the corresponding value T (6) 2 ≈ 3.99 × 10 −4 for DP5 (4). The absolute stability interval is (−3.55, 0] which is rather in normal magnitude. No extra phase-lag order is observed as bA 4 c = 13128101 9439496880 = 1 840 . See in [14] for details on phase-lag property.
In conclusion, it seems that no extra property is present. The new pair appeared in Table 1 does not possess something interesting. It is hard to believe its special performance after seeing its traditional characteristics.
Other authors have also tried recently to train coefficients of RK methods [22]. However, in that later paper, only second-and third-order methods are considered [4,5] with constant step sizes and over single problems (e.g., Van der Pol). The learning algorithm given there remains to be tested on current and stiffer cases. Our proposal for Differential evolution comes after several papers through the years [16].

Numerical Tests
We tested the following two pairs chosen from the family studied above.
DP5(4) has proven over the years to be perhaps the best pair of orders five and four. Other pairs also exist but the difference with DP5(4) is very small. We do not consider pairs that exploit the knowledge of frequency since this property is not considered here.
All the pairs were run for tolerances 10 −5 , 10 −6 , · · · , 10 −11 , and the efficiency measures of the form (3) were recorded. Notice that actually all the problems are transformed in systems of first order equations.
The problems we tested are the following.
1-5. The model problem with theoretical solution y(x) = cos(µx). This problem was run for five different selections of µ. Thus, when we use µ = 1 the problem is numbered as 1st problem. Then we choose µ = 3 and the problem is numbered as 2nd problem. In consequence when µ = 5 the problem is numbered as 3rd problem, when µ = 7 the problem is numbered as 4th problem, and when µ = 9 the problem is numbered as 5th problem. 6. The inhomogeneous problem with theoretical solution y(x) = cos(10x) + sin(10x) + sin x.

The Bessel equation
The well-known Bessel equation is verified by a theoretical solution of the form [14] y The nonlinear problem proposed by Franco and Gomez [24] follows.
with theoretical solution

Van der Pol oscillator.
The equation we solved is and no analytical solution is known. Thus, the error at grid was estimated using an eighth-order pair from [25] at tolerance 10 −14 .
We calculated seventy (i.e., 7 tolerances times 10 problems) efficiency measures for each pair. We set NEW54 as the reference pair. Then, we divide each efficiency measure of DP5(4) with the corresponding efficiency measure of NEW5(4). The results are recorded in Table 2. The two underlined numbers correspond to the ratios found in the phase of training above as the training was done for problems 2 and 4 and tolerance 10 −11 . Numbers greater than 1 are in favor of the second pair.b 7 was chosen so that the total function evaluations spend for both pairs over all 70 runs is almost equal. The rightmost column shows the mean over all tolerances for each problem. The overall average observed ratio is 1.85 meaning that DP5(4) is about 85% more expensive. This is quite remarkable since much effort has been put over the years for achieving even 10-20% of efficiency [25]. In reverse, this means that about log 10 1.85 5 ≈ 1.34 digits were gained in average at the same costs [15]. Because of the problems used for training, it is obvious that we expect better results when there is a larger linear part and a smaller nonlinear part. However, NEW5(4) outperformed DP5(4) even in the clearly nonlinear problems. We also mention that we got more or less similar results for longer integrations. Especially the results for the interval [0, 20π] are shown in Table 3. The overall average observed is 1.84 and the results slightly differ from those of the previous Table. As a final test, we included a more challenging problem which appears frequently in similar works [14,23], namely, the hyperbolic PDE, ϑu ϑx = ϑu ϑr , u(x, 0) = 0, u(0, r) = sin π 2 r 2 , is discretized by symmetric differences (with ∆r = 1/50) to the system of ODEs The 500th zero of the 20th component in the above problem is reached for x 500 = 33.50999699533, which is found by a very accurate integration at stringent tolerances. We integrated the methods to that point. The results presented as stages vs. error in a semi-log form and are given in Figure 1.  The results are very promising. Some future research may use optimization on a wider range of tolerances and model problems. Perhaps a pair spending a parameter for fulfilling the phase-lag property and then trained for periodic problems would furnish even more interesting results. Of course, application of this technique on other classes of problems is also possible, e.g., orbits.

Conclusions
We proposed the proper training the coefficients of Runge-Kutta pairs of orders five and four in order to perform best on problems with oscillatory solutions. We actually chose a couple of harmonic oscillators, an interval and a tolerance and tried to achieve an outstanding performance there. Thus, we concluded to a new pair which is found to outperform other representatives from this family in a wide range of relevant problems. This pair is supposed to be better than classical DP5(4) for problems with periodic solutions. If there are limitations remain to be clarified by applications in the future research.