7. *The Bessel equation*

The well-known Bessel equation

$$y''(x) = -y(x) \cdot \frac{1 + 400x^2}{4x^2}.$$

is verified by a theoretical solution of the form [14]

$$y(\mathbf{x}) = l0(10\mathbf{x}) \cdot \sqrt{\mathbf{x}}\_r$$

with *J*<sup>0</sup> the zeroth order Bessel function of the first kind. This equation in also integrated in the interval [0, 10*π*].

#### 8. *The Duffing equation*

Next, we choose the equation

$$\begin{array}{rcl}y''(\mathbf{x})&=&\frac{1}{500}\cdot\cos(1.01\mathbf{x})-y(\mathbf{x})-y(\mathbf{x})^3,\\y(0)&=&0.2004267280699011,y'(0)=0,\end{array}$$

with an approximate analytical solution given in [23],

$$y(\mathbf{x}) \approx \left\{ \begin{array}{l} 6 \cdot 10^{-16} \cos(11.11 \mathbf{x}) + 4.609 \cdot 10^{-13} \cos(9.09 \mathbf{x}) \\ + 3.743495 \cdot 10^{-10} \cos(7.07 \mathbf{x}) + 3.040149839 \cdot 10^{-7} \cos(5.05 \mathbf{x}) \\ + 2.469461432611 \cdot 10^{-4} \cos(3.03 \mathbf{x}) 0.2001794775368452 \cos(1.01 \mathbf{x}) \end{array} \right\}$$

We again solved the above equation in the interval [0, 10*π*].

9. *semi-Linear system.*

The nonlinear problem proposed by Franco and Gomez [24] follows.

$$\begin{array}{rcl} y''(t) &=& \left( \begin{array}{cc} -199 & -198 \\ 99 & 98 \end{array} \right) \cdot y(\mathbf{x}) + \left( \begin{array}{cc} (y\_1 + y\_2)^2 + \sin^2(10\mathbf{x}) - 1 \\ (y\_1 + 2y\_2)^2 - 10^{-6} \sin^2(\mathbf{x}) \end{array} \right), \\ & \mathbf{x} & \in & [0, 10\pi], \end{array}$$

with theoretical solution

$$y(t) = \begin{pmatrix} 2\cos(10x) - 10^{-3}\sin(x) \\ -\cos(10x) + 10^{-3}\sin(x) \end{pmatrix}.$$

10. *Van der Pol oscillator.*

The equation we solved is

$$y'' = 0.1 \cdot (1 - y(\mathbf{x})^2) y'(\mathbf{x}) - y(\mathbf{x}), \\ y(0) = -0.2, \\ y'(0) = 0, \\ \mathbf{x} \in [0, 10\pi],$$

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<sup>−</sup>11. Numbers greater than 1 are in favor of the second pair. ˆ *b*<sup>7</sup> 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 log10 1.855 <sup>≈</sup> 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,

$$\begin{array}{rcl} \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x}} &=& \frac{\partial \boldsymbol{u}}{\partial r}, \; \boldsymbol{u}(\boldsymbol{x},0) = 0, \; \boldsymbol{u}(0,r) = \sin \pi r^2 r^2, \\\ 0 &\le \quad r \le 1, \; \boldsymbol{x} \ge 0, \end{array}$$

is discretized by symmetric differences (with Δ*r* = 1/50) to the system of ODEs

$$
\begin{bmatrix} y\_1' \\ y\_2' \\ \\ y\_{50}' \end{bmatrix} = \frac{1}{2} \cdot \frac{1}{(1/50)} \begin{bmatrix} 0 & -1 & & & \\ 1 & 0 & -1 & & \\ & & 1 & 0 & -1 \\ & & & -1 & 4 & -3 \end{bmatrix} \cdot \begin{bmatrix} y\_1 \\ y\_2 \\ \\ \\ y\_{50} \end{bmatrix}.$$

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.

**Tolerances Problem 10***−***<sup>5</sup> 10***−***<sup>6</sup> 10***−***<sup>7</sup> 10***−***<sup>8</sup> 10***−***<sup>9</sup> 10***−***<sup>10</sup> 10***−***<sup>11</sup> Mean** 1 1.33 1.45 1.60 1.75 1.92 2.11 2.64 1.83 2 1.36 1.49 1.63 1.78 1.96 2.16 3.17 1.94 3 1.37 1.50 1.64 1.80 1.97 2.18 2.58 1.86 4 1.38 1.51 1.65 1.81 1.98 2.19 2.79 1.90 5 1.39 1.51 1.66 1.82 1.99 2.23 2.42 1.86 6 1.41 1.54 1.68 1.85 2.03 2.36 1.80 1.81 7 1.32 1.44 1.57 1.72 1.89 2.07 2.65 1.81 8 1.08 1.36 1.66 1.98 2.33 2.58 2.58 1.94 9 1.52 1.67 1.84 2.00 2.17 2.44 1.73 1.91 10 1.24 1.28 1.34 1.44 1.58 1.75 2.17 1.54

**Table 3.** Efficiency measures ratios of DP5(4) vs. NEW5(4) over the interval [0, 20*π*].

**Figure 1.** Results of DP5(4) vs. NEW5(4) for the Hyperbolic PDE.

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.

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

**Author Contributions:** All authors have contributed equally. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was supported by a Mega Grant from the Government of the Russian Federation within the framework of the federal project No. 075-15-2021-584.

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

