Sixth Order Numerov-Type Methods with Coefﬁcients Trained to Perform Best on Problems with Oscillating Solutions

: Numerov-type methods using four stages per step and sharing sixth algebraic order are considered. The coefﬁcients of such methods are depended on two free parameters. For addressing problems with oscillatory solutions, we traditionally try to satisfy some speciﬁc properties such as reduce the phase-lag error, extend the interval of periodicity or even nullify the ampliﬁcation. All of these latter properties come from a test problem that poses as a solution to an ideal trigonometric orbit. Here, we propose the training of the coefﬁcients of the selected family of methods in a wide set of relevant problems. After performing this training using the differential evolution technique, we arrive at a certain method that outperforms the other ones from this family in an even wider set of oscillatory problems.

dealing with the P-stability characteristic, which is important for addressing stiff oscillatory problems. Chawla [4] presented the following modified Numerov scheme, which had the advantage of being evaluated explicitly: with h a steplength that remains constant through the integration of Equation (1): The vectors z k−1 , z k and z k+1 are approximating z(t k − h), z(t k ) and z(t k + h), respectively, while v 1 ∈ R m , v 2 ∈ R m and v 3 ∈ R m are the function evaluations used by the method.
We utilize the known information at mesh according to Since we have already computed f (t k−1 , v 1 ) in the previous step, we need only to evaluate f (t k+1 , v 3 ) and f (t k , v 2 ) every step, and consequently, we spend only two stages per step.
Tsitouras then suggested a Runge-Kutta-Nyström (RKN)-style method [5]. This technique significantly lowered the cost. As a result, just four steps are required to create a sixth-order method, whereas previous implementations required six function evaluations (see [6]).
In the years that followed, our group delved thoroughly into the issue. Tsitouras developed eighth-order methods with nine steps per step in [7]. Ninth-order methods were studied in [8]. Simultaneously, a group of Spanish researchers published some highly interesting work on the same topic [9][10][11].
In the present work, we intend to present a new method for addressing the problems with periodic solutions better. Traditionally, for achieving this, we try to fulfill various properties coming from a simple test equation. The main novelty here is that we will train the available free parameters in a wide set of relevant problems. For this training, we will use the differential evolution technique. It is believed that by using this methodology, we will conclude with a method better tuned for oscillatory problems.

Theory of Two-Step Hybrid Numerov-Type Methods
For numerically addressing Equation (1), higher algebraic order methods are in great demand. We may express t, which is the independent variable, as one of the components of z. As a result, we concentrate without losing generality on the autonomous system z = f (z). Then, an s-stages hybrid Numerov method may presented as [7]: with I s ∈ R s×s the identity matrix, D ∈ R s×s , w T ∈ R s , a ∈ R s the coefficient matrices of the method and For the presentation of the coefficients, we make use of the Butcher tableau [12,13], a D w .
Method (2) can be given using matrices [8]. Since the function evaluations are computed sequentially, these methods are explicit. Thus, D is strictly lower triangular matrix. When s = 5, the associated matrices take the following form: T .
Since f (v 1 ) is known from the previous step, four function evaluations are evaluated each step. For attaining sixth algebraic order, we must cancel the associated truncation error terms (see [14]).
Seventeen parameters are shared by the scheme under examination. Namely, nine entries from the matrix D (i.e., d 31 , d 32 , · · · , d 54 ), five coefficients for vector w and 3 coefficients for vector a. However, in order to obtain 6th order, we must solve 23 condition equations (see Table 5 in [14]).
The parameters are less than the equations. This is a usual problem while developing Runge-Kutta type methods. Using simplifying assumptions is a common way to get around this issue. We proceed setting, Then we spend only the six parameters d 31 , d 32 , d 41 , d 42 , d 51 and d 52 to satisfy the above assumptions. Our profit is that all order conditions, including D · 1 and D · a, are discarded from the relevant list given in [14]. As a result, only 9 order conditions remain to be satisfied by the remaining 11 coefficients.
We select a 3 and a 4 as free parameters. The remainder of the coefficients are computed successively below through a Mathematica [15] listing presented in Figure 1.
For exhaustive information on the derivation of truncation error coefficients, see the review in [14]. Through its link with the so-called T2 rooted trees, Coleman [16] advocated using the B2 series representation of the local truncation error.
A first method from this family was given by Tsitouras [5]. We may write in Mathematica the following lines and derive the method given in there.
In Thus, we verify the efficiency of the algorithm since almost 0.01 seconds are enough for furnishing the coefficients in a Ryzen 9 3900X processor running at 3.79 GHz. Later, Franco [9] chose a 3 = − 1 5 , a 4 = − 2 5 . These were all-purpose methods. In [17], we proposed another approach for selecting a 3 and a 4 that concentrates on the method's behavior in Keplerian type orbits. There we concluded that the choice a 3 = 3 44 , a 4 = − 23 38 furnishes a method that best address the latter type of problems.

Performance of Methods in a Wide Set of Problems with Oscillating Solutions
From the above-mentioned family, we intend to develop a particular hybrid Numerovtype scheme. The resulting method has to perform best on problems with oscillating solutions. For this reason, we have chosen to test the following problems.

The Bessel equation
The wellknown Bessel equation is verified by a theoretical solution of the form, with J 0 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 with an approximate analytical solution given in [14], The three methods F6 [9], M6 [17] and T6 [5] were run for the above problems and for different numbers for steps. The results in [5,9,17] showed the superiority of the latter methods over the older schemes. The global errors over the whole mesh was recorded in Table 1. Actually, we presented the errors in the form of the accurate digits observed. A final row with the mean value is also given in Table 1.
In these 8 problems and for the 32 runs carried, it seems that T6 performed best. The question raised now is if we can do even better.

Phase-Lag and Amplification Errors
At first, we select a method of high 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 [18]. The latter approach is well suited for use in problems with periodic solutions. Thus, after considering the test problem z = −λz, and applying Method (3), we verify as phase-lag the expression: A sixth-order method shares sixth phase-lag order. Then after expanding with respect to τ = hλ, we obtain ρ = ρ 8 τ 8 + ρ 10 τ 10 + · · · The equations for eighth and tenth phase-lag order are , and , respectively. The only acceptable solution of ρ 8 = ρ 10 = 0 is a 3 = 16 15 and a 4 = 1371 245 . However, we can not use such coefficients being so far away from the interval of interest [−1, 1]. Thus, we may draw back and accept only ρ 8 = 0 by setting We name this method PL8. Another choice is the elimination of amplification errors. This is the distance from the orbit of the theoretical solution of a free oscillator. It is given as Expanding with respect to τ, we conclude to the exact form, Unfortunately, we may not satisfy simultaneously σ ≡ 1 and ρ 8 = 0 since we arrive at coefficients with indeterminate values. Thus, we may admit only σ ≡ 1 by setting a 3 = − 1 2 and a 4 = 7 11 . We name this method σ1. Another interesting property is P-stability [2,3]. Then, we have to satisfy σ ≡ 1 along Only implicit methods may address these two requirements simultaneously.

Training the Free Parameters in a Wide Set of Periodic Problems
Our current project's initial concept is based on [19]. After choosing the free parameters a 3 , a 4 , we obtain a method named NEW6 and form another column in Table 1 for it. The average value r obtained after the 32 runs may serve as a fitness measure and meant to be maximized. For the maximization process, we applied the differential evolution technique [20].
DE is an iterative procedure, and in every iteration, named generation g, we work with a "population" of individuals a , i = 1, 2, · · · , N, with N the population size.
An initial population a 4i , i = 1, 2, · · · , N is randomly created in the first step of the method. We have also set the measure r as the fitness function, i.e., the average of accurate digits after the 32 runs mentioned above. The fitness function is then evaluated for each individual in the initial population. In each generation (iteration) g, a three-phases sequential scheme updates all of the individuals involved. These phases are Differentiation, Crossover and Selection.
We used MATLAB [21] software DeMat [22] for implementing the latter technique. Indeed, we manage to produce an improvement by choosing: The coefficients of the new method in matrix forms are given below, which are suitable for double precision computations. For this method, we obtained r ≈ 7.75, which is a very impressive result. Actually, we obtained many methods with r > 7.7 since there seems to exist a small area of pairs a 3 , a 4 , where r attains high values. We also remark that for Selection (4), the amplification is σ = 1 and for phase lag ρ = O(v 8 ) holds, i.e., ρ 8 = 0, and no special property is satisfied.
We run the methods constructed for addressing periodic problems in the eight problems listed in the third section. We summarize the results in Table 2. It is clear from this Table that NEW6 performs better than all methods referred until now. Namely, F6, M6, T6, PL8 and σ1.
Other authors have also tried recently to train coefficients of RK methods [23]. However, in that later paper, only second-and third-order methods are considered [24,25] with constant step sizes and over single problems (e.g., Van der Pol). The learning al-gorithm given there remains to be tested on current and stiffer cases. Our proposal for differential evolution comes after several papers through the years [19].

Numerical Results
Method NEW6 was produced to perform best on problems 1-8 listed in Section 3. In the tests recorded in Tables 1 and 2, it was meant to outperform other methods for the intervals and steps used there.
Thus, we intend to test NEW6 in a different set of problems, intervals and number of steps. In this direction, we run again problems 1-8 to the longer interval [0, 20π]. We name these problems now 1a, 1b, · · · , 8a. In addition, we included two nonlinear problems more.
10. Two coupled oscillators with different frequencies.
The problem is characterized by the equations [27], We also integrated this problem into [0, 20π], but no analytical solution is available. For an estimation of the error in the grid points, we used a Runge-Kutta-Nyström method [28] with very stringent tolerance.
Finally, we consider the linearized wave equation, which is a rather large-scale problem [14], with the theoretical solution We semi-discretisize ϑ 2 u ϑx 2 with fourth order symmetric differences at internal points and one sided differences of the same order at the boundaries (including the knowledge of ϑu ϑx there) and conclude with the system: By choosing ∆x = 5, we arrive at a constant coefficients system with N = 20. The results for this problem were dominated by the semi-discretization errors.
We run these 11 problems for various numbers of steps and tabulated the results in Table 3. There, we included results with other state-of-the-art methods considered in the area of sixth-order Numerov-type (i.e., including off step points) methods. It is obvious from there that NEW6 outperformed all other methods from the literature by a considerable distance.

Conclusions
The main points of our research was the following.
• We considered a family of a sixth-order hybrid two-step scheme that shares the lowest number of stages, and the main novelty is suggesting a method for selecting proper free parameters. • The parameters of the new method were chosen after testing their performance in a large set of periodic problems. • The best choice was found using the differential evolution method. In a wide range of problems with oscillating solutions, the developed scheme significantly outperformed other methods from the same or other families. • The presented method is tuned for problems with periodic solutions,F especially when these problems share a large linear part.