A Singularly P-Stable Multi-Derivative Predictor Method for the Numerical Solution of Second-Order Ordinary Differential Equations

In this paper, a symmetric eight-step predictor method (explicit) of 10th order is presented for the numerical integration of IVPs of second-order ordinary differential equations. This scheme has variable coefficients and can be used as a predictor stage for other implicit schemes. First, we showed the singular P-stability property of the new method, both algebraically and by plotting the stability region. Then, having applied it to well-known problems like Mathieu equation, we showed the advantage of the proposed method in terms of efficiency and consistency over other methods with the same order.


Introduction
The most natural mathematical form of many general laws of nature (in physics, chemistry, biology, and astronomy) is found in the language of differential equations. Differential equations are generally categorized into two groups: ODEs and PDEs. Traces of ordinary differential equations can be found in various fields of mathematics, natural, or social sciences. Geometry, various engineering fields including analytical mechanics and electrical engineering, geology, physics, chemistry (in the analysis of nuclear chain reactions), biology (in the modeling of infectious diseases and genetic changes), ecology (in population modeling), and economy (in the modeling of dividend and stock price changes) are some of the scientific branches in which the ordinary differential equations play an essential role [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Since most of the differential equations that provide a relatively accurate model of the target phenomena have a complex and nonlinear form, finding an analytical solution for these problems is usually very difficult or even impossible. The lack of analytical solutions for these complex and nonlinear equations has led to the introduction and development of numerical solutions [16][17][18][19][20][21]. Speed, accuracy, and precision are the most common evaluation metrics used to evaluate these numerical solutions. Due to significant advances in the processing capabilities and speed of the processors and computers in the late 20th century, numerical solutions became more prevalent and continue to expand [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. The main topic presented in this paper is about the approximate solution of the following equation: y = f (x, y), y(x 0 ) = y 0 , y (x 0 ) = y 0 , (1) where f is a sufficiently continuous and derivative function of any order. Approximate methods for solving Equation (1) are divided into two main categories: • Methods with unknown constant coefficients, • Methods with variable coefficients (coefficients are a function of the frequency of the problem).
Computational schemes involving a parameter proposed by Neta [24,37], Steifel and Bettis [38], and Alolyan et al. [4] belong in the first class. Mehdizadeh Khalsaraei and Shokri [22], Lambert and Watson [14], Neta [25], Van Daele et al. [40], Vigo-Aguiar et al. [41,42], Franco [16], Martin-Vaquero and Vigo-Aguiar [15], Wang et al. [43] have designed schemes to solve problems of the second class. Because the methods with fixed coefficients are mainly not able to control the quality of the answers the methods with variable coefficients produce higher quality answers (i.e., the solutions produced from the approximate methods of oscillating and periodic quality memory are the main answers). Therefore, in this paper, we produce a new method from the category of methods with variable coefficients. Given that in most cases the answers to Equation (1) are oscillating and periodic, approximate methods that are generated using trigonometric coefficients (functions of sine and cosine) in addition to production high-order algebraic methods also therefore maintain the qualitative properties of the answer. Therefore, with less operation volume, they produce better answers than other methods. Of course, this is shown in one of the most important examples called the Matheiu equation. This equation, which is one of the important ones in physical sciences, has an oscillating answer and is one of the stiff problems in numerical analysis that every method is not able to solve with the desired accuracy. We have solved this equation with the generated method and obtained better solutions than would be obtained by employing other methods. This method will be discussed in the section on numerical results.

Preliminaries
To check the characteristics of the 2k-step symmetric scheme below the test equation of Lambert and Watson is then used in all articles to examine the properties of linear multi-step methods. In addition, the 2k-step multiderivative symmetric multistep methods are defined as follows: After applying Method (2) in Equation (3), we have k ∑ j=1 A j (ν) y n+j + y n−j + A 0 (ν)y n = 0.
A polynomial is now generated using (5) as follows: wherein the position of roots of the polynomial (6) will indicate the main characteristics of method (1).

Definition 2.
If the periodicity interval is (0, ∞) − S, the P-stable scheme (4) is called singularly P-stable, in which S is a finite set of points. (5), the phase-lag is given by

Definition 3. For any multistep scheme with difference Equation
Therefore, the order of the phase-lag will be δ if PL = O ν δ+1 as ν → ∞.

Theorem 1.
The symmetric 2k-step scheme with difference Equation (5) has phase-lag order δ and phase-lag constant ζ given by Proof. Please see [36].

Development
To approximate the solution of Equation (1), we propose a new method in the following form: where k = 4, f n±i = y n±i , g n±i = y (4) n±i , and w n±i = y (6) n±i . In addition, and α i , β i , γ i and η i , i = 0, 1, 2, 3, 4 are real constant unknown coefficients. It is natural that all the properties and characteristics of any numerical method are deduced from its coefficients. However, the coefficients must be carefully generated in order to achieve superior properties. In this paper, several goals are pursued simultaneously so that we can produce a better and more accurate method. Not only do we want to control the interval of periodicity so that (0, +∞) is produced, but we also intend to increase LTE. Applying the new scheme (9), the characteristic equation for the scheme (9) will be as in which where v = ωh. Now, to control the interval of periodicity, we will generate the coefficients so that the characteristic equation of the method is as follows: To this end, the polynomials A 1 (v), A 2 (v) and A 3 (v) of the characteristic function must be equal to zero, which produce three equations. Hence, we have to have Now, by applying the direct formula for the computation of the phase-lag (8) for k = 4 and for A j , j = 0, 1, 2, 3, 4. This leads to the following equation: Since the new method has nine free parameters, we need six more independent linear equations to obtain unique coefficients. Suppose that the phase-lag and its first five derivatives are zero. Then, we have The relations (12) and (13) are nine equations with nine unknowns. Using Maple 17, we have solved this system and the obtained coefficients. The coefficients of the new scheme are given in Appendix A. Whenever, for some |v|, the denominator of the coefficients tend to zero, the following expansions can be considered: The behavior of the coefficients is shown in Figures 1-5.

Theorem 2 ([10]). A necessary condition for convergence of a linear multistep method is that it
should be zero-stable.

Theorem 3 ([10]). The order of a convergent linear multistep method is at least 1.
This condition is called the condition of consistency. In the other words, a linear multistep method is called consistent if it is of order p ≥ 1. Note that the new method proposed in this paper is of algebraic orders ten and, as a result, is consistent. Note that, since the algebraic order of the new method is ten, the coefficient of h 0 in its local truncation error is zero. Therefore, +1 is a root of the first characteristic polynomial (that is, its principal root). Similarly, it is easy to see that the new method satisfies zero-stability.

Theorem 4 ([10]). A linear multistep method is convergent if and only if it satisfies the zerostability and consistency.
According to Theorem 1, the new method is convergent.

Periodicity Analysis
In this section, we discuss and analyze the stability and periodicity properties of the new method (9). To study the stability of a symmetric multistep method, we apply our new method to the test equation Application of the new method (9) to the scalar test (14) leads to the following difference equation Given this, to produce the stability polynomial of the method (to draw the stability region), two frequencies ω and φ have been used, which in the general case ω = φ. Therefore, the generated method can be easily used for two-frequency problems. The method is called P-stable when all parts of the first quarter of the coordinate system are colored. More precisely, the colored parts in the figure related to the stability zone ( Figure 6) show the stability region of the method, while the white parts show the instability region of the method. Now, since the main purpose of this paper is to investigate problems with a frequency, it is sufficient to consider the bisector of the first quarter of the Cartesian coordinate system by equating the frequencies ω and φ as ω = φ. Therefore, the term singular P-stability is used instead of P-stability. It is a singularly P-stable method in which the bisector of the first quarter of the coordinate system is completely inside the colored area. According to Figure 6, the new method presented in this paper is singularly P-stable. Of course, this is proved algebraically in the next theorem.
Theorem 5. The new tenth algebraic order explicit eight-step scheme (9) is singularly P-stable.
Proof. If we assume s = v, the characteristic equation of (9) is given by: Since cos(4v) = 8 cos(v) 4 − 8 cos(v) 2 + 1, the characteristic equation of the new method may be rewritten as Then, by using simple calculation, the characteristic roots are obtained as Thus, when s = v, the periodicity interval of the new scheme is (0, ∞). Therefore, the new scheme is P-stable. Thus, the method is singularly P-stable.

Example 2.
Consider the following pseudo-test problem, which is a highly volatile: where y(0) = 1, y (0) = 0 and ω = 5. The exact answer to this problem is y(x) = 1 3 (cos(x) + cos(3x) + cos(5x)). The oscillating behavior of the answer to this equation is shown in Figure 10 (left). We have approximated this equation in three ways. These methods are: the new method, Simos [34] and method (2.11) in Ananthakrishnaiah [6]. The absolute errors for this problem are shown in Figures 10-12.  . The absolute errors for the test-like problem (16) with h = π 40 using the methods: Ananthakrishnaiah [6] (left) and Simos [34] (right) for Example 2. Example 3. We consider the following almost periodic problem studied by Stiefel and Bettis [38]: where i 2 = −1, z(0) = 1, z (0) = 0.995i and z ∈ C. Here, a theoretical solution represents a motion of the perturbed circular orbit in the complex plane and can be written as z(x) = (1 − 0.0005i) exp(ix). This problem has been solved numerically by using the step size of h = π, π/2, π/3, · · · , π/6 in the range of 0 < x < 40π (which corresponds to 20 orbits of the point z(x)). For comparison purposes, Table 1 presents the errors in the computed values of z(40π) obtained by Ananthakrishnaiah [6], Simos [34], Wang [43], and the new method.
where B = 0.002, ω = 1.01 and x ∈ [0, 10π]. We use the following exact solution for (18) from [24], To integrate this equation by a multiderivative method, one needs the values of y , which occur in calculating y (4) . These higher order derivatives can all be expressed in terms of y(x) and y (x) through (18), i.e., The absolute errors at x = π, 2π, 4π, 6π, 8π, and 10π for the new method, in comparison with methods of Hairer [19], (2.12) in Ananthakrishnaiah [6], Simos [34] and Neta and Fukushima [26] are also listed in Table 2 for comparison with h = π 5 . Example 5. In our last example, we consider the Mathieu differential equation where y(0) = 1, y (0) = 0 and α = 0.1. The frequency is π/5 (see Gautschi [18]). Mathieu equation appears in physical problems involving elliptical shapes or periodic potentials. We will solve the equation for 0 ≤ x ≤ 50 and tabulate at several points. The results are compared to the values obtained from Maple using the function MathieuC(100,5,x). See the plot of the function MathieuC(100,5,x) on the range 0 ≤ x ≤ 50 in Figure 13 (left). In order to integrate this equation by an Obrechkoff method, one needs the values of y , which occur in calculating y (4) and y (6) . For this purpose, we have used the differentiation formula given by: (305y n+1 − 544y n + 239y n−1 ) + h 1980 −5728y n − 571y n−1 + 119y n+1 The absolute error is computed with h = 0.01 at the end point of integration using the methods: Ananthakrishnaiah [6], Simos [34], Wang [43], and the new method, and these are shown in Figures 13-15.

Conclusions
In this paper, by vanishing of the phase-lag and its derivatives up to five, we have developed a new eight-step singularly P-stable multiderivative method of algebraic order ten. It presents improved characteristics over other multistep multiderivative methods in the literature: reduced local truncation error, reduced phase-lag error, increased interval of periodicity, increased periodicity region, and reduced CPU time can be seen in the numerical results. The conclusion of the research presented in this paper is that the new obtained method produces more efficient results than do known ones for the numerical solution of second-order initial value problems.
Author Contributions: All authors planned the scheme, developed the mathematical modeling and examined the theory validation. The manuscript was written through the contribution by all authors. All authors discussed the results, reviewed the final version of the manuscript. All authors have read and agreed to the published version of the manuscript. Acknowledgments: It is our duty to express our gratitude to the esteemed reviewers who read this article patiently and carefully and helped us to complete this article with their appropriate and insightful comments.