High Order Two-Derivative Runge-Kutta Methods with Optimized Dispersion and Dissipation Error

: In this work we consider explicit Two-derivative Runge-Kutta methods of a speciﬁc type where the function f is evaluated only once at each step. New 7th order methods are presented with minimized dispersion and dissipation error. These are two methods with constant coefﬁcients with 5 and 6 stages. Also, a modiﬁed phase-ﬁtted, ampliﬁcation-ﬁtted method with frequency dependent coefﬁcients and 5 stages is constructed based on the 7th order method of Chan and Tsai. The new methods are applied to 4 well known oscillatory problems and their performance is compared with the methods in that of Chan and Tsai.The numerical experiments show the efﬁciency of the derived methods.


Introduction
We consider systems of first order ODEs of the form y (x) = f (x, y(x)), x ∈ [x 0 , X], y(x 0 ) = y 0 (1) whose solutions have oscillatory or periodic behaviour. The well known Runge-Kutta (RK) methods are used to solve this problem for more than a century. Chan and Tsai [1] considered methods where only the first derivative of f (x, y(x)) (or second derivative of y) is involved; these methods are called Two-Derivative Runge-Kutta (TDRK) methods. The idea dates back to 1972 [2] where Kastlunger and Wanner introduced Runge-Kutta processes with multiple nodes. These methods use evaluations of f (x, y(x)) and its derivatives of f (x, y(x)) at the intermediate points. The advantage of TDRK methods is that they attain higher order with fewer stages than RK methods. Several authors considered numerical integration of ODEs using methods with properties such as high dispersive or dissipative order, phase-fitted, amplification-fitted, trigonometrically fitted methods, some recent articles are . These methods have either constant or variable (frequency dependent) coefficients and are suitable for the numerical integration of systems with periodic or oscillatory behavior of the solution.
Recently TDRK methods with special properties have been considered in order to integrate ODEs with periodic or oscillatory behavior of the solution. These are methods with either constant or variable (frequency dependent) coefficients.
Fang et al. [24] present trigonometrically-fitting conditions for TDRK methods of the special type where the function f is evaluated only once at each step and constructed Trigonometrically-Fitted (TFTDRK) methods. Considering the general case of TDRK methods that use several f evaluations at each step the authors in [25] derived trigonometricallyfitting conditions generalizing the conditions given in [24]. Also considering the general case the authors in [26] present 5th order TDRK methods with constant coefficients and in [27] constructed modified TFTDRK methods.
Ahmad et al. [28] presented a phase-fitted and amplification-fitted TDRK method of higher order 6 with 4 stages. Fang et al. in [29] constructed two 6th order TDRK methods with 4 stages and increased phase-lag and dissipation order for the Schrödinger equation. The authors in [30] present 5th order TDRK methods with 3 and 4 stages with minimum dispersion and dissipation error. Also in [31] they present TDRK methods of order 6 with 4 stages.
In this work we construct special explicit TDRK methods of algebraic order 6 and 7, dispersion order 8 to 12, dissipation order 7 and 9 using 4, 5 and 6 stages, as well as, a phase-fitted and amplification fitted method. In Section 2 we summarize some results on TDRK methods, Section 3 is devoted to the construction of the new methods. Numerical results are presented in Section 4 using four well known test problems.

Two-Derivative Runge-Kutta Methods
A TDRK method with s stages is of the form where y (x, y(x)) = g(x, y(x)) = f x (x, y(x)) + f y (x, y(x)) f (x, y(x)), and h is the step size. The associated Butcher tableau is c AÂ b TbT where A andÂ are square matrices of size s, b andb, c are vectors of length s.

Special TDRK Methods
Chan and Tsai [1] considered explicit methods of a specific type where the function f is evaluated only once at each step the associated Butcher tableau is This is a method of the form (2) with A and b a i1 = c i , a ij = 0 for j = 2, . . . , s, i = 1, . . . , s b 1 = 1, b i = 0 for i = 2, . . . , s Let C be a diagonal matrix with diagonal the vector c and e = (1, 1, . . . , 1) a vector of size s. The order conditions for order up to seven are given in [1] based on Butcher's the algebraic theory of trees [32] (see also Hairer [33]). The order conditions for order 2, 3 and 4 areb the order conditions for orders 5, 6 and 7 are .
As in the case of RK methods simplifying assumptions for TDRK methods are The first is the row condition A e = C e, the second isÂ e = (C 2 e)/2 and the third iŝ A C e = (C 3 e)/6.

Stability, Dispersion and Dissipation
The stability function of a TDRK method is We follow the theoretical framework of Van der Houwen and Sommeijer [34] in order to analyze the behavior of numerical methods for oscillatory problems. In [34] the definitions of dispersion (or phase-lag) error φ(v) and the dissipation (or amplification) error α(v) are given via the stability on the imaginary axis.
where v = wh and w is the frequency of the specific problem. The first is the error in the phase of the numerical solution and the second is the error in the numerical damping.
A TDRK method is said to have dispersion order p and dissipation order q if For explicit TDRK methods the real and the imaginary part of the stability function R(iv) are polynomials in v 2 of degree s

Methods with 4 Stages
Chan and Tsai [1] have shown that 6th order can be attained, they present 3 such methods. The authors in [31] presented two 6th order methods with 4 stages and optimized dispersion and dissipation order. Here we present again the construction of these methods and we shall use them later for the numerical experiments.
The coefficients c 3 and c 4 are used to increase the dispersion and the dissipation order. The following methods are derived, the first has phase-lag order 10 and dissipation order 7, this method has been also given in [29] as a 5th order method. The second method has phase-lag order 8 and dissipation order 9. The phase-lag and the dissipation (8) of these method are The stability functions (7) are The methods presented in [1] are of 6th phase-lag order and 7th dissipation order:

Method with Constant Coefficients
Chan and Tsai in [1] presented two 7th order methods with 5 stages. In this section we also consider methods with 5 stages, and follow the next procedure to derive a method of order 7.

2.
Using third simplifying assumption (6) (except for the second component) we eliminateâ i,i−2 for i = 3, 4, 5. We have to setb 2 = 0. Under this assumption the second condition for orders 5, 6 and 7 follow from the first condition.

Method with Variable Coefficients
We also derive a phase-fitted and amplification-fitted method by modifying the first 7th order 5 stages method in [1]. Asking for φ(v) = 0 and α(v) = 0 we modify the coefficientsb 1 andb 5 as followŝ The We see that the modified method reduces to the original method as v → 0.

Methods with 6 Stages
As we have seen in the construction of the 7th order method using 5 stages there is only one free parameter c 2 which we determine by nullifying the next term of the amplification error. In order to have more flexibility for further optimization of the method we consider 6 stages. We impose the first two simplifying assumptions and follow the procedure:
Solve the remaining conditions forâ i,j , i = 3, 4, 5, j = 1, . . . , i − 2 andâ 64 . The phase-lag order of this method is 8 and the dissipation order is 7. We find the remaining coefficientsâ 61 ,â 62 ,â 63 requiring the method to have phase-lag and dissipation order 12 and 9. The resulting method is:

Numerical Results
The methods used to illustrate the efficiency of the new method are The harmonic oscillator The exact solution is In Figure 1 we present the maximum absolute error of the solution in the interval [0, 1000]. For this problem we use w = 1. The fitted method has maximum absolute error less than 10 −14 and for this reason it is not in the figure. The 3 methods with minimized phase-lag and amplification error have clearly superior performance compared with the methods in [1].
The exact solution is y(x) = cos (10x) + sin (10x) + sin (x). For this problem we use w = 10. In Figure 2 we present the maximum absolute error of the solution in the interval [0, 100]. Results again show the high accuracy of the fitted method, performance of all methods tested is similar to that for Problem 1.
For this problem we use w = 5. In Figure 3 we present the maximum absolute error of the solution in the interval [0, 100], we see that the new methods produce less absolute error.

Problem 4
The almost periodic orbit problem studied by Stiefel and Bettis [38]: The exact solution is For this problem we use w = 1. In Figure 4 we present the maximum absolute error of the solution in the interval [0, 1000].

Conclusions
In this work two new methods are presented of order 7 with 5 and 6 stages and minimized dispersion and dissipation error, as well as a modified phase-fitted, amplificationfitted method with frequency dependent coefficients and 5 stages. The new methods have been tested on 4 oscillatory problems. The numerical results show that the new methods attain higher accuracy at the same cost.