**2. 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

$$A \cdot \mathfrak{e} = \mathfrak{e}, \mathfrak{e} = [1, 1, \mathfrak{h}, \dots, 1] \in \mathbb{R}^s \tag{2}$$

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*<sup>7</sup> 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.

*<sup>b</sup>*<sup>3</sup> <sup>=</sup> *<sup>c</sup>*4(<sup>5</sup> <sup>−</sup> <sup>10</sup>*c*5) + <sup>5</sup>*c*<sup>5</sup> <sup>−</sup> <sup>3</sup> 60(*c*<sup>3</sup> − 1)*c*3(*c*<sup>3</sup> − *c*4)(*c*<sup>3</sup> − *c*5) , *<sup>b</sup>*<sup>4</sup> <sup>=</sup> <sup>5</sup>*c*3(2*c*<sup>5</sup> <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>5</sup>*c*<sup>5</sup> <sup>+</sup> <sup>3</sup> 60(*c*<sup>4</sup> − 1)*c*4(*c*<sup>3</sup> − *c*4)(*c*<sup>4</sup> − *c*5) , *<sup>b</sup>*<sup>5</sup> <sup>=</sup> <sup>5</sup>*c*3(2*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>5</sup>*c*<sup>4</sup> <sup>+</sup> <sup>3</sup> 60(*c*<sup>5</sup> − 1)*c*5(*c*<sup>3</sup> − *c*5)(*c*<sup>5</sup> − *c*4) , *<sup>b</sup>*<sup>6</sup> <sup>=</sup> <sup>5</sup>*c*3(*c*4(6*c*<sup>5</sup> <sup>−</sup> <sup>4</sup>) <sup>−</sup> <sup>4</sup>*c*<sup>5</sup> <sup>+</sup> <sup>3</sup>) <sup>−</sup> <sup>20</sup>*c*4*c*<sup>5</sup> <sup>+</sup> <sup>15</sup>*c*<sup>4</sup> <sup>+</sup> <sup>15</sup>*c*<sup>5</sup> <sup>−</sup> <sup>12</sup> <sup>60</sup>(*c*<sup>3</sup> <sup>−</sup> <sup>1</sup>)(*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>)(*c*<sup>5</sup> <sup>−</sup> <sup>1</sup>) , ˆ *b*<sup>3</sup> = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ ⎛ ⎜⎜⎝ 10(6ˆ *<sup>b</sup>*<sup>7</sup> <sup>−</sup> <sup>1</sup>)*c*<sup>2</sup> 3*c*4 <sup>+</sup>*c*3(−8ˆ *b*7(7*c*<sup>4</sup> + 1) +8*c*<sup>4</sup> + 1) + 2(8ˆ *b*<sup>7</sup> − 1)*c*<sup>4</sup> ⎞ ⎟⎟⎠ 5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) <sup>−</sup>20*c*4*c*<sup>5</sup> <sup>+</sup> <sup>15</sup>*c*<sup>4</sup> <sup>+</sup> <sup>15</sup>*c*<sup>5</sup> <sup>−</sup> <sup>12</sup> 5(*c*3−1) 10*c*<sup>2</sup> 3*c*4 −*c*3(8*c*<sup>4</sup> + 1) + 2*c*<sup>4</sup> <sup>−</sup>12ˆ *b*7(*c*<sup>4</sup> − 1)(*c*<sup>5</sup> − 1) + 2*c*4(3*c*<sup>5</sup> − 2) − 4*c*<sup>5</sup> + 3 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ <sup>12</sup>*c*3(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*4)(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*5) , ˆ *b*<sup>4</sup> = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ − ⎛ ⎜⎜⎝ 10(6ˆ *<sup>b</sup>*<sup>7</sup> <sup>−</sup> <sup>1</sup>)*c*<sup>2</sup> 3*c*4 <sup>+</sup>*c*3(−8ˆ *b*7(7*c*<sup>4</sup> + 1) + 8*c*<sup>4</sup> + 1) +2(8ˆ *b*<sup>7</sup> − 1)*c*<sup>4</sup> ⎞ ⎟⎟⎠ 5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) <sup>−</sup>20*c*4*c*<sup>5</sup> <sup>+</sup> <sup>15</sup>*c*<sup>4</sup> <sup>+</sup> <sup>15</sup>*c*<sup>5</sup> <sup>−</sup> <sup>12</sup> 5(*c*4−1) 10*c*<sup>2</sup> 3*c*4 −*c*3(8*c*<sup>4</sup> + 1) + 2*c*<sup>4</sup> +12ˆ *b*7(*c*<sup>3</sup> − 1)(*c*<sup>5</sup> − 1) − 2*c*3(3*c*<sup>5</sup> − 2) + 4*c*<sup>5</sup> − 3 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ <sup>12</sup>*c*4(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*4)(*c*<sup>4</sup> <sup>−</sup> *<sup>c</sup>*5) , *b*<sup>5</sup> = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ ⎧ ⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎩ ⎛ ⎜⎜⎝ 10(6ˆ *<sup>b</sup>*<sup>7</sup> <sup>−</sup> <sup>1</sup>)*c*<sup>2</sup> 3*c*4 <sup>+</sup>*c*3(−8ˆ *b*7(7*c*<sup>4</sup> + 1) + 8*c*<sup>4</sup> + 1) +2(8ˆ *b*<sup>7</sup> − 1)*c*<sup>4</sup> ⎞ ⎟⎟⎠ 5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) <sup>−</sup>20*c*4*c*<sup>5</sup> <sup>+</sup> <sup>15</sup>*c*<sup>4</sup> <sup>+</sup> <sup>15</sup>*c*<sup>5</sup> <sup>−</sup> <sup>12</sup> 5(*c*5−1) 10*c*<sup>2</sup> 3*c*4 −*c*3(8*c*<sup>4</sup> + 1) + 2*c*<sup>4</sup> ⎫ ⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎭ <sup>−</sup>12ˆ *b*7(*c*<sup>3</sup> − 1)(*c*<sup>4</sup> − 1) + 2*c*3(3*c*<sup>4</sup> − 2) − 4*c*<sup>4</sup> + 3 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ <sup>12</sup>*c*5(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*5)(*c*<sup>4</sup> <sup>−</sup> *<sup>c</sup>*5) , ˆ *b*<sup>6</sup> = − 10(6ˆ *<sup>b</sup>*<sup>7</sup> <sup>−</sup> <sup>1</sup>)*c*<sup>2</sup> <sup>3</sup>*c*<sup>4</sup> <sup>+</sup> *<sup>c</sup>*3(−8ˆ *b*7(7*c*<sup>4</sup> + 1) + 8*c*<sup>4</sup> + 1) + 2(8ˆ *b*<sup>7</sup> − 1)*c*<sup>4</sup> · (5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) − 20*c*4*c*<sup>5</sup> + 15*c*<sup>4</sup> + 15*c*<sup>5</sup> − 12) & 60(*c*<sup>3</sup> − 1)(*c*<sup>4</sup> − 1)(*c*<sup>5</sup> − 1) 10*c*<sup>2</sup> <sup>3</sup>*c*<sup>4</sup> − *c*3(8*c*<sup>4</sup> + 1) + 2*c*<sup>4</sup> ,

ˆ

*<sup>a</sup>*3,2 <sup>=</sup> *<sup>c</sup>*<sup>2</sup> 3 2*c*<sup>2</sup> , *<sup>a</sup>*4,2 <sup>=</sup> *<sup>c</sup>*<sup>2</sup> <sup>4</sup>(3*c*<sup>3</sup> − 2*c*4) 2*c*2*c*<sup>3</sup> , *<sup>a</sup>*4,3 <sup>=</sup> *<sup>c</sup>*<sup>2</sup> <sup>4</sup>(*c*<sup>4</sup> − *c*3) *c*2 3 , *<sup>a</sup>*5,2 <sup>=</sup> *<sup>c</sup>*<sup>5</sup> 15*c*<sup>2</sup> <sup>3</sup>*c*4(2*c*<sup>5</sup> − 1) + *c*<sup>3</sup> *c*4 <sup>6</sup> <sup>−</sup> <sup>20</sup>*c*<sup>2</sup> 5 + (3 − 5*c*5)*c*<sup>5</sup> + 2*c*4*c*5(5*c*<sup>5</sup> − 3) <sup>2</sup>*c*2*c*3(5*c*3(2*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>5</sup>*c*<sup>4</sup> <sup>+</sup> <sup>3</sup>) , *a*5,3 = − *c*5(*c*<sup>3</sup> − *c*5) - 10*c*<sup>2</sup> <sup>3</sup>*c*4(2*c*<sup>5</sup> − 1) + *c*<sup>3</sup> - <sup>−</sup>5*c*<sup>2</sup> <sup>4</sup>(4*c*<sup>5</sup> − 3) +*c*4(4 − 15*c*5) + 2*c*<sup>5</sup> + 2*c*<sup>2</sup> <sup>4</sup>(5*c*<sup>5</sup> − 3) 2*c*<sup>2</sup> <sup>3</sup>(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*4)(5*c*3(2*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>5</sup>*c*<sup>4</sup> <sup>+</sup> <sup>3</sup>) , *<sup>a</sup>*5,4 <sup>=</sup> (5*c*<sup>3</sup> <sup>−</sup> <sup>2</sup>)*c*5(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*5)(*c*<sup>4</sup> <sup>−</sup> *<sup>c</sup>*5) 2*c*4(*c*<sup>3</sup> − *c*4)(5*c*3(2*c*<sup>4</sup> − 1) − 5*c*<sup>4</sup> + 3) , *<sup>a</sup>*6,2 <sup>=</sup> <sup>15</sup>*c*<sup>2</sup> <sup>3</sup>*c*4(2*c*<sup>5</sup> − 1) + *c*3(*c*4(16 − 30*c*5) − 5*c*<sup>5</sup> + 3) + 2*c*4(5*c*<sup>5</sup> − 3) 2*c*2*c*3(5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) − 20*c*4*c*<sup>5</sup> + 15*c*<sup>4</sup> + 15*c*<sup>5</sup> − 12) , *a*6,3 = − (*c*<sup>3</sup> − 1) · ⎛ ⎝ <sup>−</sup>*c*<sup>2</sup> 3 5*c*<sup>2</sup> <sup>4</sup>(4*c*<sup>5</sup> <sup>−</sup> <sup>3</sup>) + <sup>20</sup>*c*4*c*<sup>2</sup> <sup>5</sup> + *c*<sup>4</sup> − 2 +*c*<sup>3</sup> *c*2 <sup>4</sup>(25*c*<sup>5</sup> − 16) + *c*<sup>4</sup> 40*c*<sup>2</sup> <sup>5</sup> − 45*c*<sup>5</sup> + 16 − 2 5*c*<sup>2</sup> <sup>5</sup> − 7*c*<sup>5</sup> + 3 10*c*<sup>3</sup> <sup>3</sup>*c*4(2*c*<sup>5</sup> <sup>−</sup> <sup>1</sup>) + <sup>2</sup>*c*<sup>2</sup> <sup>4</sup>(3 − 5*c*5)*c*<sup>5</sup> ⎞ ⎠ 2*c*<sup>2</sup> <sup>3</sup>(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*4)(*c*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*5)(5*c*3(*c*4(6*c*<sup>5</sup> <sup>−</sup> <sup>4</sup>) <sup>−</sup> <sup>4</sup>*c*<sup>5</sup> <sup>+</sup> <sup>3</sup>) <sup>−</sup> <sup>20</sup>*c*4*c*<sup>5</sup> <sup>+</sup> <sup>15</sup>*c*<sup>4</sup> <sup>+</sup> <sup>15</sup>*c*<sup>5</sup> <sup>−</sup> <sup>12</sup>) , *<sup>a</sup>*6,4 <sup>=</sup> (*c*<sup>3</sup> <sup>−</sup> <sup>1</sup>)(*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>) 5*c*<sup>3</sup> *<sup>c</sup>*<sup>4</sup> <sup>−</sup> <sup>4</sup>*c*<sup>2</sup> <sup>5</sup> + 5*c*<sup>5</sup> − 2 − 2 *<sup>c</sup>*<sup>4</sup> <sup>−</sup> <sup>5</sup>*c*<sup>2</sup> <sup>5</sup> + 7*c*<sup>5</sup> − 3 2*c*4(*c*<sup>3</sup> − *c*4)(*c*<sup>4</sup> − *c*5)(5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) − 20*c*4*c*<sup>5</sup> + 15*c*<sup>4</sup> + 15*c*<sup>5</sup> − 12) , *<sup>a</sup>*6,5 <sup>=</sup> (*c*<sup>3</sup> <sup>−</sup> <sup>1</sup>)(*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>)(*c*<sup>5</sup> <sup>−</sup> <sup>1</sup>)(5*c*3(2*c*<sup>4</sup> <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>5</sup>*c*<sup>4</sup> <sup>+</sup> <sup>3</sup>) *c*5(*c*<sup>3</sup> − *c*5)(*c*<sup>4</sup> − *c*5)(5*c*3(*c*4(6*c*<sup>5</sup> − 4) − 4*c*<sup>5</sup> + 3) − 20*c*4*c*<sup>5</sup> + 15*c*<sup>4</sup> + 15*c*<sup>5</sup> − 12) , *<sup>b</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>b</sup>*<sup>3</sup> <sup>−</sup> *<sup>b</sup>*<sup>4</sup> <sup>−</sup> *<sup>b</sup>*<sup>5</sup> <sup>−</sup> *<sup>b</sup>*6, <sup>ˆ</sup> *<sup>b</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>ˆ</sup> *<sup>b</sup>*<sup>3</sup> <sup>−</sup> <sup>ˆ</sup> *<sup>b</sup>*<sup>4</sup> <sup>−</sup> <sup>ˆ</sup> *<sup>b</sup>*<sup>5</sup> <sup>−</sup> <sup>ˆ</sup> *<sup>b</sup>*<sup>6</sup> <sup>−</sup> <sup>ˆ</sup> *b*7, *a*<sup>21</sup> = *c*2, *a*<sup>31</sup> = *c*<sup>3</sup> − *a*32, *a*<sup>41</sup> = *c*<sup>4</sup> − *a*<sup>42</sup> − *a*43, *a*<sup>51</sup> = *c*<sup>5</sup> − *a*<sup>52</sup> − *a*<sup>53</sup> − *a*54, *a*<sup>61</sup> = *c*<sup>6</sup> − *a*<sup>62</sup> − *a*<sup>63</sup> − *a*<sup>64</sup> − *a*65,

and finally the FSAL property holds

$$a\_{\overline{\gamma}} = b\_{\overline{\gamma}'} \mathbf{j} = 1\_{\prime} 2\_{\prime} \cdots \mathbf{j}\_{\prime} \mathbf{6}\_{\cdot}$$

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*<sup>6</sup> 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.

## **3. 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

$$y'' = -\mu^2 y, \; y(0) = 1, \; y'(0) = 0, \pi \in [0, 10\pi]\_\prime$$

with theoretical solution *y*(*x*) = cos *μx*. This problem can be transformed to a first-order system,

$$
\begin{bmatrix} y\_1'\\ y\_2' \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -\mu^2 & 0 \end{bmatrix} \begin{bmatrix} y\_1\\ y\_2 \end{bmatrix} \prime \ y\_1(0) = 1, y\_2(0) = 0, \mathbf{x} \in [0, 10\pi]\_\prime
$$

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−<sup>11</sup> 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*<sup>1</sup> and *k*<sup>2</sup> of stages (i.e., function evaluations) needed and the global errors *g*<sup>1</sup> and *g*<sup>2</sup> observed over the grid (mesh) in the interval of integration, respectively (i.e for both selections of *μ*). Then, we form two efficiency measures

$$
u\_{\rangle} = k\_{\rangle} \cdot g\_{\rangle}^{1/5}, \, j = 1, 2 \tag{3}$$

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*<sup>1</sup> + *u*<sup>2</sup> which is meant to be minimized. Thus, the fitness function consists of two runs of an Initial Value Problem. The value *u*<sup>1</sup> + *u*<sup>2</sup> changes according to the selection of the free parameters *c*2, *c*3, *c*4, *c*5, and ˆ *b*7. We actually do not care at this stage for ˆ *b*7, as this coefficient actually affects only the tolerance. Indeed we may choose

$$
\check{b}\_7 = \lambda \cdot \hat{b}\_{7'} \lambda \neq 0\_\prime
$$

and set a new ˜ *b* = *λ*ˆ *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*<sup>1</sup> + *u*<sup>2</sup> 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* (*g*) <sup>2</sup> , *c* (*g*) <sup>3</sup> , *c* (*g*) <sup>4</sup> , *c* (*g*) <sup>5</sup> , <sup>ˆ</sup> *b* (*g*) 7 *i* , *i* = 1, 2, ··· , *P* with *P* the population size. An initial population *c* (0) <sup>2</sup> , *c* (0) <sup>3</sup> , *c* (0) <sup>4</sup> , *c* (0) <sup>5</sup> , <sup>ˆ</sup> *b* (0) 7 *i* , *i* = 1, 2, ··· , *P* is randomly created in the first step of the method. We have also set as fitness function the measure *u*<sup>1</sup> + *u*<sup>2</sup> 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,

$$c\_2 = \frac{6618}{21991},\ c\_3 = \frac{3679}{11497},\ c\_4 = \frac{25691}{30789},\ c\_5 = \frac{5444}{5589},\ b\_7 = \frac{11}{400}.$$

The resulting pair is presented in Table 1.


**Table 1.** Coefficients of NEW5(4) pair, accurate for double precision computations.

For the above selection of free parameters, we got

$$\mu\_1^{\text{NEW54}} = 88.37, \; \mu\_2^{\text{NEW54}} = 284.89,$$

while for DP5(4) we observed

*u*DP54 <sup>1</sup> = 279.28, *<sup>u</sup>*DP54 <sup>2</sup> = 797.55

i.e.,

$$\frac{\mu\_1^{\text{DP54}}}{\mu\_1^{\text{NEW54}}} + \frac{\mu\_2^{\text{DP54}}}{\mu\_2^{\text{NEW54}}} = \frac{279.28}{88.37} + \frac{797.55}{284.89} \approx 3.16 + 2.80 = 5.96\%$$

meaning that DP5(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 *<sup>T</sup>*(6) <sup>2</sup> <sup>≈</sup> 2.82 <sup>×</sup> <sup>10</sup>−<sup>4</sup> which is a little smaller than the corresponding value *<sup>T</sup>*(6) <sup>2</sup> <sup>≈</sup> 3.99 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 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* = <sup>13128101</sup> <sup>9439496880</sup> <sup>=</sup> <sup>1</sup> <sup>840</sup> . 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].

## **4. 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<sup>−</sup>5, 10−6, ··· , 10<sup>−</sup>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*

$$y''(\mathbf{x}) = -\mu^2 y(\mathbf{x}), \ y(0) = 1, y'(0) = 0, \mathbf{x} \in [0, 10\pi]\_\*$$

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*

$$y'''(\mathbf{x}) = -100y(\mathbf{x}) + 99\sin\mathbf{x}, \\ y(0) = 1, \\ y'(0) = 11, \mathbf{x} \in [0, 10\pi]\_q$$

with theoretical solution *y*(*x*) = cos(10*x*) + sin(10*x*) + sin *x*.
