Variable Step Exponentially Fitted Explicit Sixth-Order Hybrid Method with Four Stages for Spring-Mass and Other Oscillatory Problems

: For the numerical integration of di ﬀ erential equations with oscillatory solutions an exponentially ﬁtted explicit sixth-order hybrid method with four stages is presented. This method is implemented using variable step-size while its derivation is accomplished by imposing each stage of the formula to integrate exactly (cid:110) 1, t , t 2 , . . . , t k , exp ( ± µ t ) (cid:111) where the frequency µ is imaginary. The local error that is employed in the step-size selection procedure is approximated using an exponentially ﬁtted explicit fourth-order hybrid method. Numerical comparisons of the new and existing hybrid methods for the spring-mass and other oscillatory problems are tabulated and discussed. The results show that the variable step exponentially ﬁtted explicit sixth-order hybrid method outperforms the existing hybrid methods with variable coe ﬃ cients for solving several problems with oscillatory solutions.


Introduction
Computation of the solutions of the special second order initial value problems y (t) = f (t, y(t)), y(t 0 ) = y 0 , y (t 0 ) = y 0 in which the first derivative does not appear explicitly has spawned many numerical algorithms and approaches. These problems often emerge in engineering and applied sciences and may be directly solved using numerical methods such as Runge-Kutta-Nystrom methods, multistep methods and hybrid methods.
In general, numerical methods are divided into two categories: (1) methods with constant coefficients and (2) methods with variable coefficients. In the development of methods with constant coefficients, one has to consider the algebraic order, the phase-lag order, the dissipation order, the strategy of minimizing the error constant and the size of interval of periodicity or the interval of absolute stability of the resulting methods (e.g., [1][2][3][4][5]). Meanwhile, development of methods with variable coefficients involves the usage of various techniques such as phase-fitting, amplification-fitting, trigonometric-fitting and exponential-fitting (e.g., [6][7][8]). These techniques are used to get methods specially adapted to the certain behavior or the structure of the solution of the problem, hence the methods with variable coefficients usually depend on the frequency of the problem. For variable step-size implementation where the step-size varies throughout the integration, accuracy and the computational cost are very much depending on the formulation of the numerical methods and the step-size control algorithm. In the step-size control algorithm, a certain procedure is usually employed Symmetry 2020, 12, 387 2 of 12 to keep the local error smaller than the desired tolerance. The local error may be approximated by a pair of embedded formulas with different algebraic orders (e.g., [9,10]) or phase lag orders (e.g., [11]). In previous work, it is observed that more attention is paid to the algebraic orders and phase-lag properties than the interval of absolute stability in the development of embedded formulas for variable step-size implementation.
In this paper, we attempt to derive an exponentially fitted explicit sixth-order hybrid method by taking into consideration the strategy of maximizing the interval of absolute stability. The method is implemented using a variable step-size in which a lower order embedded formula of hybrid method is employed for the local error estimation.
Our interest is confined to the class of explicit hybrid methods proposed by Franco [12].
where h is the step-size while f n−1 and f n represent f (t n−1 , y n−1 ) and f (t n , y n ) respectively. These methods require s -1 function evaluations or stages at each step of integration. In Butcher tableau notation, these methods can be represented by Symmetry 2020, 12, x FOR PEER REVIEW 2 of 13 error may be approximated by a pair of embedded formulas with different algebraic orders (e.g., [9,10]) or phase lag orders (e.g., [11]). In previous work, it is observed that more attention is paid to the algebraic orders and phase-lag properties than the interval of absolute stability in the development of embedded formulas for variable step-size implementation.
In this paper, we attempt to derive an exponentially fitted explicit sixth-order hybrid method by taking into consideration the strategy of maximizing the interval of absolute stability. The method is implemented using a variable step-size in which a lower order embedded formula of hybrid method is employed for the local error estimation.
Our interest is confined to the class of explicit hybrid methods proposed by Franco [12].
where h is the step-size while respectively. These methods require s -1 function evaluations or stages at each step of integration. In Butcher tableau notation, these methods can be represented by Order conditions for this class of methods are as listed in [13] while the leading term of the local truncation error for a pth-order hybrid method is defined to be   where ti  T2 and (ti) = p + 2. The error constant is the quantity is the number of trees of order p + 2 [12].

Stability Analysis
The stability analysis is presented for hybrid methods with constant coefficients and for hybrid methods with coefficients depending on the multiplication of frequency and the step-size.
Consider the standard equation Theory of linear stability analysis requires that the characteristic polynomial (Equation (4)) must be Schur stable to guarantee the stability of the corresponding method. Thus, all roots of Order conditions for this class of methods are as listed in [13] while the leading term of the local truncation error for a pth-order hybrid method is defined to be where t i ∈ T 2 and ρ(t i ) = p + 2. The error constant is the quantity where n p+2 is the number of trees of order p + 2 [12].

Stability Analysis
The stability analysis is presented for hybrid methods with constant coefficients and for hybrid methods with coefficients depending on the multiplication of frequency and the step-size.
Consider the standard equation with exact solution y(t) = C 1 exp(iλt) + C 2 exp(−iλt). Let H = λh, and e = (1 1 . . . 1) T . Applying the hybrid methods defined in Equation (1) with constant coefficients to the differential Equation (2) yields the following equation The characteristic polynomial, which determines the solution of Equation (3), is Theory of linear stability analysis requires that the characteristic polynomial (Equation (4)) must be Schur stable to guarantee the stability of the corresponding method. Thus, all roots of Equation (4) must lie inside a unit circle in a complex plane. Hence, the following definitions give conditions for the interval of periodicity and the interval of absolute stability of the methods. Definition 1. For hybrid methods corresponding to Equation (4), the interval (0, H p ) is called the interval of periodicity if P(H 2 ) = 1 and S(H 2 ) < 2 for all H ∈ (0, H p ).

Definition 2.
For hybrid methods corresponding to Equation (4), the interval (0, H a ) is called the interval of absolute stability if P(H 2 ) < 1 and S( If the coefficients of hybrid methods defined in Equation (1) are functions of v = µh where µ is the frequency of the problem and h is the step-size, then the interval of absolute stability is generalized to the region of stability. Setting µ = iw, w is real to match the oscillatory behavior of the exact solution of Equation (2), and θ = wh, we arrive at the following definition.

Exponentially Fitted Sixth-Order Method
Consider coefficients of a class of four-stage explicit hybrid methods given by Table 1: Table 1. Coefficients of a class of four-stage explicit hybrid methods.
Setting c 5 = 1 and solving equations of conditions for a sixth-order hybrid method, we get The free parameter c 3 is chosen to maximize the interval of absolute stability. Applying the transformation ζ := 1+ζ 1−ζ to the characteristic polynomial (Equation (4)) yields the equation: Symmetry 2020, 12, 387

of 12
This transformation maps the interior of a unit circle |ζ| < 1 to the left half plane Re (ζ) < 0 while the boundary of the unit circle to the imaginary axis. For stability, the points (c 3 , H) should satisfy the inequalities arising from the Routh-Hurwitz criterion. A graph with a region of feasible solutions to the Routh-Hurwitz criterion is depicted. In order to extend the interval of absolute stability, the point (c 3 , H) with the biggest possible interval (0, H) apparent in the graph is chosen. Here, we choose c 3 = 3 4 . Consequently, all other coefficients are The interval of absolute stability is (0, 4.42) while the error constant E = 2.08 × 10 −3 [14]. Using these coefficients, the characteristic polynomial (Equation (4)) is a Schur polynomial if H < 4.42. The Schur polynomial is symmetric and a basis of all symmetric polynomials.
Next, we assume that the coefficients b i (i = 1, . . . , 5), a 31 , a 32 , a 42 , a 43 , a 53 , a 54 are functions of θ while other coefficients remain constants. Then, we associate each stage formula of four-stage explicit hybrid methods with the linear operator L i [y(t)] as follows:    Based on the Taylor series expansions, it is obvious that the coefficients of the exponentially fitted sixth-order hybrid method (denoted as EXH6) will be equal to their original constant values if θ = 0. The region of stability of the method is shown in Figure 1.
Based on the Taylor series expansions, it is obvious that the coefficients of the exponentially fitted sixth-order hybrid method (denoted as EXH6) will be equal to their original constant values if  = 0. The region of stability of the method is shown in Figure 1.

Exponentially Fitted Fourth-Order Method
This section discusses the derivation of the exponentially fitted fourth-order method to be used as a tool for local error estimation. Consider a class of three-stage explicit hybrid method represented by Table 2: Table 2. Coefficients of a class of three-stage explicit hybrid methods. The fourth-order method has the same c and A values as the sixth-order method with constant coefficients described in Section 3.1. Using the order conditions for a fourth-order explicit hybrid method as listed in [13], we obtain b 1 = 0, b 2 = 19 27 , b 3 = 4 27 and b 4 = b 3 [14]. Assume that b i are functions of θ. Setting the following linear operator to zero for y(t) = t 2 , y(t) = t 3 and y(t) = exp(±µt), we have In the variable step-size setting, the local error estimation is calculated using the formula where y n+1 and y n+1 are solutions obtained from the exponentially fitted sixth-order hybrid method and the exponentially fitted fourth-order hybrid method respectively. If "tol" is the user-specified tolerance, then the step-size is selected by this procedure: If LTE < tol, then the step is accepted and the new step-size is unchanged If LTE ≥ tol, then the step is rejected and the new step-size is computed using the formula

Results
In our numerical experiments, the new and existing methods were coded using Microsoft Visual C++ software and abbreviated as follows.
EXH6: The exponentially fitted explicit sixth-order hybrid method with four stages derived in this paper. This method was implemented in variable step.
TSI6: The exponentially fitted explicit sixth-order hybrid method with four stages proposed by Tsitouras and Simos [15]. This method was implemented using the same step-size selection procedure as that for EXH6. The rational approximations of the coefficients in [15] were implemented to control the local error.
EEHM6 (4): The exponentially fitted explicit sixth-order hybrid method proposed in [16] and implemented using the step-size selection procedure as discussed in [16].
where B = 1 500 and v = 1.01. The approximate solution computed by the Galerkin method with a precision of 10 −12 is given by For this problem, we chose w = 1.
where r(t) is the radial motion over time t, k is the stiffness of the leg spring, 0 is the rest length and g is the gravitational acceleration. It is also noted that . ϕ = ω 1+ρ denotes the angular velocity with the touch-down condition ω = − g/ 0 and the spring amplitude ρ. Assume that the initial values are given by The exact solution of this problem is We solved the problem by choosing as an example k = 11, g = 9.81, 0 = 1, m = 80 and ρ = 0.001 for 0 ≤ t ≤ 100. For all methods, the selected w was √ 9.633357907.

Discussion and Conclusions
In this paper, the development of the sixth-order explicit hybrid method with an extended interval of absolute stability is described along with the exponential-fitting technique used to obtain the exponentially-fitted sixth-order hybrid method which is denoted by EXH6. In the technique, each stage formula is imposed to exactly integrate 1, t, t 2 , . . . , t k , exp(±µt) where µ is imaginary. This method is compared with the existing methods in [15] and [16] for several oscillatory problems. The derivation of TSI6 involved nullifying the phase-lag along with the phase-lag first and second derivatives [15] while the derivation of EEHM6(4) involved the usage of the exponential-fitting technique to the sixth-order hybrid method derived in [12] having an interval of periodicity (0, 2.75).
Of all methods, EXH6 gives high accuracy with a reasonable number of function evaluations for all tolerances considered in Problem 1 (the perturbed system) and has the smallest number of function evaluations for solving Problem 2 (the linear oscillatory problem). For Problem 3, EEHM6(4) has the smallest number of function evaluations followed by EXH6, while EEHM6(4) and EXH6 both perform well for Problem 4 (the nonlinear oscillatory problem). Application of the methods to the spring-mass model for running in Problem 5 shows that EXH6 can reach up to fifteenth order of accuracy while EEHM6(4) and TSI6 can only reach up to the third order and eleventh order of accuracy even for stringent tolerances. For tolerances 10 −10 and 10 −12 , the accuracy of EXH6 is slightly bigger than that for tolerance 10 −8 due to round off errors. Furthermore, TSI6 has the biggest number of function evaluations compared to the EXH6 and EEHM6(4) methods for most of the tolerances in all problems considered. This fact indicates that TSI6 needs more computational cost than that needed by EXH6 and EEHM6(4). In conclusion, the new method outperforms the existing hybrid methods for solving several special second order problems. Hence, these results offer evidence that variable step methods with variable coefficients are more efficient if the derivation of the methods takes into account the strategy of maximizing the interval of stability.