A Class of Two-Derivative Two-Step Runge-Kutta Methods for Non-Stiff ODEs

In this paper, a new class of two-derivative two-step Runge-Kutta (TDTSRK) methods for the numerical solution of non-stiff initial value problems (IVPs) in ordinary differential equation (ODEs) is considered. The TDTSRK methods are a special case of multi-derivative Runge-Kutta methods proposed by Kastlunger and Wanner (1972). The methods considered herein incorporate only the first and second derivatives terms of ODEs. These methods possess large interval of stability when compared with other existing methods in the literature. The experiments have been performed on standard problems, and comparisons were made with some standard explicit Runge-Kutta methods in the literature.


𝑔 (𝑌
) . ( The compact form of ( 2) is where and  is a scalar and 0 is a null matrix.
The Butcher tableau of the method in ( 2) is The  is strictly lower trigular matrix while  is a null matrix.The method in (2) is an extension of the methods in [2-4, 7, 9-13, 19], and a subclass of the methods in [3,5,8,15,17].Some of advantages of the TDTSRK methods over the classical Runge-Kutta (RK) methods have been highlighted in [5,8,15,17].In Section 2, the order conditions and the stability analysis of the TDTSRK methods are stated.
In Section 3, we derive TDTSRK methods using the stated order conditions and in Section 5, numerical results are presented.

The Order Condition and Stability Properties of the TDTSRK Methods
Butcher and Tracogna [3] have shown that the order conditions for TSRK method of order  can be tabulated conveniently using mapping of all rooted tress.In the spirit of [3], Turac and Ozi [18] obtained the order conditions of TDTSRK.Here, we use an equivalent set of order conditions, which may be solved directly using the strategies described in [19] and it follows from the results of [19] that the necessary and sufficient conditions which the method in (2) must satisfy in order to have methods of order  and stage order  are ℎ(  ) = ∑ ∞  Ĉ We investigate the stability properties of the TDTSRK methods using the standard test equation: where  is a complex parameter.Casting the method in (2) in general linear method (GLM) format gives ) . (7) Applying the GLM in (7) to (7) yields where () is the stability matrix, Θ = (e − u u B B) and The stability polynomial, (, ) of the GLM in (7) is where  is the identity matrix of the same order with matrix ().

The TDTSRK Methods with Large Region of Absolute Stability
In this section, we consider methods with stage  = 1,  = 2 and  = 3 of order two, four, and six, respectively.In (2),  is strictly lower triangular matrix while  is null matrix.
. .TDTSRK Method of Order  = 2,  = 2,  = 1.In the spirit of [11,19] we obtain the coefficients of TDTSRK method of order  = 2 as The tableau of the method of order  = 2 is The stability polynomial of this method is where 728 . ( The interval of absolute stability is approximately equal to [−5.53, 0]. . .TDTSRK Method of Order  = 4,  = 2,  = 2. Fixing  = 4 and solving the arising system of order conditions and stage order conditions   = 0,  = 1(1)3 and Ĉ = 0,  = 1, 2 with respect to  11 ,  21 ,  11 ,  21 ,  1 ,  2 and  1 in ( 6) we obtain a twelve-parameter family of methods of order  = 4 and stage order  = 2.The coefficients of the methods are given by The coefficients of the resulting TDTSRK method are The stability polynomial of this method takes the form where The interval of absolute stability of the fourth-order TDTSRK method is approximately equal to [−14.68, 0].
We conclude this section with the construction of TDTSRK methods of order  = 6 and stage order  = 4. Solving the system of order and stage order condition   , where  = 1(1)6 and Ĉ = 0 where  = 1(1)4 with respect to  11 , The stability polynomial of this method is where The interval of absolute stability of the sixth order TDTSRK method is approximately equal to [−30.24, 0].See Figure 1.
The stability plots of TDTSRK (2) in Figure 1 show that the interval of absolute stability of the method in (2) is larger than that of the TDTSRK methods in RK [4], TDRK [5], and TSRK [10,18].This serves as an advantage over other existing methods and this justifies the inclusion of the second derivative term of the ODEs (1) and the generalization of TSRK method in [10].

Numerical Experiment
In this section, we solve some initial value problems, and our results were compared with the results of other researchers in the literature.The fourth-order methods used for comparing are as follows: (i).TDRK methods in [5], (ii).SDTSRK methods in [15], (iii).ESDTSRK methods in [20], (iv).TDTSRK methods (new) in Section 3.2, and (v).TDTSRK methods in [18].
To start up the algorithms in (i)-(v) we use the initial value  0 and compute the value of  1 from the exact solution or onestep explicit method of order four.For easy implementation, the methods in (i)-(v) are designed to have the capacity of varying step size.Variable stepsize strategy proposed in [4] is applied herein.The non-stiff and mildly stiff problems solved are as follows: Problem .The electrical circuit problem [21], Problem .Nonlinear chemical problem Problem .The Kaps problem (see [5,18]) This problem is mildly stiff.We take  = 100 and the results are given in Figure 4.The third-order error estimator formula for the methods in (i)-(v) is The error estimators for order  = 2 and  = 6 in this paper are given in [20].In Figures 2(a TOL: Tolerance used, : number of function evaluations, CPU time: is the computational time and is measured in seconds, ‖  −   ‖ ∞ : Is the maximum error and is obtained from the difference between error estimate formula,   and the output method,   of TDTSRK method in Section 3.2, ‖  −   ‖ ∞ : is the maximum global error and is the difference between the exact and the numerical solution, where   and   are the exact and the numerical solutions, respectively. The tests compare five different fourth-order methods with three different tolerances,  = 10 − ,  = 1, 2, 3 for Problems 1-3 and for every tolerance the same initial step size is used.When carrying out a comparison among numerical methods, the criterion to be used is very important.So, if the exact solution is known, we have decided to employ the usual test based on computing of the maximum global error over the whole integration interval, because it gives a more significant measure of the efficiency.Figure 2 depicts the efficiency curves for the tested methods.Figures 2(a The plots in Figures 2(a) and 2(b) show that the fourthorder TDTSRK methods compared favorably with other existing methods in [5,18] but outpserform the TDTSRK methods in [15,20] for Problem 1 in terms of accuracy and computational time.We observed that as the value assigned to  becomes smaller the accuracy becomes better.
Similarly, Figure 3 depicts the efficiency curves of the tested fourth-order methods on Problem 2. Figure 3    method in [5], but outperformed the fourth-order TDTSRK methods in [15,18,20] for Problem 2 in terms of accuracy and computational time.
Finally, Figure 4(a) shows the logarithm of number of function evaluations against the maximum norm of error (log 10 ‖  −   ‖ ∞ ), while Figure 4(b) depicts the logarithm of the maximum norm of error (log 10 ‖  −   ‖ ∞ ) against the computational effort measured by the CPU time required by each method for Problem 3 for  = 10 − ,  = 1, 2, 3.The plots in Figures 4(a) and 4(b) indicate that our fourth-order TDRK method [5] outperforms the fourth-order methods in [15,18,20] for Problem 3 in terms of accuracy and computational time, but the new TDTSRK method is better than the fourth-order TDTSRK methods in [15,18,20] for Problem 3 in terms of accuracy and computational time.

Conclusion
We have developed and implemented our methods alongside other well-known schemes of the same order as the schemes we developed.Some initial value problems (IVPs) were used to test the efficiency of our methods.The accuracy of our method can be shown from the numerical examples.
) and 4(a) show the logarithm of the number of function evaluations against the maximum global error (log 10 ‖  −  ‖ ∞ ), while Figures 2(b) and 4(b) show the logarithm of the maximum global error (log 10 ‖  −   ‖ ∞ ) against the computational effort measured by the CPU time required by each method for Problem 1 and Problem 3.