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

Authors: I. B. Aiguobasimwin; R. I. Okuonghae
Journal: Journal of Applied Mathematics
 (2019)

---

## Abstract

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.


---

## Body

## 1. Introduction

In this paper we propose a class of the TDTSRK methods for the numerical solution of the IVPs in ODEs,(1)y′=fy,x∈x0,X,yx0=y0,where g=y′′=fyf and f:RN→RN and g:RN→RN. Examples of problems leading to ODEs (1) are in [1–18]. The new TDTSRK is(2)Yin=1-uiyn-1+uiyn-2+h∑j=1i-1aijfYjn+h∑j=1sbijfYjn-1+h2∑j=1sb¯ijgYjn-1,i=11s,yn=1-θyn-1+θyn-2+h∑j=1svjfYjn+wjfYjn-1+h2∑j=1sv-jgYjn+w-jgYjn-1,θ∈-1,1,n=2,3,…,N.The h is the step size and c=[c1,⋯cs]T is the abscissa value. The abscissa value {ci}i=1s represents the positions of the the internal stages. The y′′(x)=g=fyf is the second derivative form of ODEs (1). The yn is an approximation to the exact solution y(xn). The past and present stage values are Yi[n-1]=y(xn-2+cih) and Yi[n]=y(xn-1+cih) respectively, while the first and second derivatives of the past stage values are f(Yj[n-1]), and g(Yj[n-1]) respectively. The present stage value and its derivatives are(3)Yn=Y1n⋮Ysn,FYn=fY1n⋮fYsn,GYn=gY1n⋮gYsn.The compact form of (2) is(4)Yn=e-u⊗Iyn-1+u⊗Iyn-2+hA⊗IFYn+B⊗IFYn-1+h2B¯⊗IGYn-1,yn=1-θyn-1+θyn-2+hV⊗IFYn+W⊗IFYn-1+h2V¯⊗IGYn+W¯⊗IGYn-1,where e=[1,1…,1]T∈Rs and c=[c1,…,cs]T, A=[aij], U=[uij], A¯=0, B=[bij], B¯=[b¯ij], V=[vj], W=[wj], V¯T=[v¯j], W¯=[w¯j], and θ is a scalar and 0 is a null matrix. The Butcher tableau of the method in (2) is(5)
c
u
A
A
¯
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
c
1
u
1
0
0
0
⋯
0
0
⋯
0
b
11
⋯
b
1
s
b
¯
11
⋯
b
¯
1
s
c
2
u
2
a
21
0
0
⋯
0
0
⋯
0
b
21
⋯
b
2
s
b
¯
21
⋯
b
¯
2
s
c
3
u
3
a
31
a
32
0
⋯
0
0
⋯
0
b
31
⋯
b
3
s
b
¯
31
⋯
b
¯
3
s
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
⋯
⋮
⋮
⋯
⋮
⋮
⋯
⋮
c
s
u
s
a
s
1
a
s
2
a
s
3
⋯
0
0
⋯
0
b
s
1
⋯
b
s
s
b
¯
s
1
⋯
b
¯
s
s
θ
v
1
v
2
v
3
⋯
v
s
v
¯
1
⋯
v
¯
s
w
1
⋯
w
s
w
¯
1
⋯
w
¯
s
.The A is strictly lower trigular matrix while A¯ 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.


## 2. 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 p 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 q and stage order p are hd(xn)=∑j∞‍C^jy(j)(xn-1)hj and hd~(xn)=∑j∞‍Cjy(j)(xn-1)hj with(6)C=A+Be+B¯e-u;j=1,Cj=cjj!--1jj!u-Acj-1j-1!-Bc-ej-1j-1!-B¯c-ej-2j-2!;j=21q,VT+WTe+V¯T+W¯Te=1+θ;j=1,C^j=1j--1jj!θ-VTcj-1j-1!-V¯Tcj-2j-2!-WTc-ej-1j-1!-W¯Tc-ej-2j-2!;j=21p.We investigate the stability properties of the TDTSRK methods using the standard test equation:(7)y′=λy,where λ is a complex parameter. Casting the method in (2) in general linear method (GLM) format gives7
b
Y
n
y
n
y
n
-
1
h
f
Y
n
h
2
g
Y
n
=
0
A
¯
e
-
u
u
B
B
¯
v
T
v
¯
T
1
-
ϑ
ϑ
w
T
w
¯
T
0
0
1
0
0
0
I
0
0
0
0
0
0
I
0
0
0
0
h
f
Y
n
h
2
g
Y
n
y
n
-
1
y
n
-
2
h
f
Y
n
-
1
h
2
g
Y
n
-
1
.Applying the GLM in 7b to (7) yields(8)Mz=P+zV+z2V¯I-z2A¯-1Θ,where M(z) is the stability matrix, Θ=e-uuBB¯ and(9)P=1-ϑϑwTw¯T100000000000;V=vT0I0;V¯=v¯T00I.The stability polynomial, p(w,z) of the GLM in 7b is(10)pw,z=det⁡wI-Mz,where I is the identity matrix of the same order with matrix M(z).


## 3. The TDTSRK Methods with Large Region of Absolute Stability

In this section, we consider methods with stages=1, s=2 and s=3 of order two, four, and six, respectively. In (2), A is strictly lower triangular matrix while A¯ is null matrix.


### 3.1. TDTSRK Method of Orderp=2, q=2, s=1

In the spirit of [11, 19] we obtain the coefficients of TDTSRK method of order p=2 as(11)b11=c1+u1,b¯11=122c1-c12+u1-2c1u1,w1=1-v1,w¯1=123-2c1-2v¯1-2v1.The tableau of the method of order p=2 is(12)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
1
25
26
0
0
51
26
1
52
0
-
14
40
-
2
25
27
20
127
140
.The stability polynomial of this method is(13)pw,z=w4-w3-p3zw3+p2zw2+p1zw,where(14)p3z=-1013z520-31z21820,p2z=146z65+11z2280,p1z=-135z104-635z2728.The interval of absolute stability is approximately equal to [-5.53,0].


### 3.2. TDTSRK Method of Orderp=4, q=2, s=2

Fixingp=4 and solving the arising system of order conditions and stage order conditions Cj=0,j=1(1)3 and C^j=0,j=1,2 with respect to b11,b21,b11,b¯21,w1,w2 and w¯1 in (6) we obtain a twelve-parameter family of methods of order p=4 and stage order q=2. The coefficients of the methods are given by(15)b11=-b12+u1,b¯11=-b12-b¯12+u12,b21=1-a21-b22+u2,b¯21=32-a21-b22-b¯22+u22,w1=32-12v¯2-5v2,w2=-12+12v¯2-v1+4v2,w¯1=712-5v¯2-2v2,w¯2=1712-v¯1-8v¯2-4v2.The coefficients of the resulting TDTSRK method are(16)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
0
1
9
0
0
0
0
-
52
99
7
11
-
707
990
2
15
1
1
15
17
26
0
0
0
601
4290
3
11
-
371
8580
13
20
0
1
9
1
6
-
1
10
-
1
8
13
6
-
13
9
7
8
37
20The stability polynomial of this method takes the form(17)pw,z=w6-w5-p5zw5+p4zw4+p3zw3+p2zw2where(18)p5z=-8z4455-4696z219305-17z31170,p4z=-7607z8910-3890911z2231660-1043z3720-24877z458500,p3z=-13z90+18021z2772200+3914887z34633200-528301z41716000,p2z=139z22970+21977z3231660-122381z45148000.The interval of absolute stability of the fourth-order TDTSRK method is approximately equal to [-14.68,0].


### 3.3. TDTSRK Method of Orderp=6, q=4, s=3

We conclude this section with the construction of TDTSRK methods of orderp=6 and stage order q=4. Solving the system of order and stage order condition Cj, where j=1(1)6 and C^j=0 where j=1(1)4 with respect to b11,b¯11,b12,b¯12, b21,b¯21,b22,b¯22,b31,b¯31, b32,b¯32,w1,w2,w3,w¯1,w¯2 and w¯3, we obtain a twenty-one-parameter family of methods of order p=6 and stage order q=4; see appendix for the coefficients of TDTSRK methods of order p=6. The Butcher’s tableau is(19)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
0
1
15
0
0
0
0
0
0
14
165
-
36
55
7
11
41
990
-
91
495
-
2
15
1
2
1
25
2
13
0
0
0
0
0
-
237
325
31
156
11
12
-
4613
46800
-
25151
46800
1
9
1
1
20
1
13
1
17
0
0
0
0
13853
2210
-
6966
1105
19
20
43187
26520
14587
13260
5
6
0
1
8
1
2
1
6
-
1
3
5
9
-
1
11
731
165
-
863
110
4783
1320
809
990
7
6
-
503
330The stability polynomial of this method is(20)pw,z=w8-w7-p7zw7+p6zw6+p5zw5+p4zw4+p3zw3,where(21)p7z=-56873z28600-97561z2171600-71011z31640925+24z412155,p6z=5926769z5148000+730039487z2262548000-852286207z33850704000-632475763z43057912000+12819901z5147683250,p5z=-64433z396000-25121949287z219253520000-89071610501z39241689600+2262450564601z4693126720000+34222229413z5519845040000-427138145153z61559535120000,p4z=-33582133259z257760560000+1662083701181z3231042240000+1341594268543z4259922520000+301196593151z5319904640000+260210961613z61039690080000,p3z=-1729648873z328600+87855429059z4173281680000+170828953577z51039690080000+43438448533z6779767560000.The interval of absolute stability of the sixth order TDTSRK method is approximately equal to [-30.24,0]. See Figure 1.Figure 1
The region of absolute stability (closed curve) of TDTSRK methods in (7) for order p=2, p=4, and p=6.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].


## 3.1. TDTSRK Method of Orderp=2, q=2, s=1

In the spirit of [11, 19] we obtain the coefficients of TDTSRK method of order p=2 as(11)b11=c1+u1,b¯11=122c1-c12+u1-2c1u1,w1=1-v1,w¯1=123-2c1-2v¯1-2v1.The tableau of the method of order p=2 is(12)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
1
25
26
0
0
51
26
1
52
0
-
14
40
-
2
25
27
20
127
140
.The stability polynomial of this method is(13)pw,z=w4-w3-p3zw3+p2zw2+p1zw,where(14)p3z=-1013z520-31z21820,p2z=146z65+11z2280,p1z=-135z104-635z2728.The interval of absolute stability is approximately equal to [-5.53,0].


## 3.2. TDTSRK Method of Orderp=4, q=2, s=2

Fixingp=4 and solving the arising system of order conditions and stage order conditions Cj=0,j=1(1)3 and C^j=0,j=1,2 with respect to b11,b21,b11,b¯21,w1,w2 and w¯1 in (6) we obtain a twelve-parameter family of methods of order p=4 and stage order q=2. The coefficients of the methods are given by(15)b11=-b12+u1,b¯11=-b12-b¯12+u12,b21=1-a21-b22+u2,b¯21=32-a21-b22-b¯22+u22,w1=32-12v¯2-5v2,w2=-12+12v¯2-v1+4v2,w¯1=712-5v¯2-2v2,w¯2=1712-v¯1-8v¯2-4v2.The coefficients of the resulting TDTSRK method are(16)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
0
1
9
0
0
0
0
-
52
99
7
11
-
707
990
2
15
1
1
15
17
26
0
0
0
601
4290
3
11
-
371
8580
13
20
0
1
9
1
6
-
1
10
-
1
8
13
6
-
13
9
7
8
37
20The stability polynomial of this method takes the form(17)pw,z=w6-w5-p5zw5+p4zw4+p3zw3+p2zw2where(18)p5z=-8z4455-4696z219305-17z31170,p4z=-7607z8910-3890911z2231660-1043z3720-24877z458500,p3z=-13z90+18021z2772200+3914887z34633200-528301z41716000,p2z=139z22970+21977z3231660-122381z45148000.The interval of absolute stability of the fourth-order TDTSRK method is approximately equal to [-14.68,0].


## 3.3. TDTSRK Method of Orderp=6, q=4, s=3

We conclude this section with the construction of TDTSRK methods of orderp=6 and stage order q=4. Solving the system of order and stage order condition Cj, where j=1(1)6 and C^j=0 where j=1(1)4 with respect to b11,b¯11,b12,b¯12, b21,b¯21,b22,b¯22,b31,b¯31, b32,b¯32,w1,w2,w3,w¯1,w¯2 and w¯3, we obtain a twenty-one-parameter family of methods of order p=6 and stage order q=4; see appendix for the coefficients of TDTSRK methods of order p=6. The Butcher’s tableau is(19)
c
u
A
0
B
B
¯
θ
v
T
v
¯
T
w
T
w
¯
T
=
0
1
15
0
0
0
0
0
0
14
165
-
36
55
7
11
41
990
-
91
495
-
2
15
1
2
1
25
2
13
0
0
0
0
0
-
237
325
31
156
11
12
-
4613
46800
-
25151
46800
1
9
1
1
20
1
13
1
17
0
0
0
0
13853
2210
-
6966
1105
19
20
43187
26520
14587
13260
5
6
0
1
8
1
2
1
6
-
1
3
5
9
-
1
11
731
165
-
863
110
4783
1320
809
990
7
6
-
503
330The stability polynomial of this method is(20)pw,z=w8-w7-p7zw7+p6zw6+p5zw5+p4zw4+p3zw3,where(21)p7z=-56873z28600-97561z2171600-71011z31640925+24z412155,p6z=5926769z5148000+730039487z2262548000-852286207z33850704000-632475763z43057912000+12819901z5147683250,p5z=-64433z396000-25121949287z219253520000-89071610501z39241689600+2262450564601z4693126720000+34222229413z5519845040000-427138145153z61559535120000,p4z=-33582133259z257760560000+1662083701181z3231042240000+1341594268543z4259922520000+301196593151z5319904640000+260210961613z61039690080000,p3z=-1729648873z328600+87855429059z4173281680000+170828953577z51039690080000+43438448533z6779767560000.The interval of absolute stability of the sixth order TDTSRK method is approximately equal to [-30.24,0]. See Figure 1.Figure 1
The region of absolute stability (closed curve) of TDTSRK methods in (7) for order p=2, p=4, and p=6.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].


## 4. 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 valuey0 and compute the value of y1 from the exact solution or one-step 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 1.
The electrical circuit problem [21], (22)y1′x=-4y1+3y2+6,y10=0,y1x=-3.375e-2x+1.875e-0.4x+1.5,y2′x=-2.4y1+1.6y2+3.6,y20=0,y2x=-2.25e-2x+2.25e-0.4x,x∈0,1.Problem 2.
Nonlinear chemical problem(23)y1′x=-y1,y10=1,y2′x=y1-y22,y20=0,y3′x=y22,y30=0,x=0,5.Problem 3.
The Kaps problem (see [5, 18]) (24)y1′x=-y11+y1+y2,y10=1,y1x=e-x,y2′x=τy12-y2-2y2,y20=1,y2x=e-2x,x=0,1.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 (25)y¯n=yn-1+h111fY1n+122fY2n-1366fY1n-1+3533fY2n-1+h212gY1n-833gY2n-1,c=0,1T,p=3.The error estimators for order p=2 and p=6 in this paper are given in [20]. In Figures 2(a), 2(b), 3(a), 3(b), 4(a), and 4(b), the notations used are as follows:Figure 2
Efficiency Plots for Problem1.
(a) (b)Figure 3
Efficiency Plots for Problem2.
(a) (b)Figure 4
Efficiency Plots for Problem3, τ=100.
(a) (b)TOL: Tolerance used, nfe: number of function evaluations, CPU time: is the computational time and is measured in seconds,y ¯ n - y n ∞: Is the maximum error and is obtained from the difference between error estimate formula, y¯n and the output method, yn of TDTSRK method in Section 3.2,y e x a c t - y n ∞: is the maximum global error and is the difference between the exact and the numerical solution, where yexact and yn are the exact and the numerical solutions, respectively.The tests compare five different fourth-order methods with three different tolerances,Tol=10-j,j=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) and 4(a) show the logarithm of the number of function evaluations against the maximum global error (log10yexact-yn∞), while Figures 2(b) and 4(b) show the logarithm of the maximum global error (log10yexact-yn∞) against the computational effort measured by the CPU time required by each method for Problem 1 and Problem 3.The plots in Figures2(a) and 2(b) show that the fourth-order 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 Tol becomes smaller the accuracy becomes better.Similarly, Figure3 depicts the efficiency curves of the tested fourth-order methods on Problem 2. Figure 3(a) shows the logarithm of number of function evaluations against the maximum error (log10y¯n-yn∞), while Figure 3(b) depicts the logarithm of the maximum error (log10y¯n-yn∞) against the computational effort measured by the CPU time required by each method for Problem 2 for Tol=10-j,j=1,2,3. The plots in Figures 3(a) and 3(b) indicate that our fourth-order TDTSRK method compared favorably with the fourth-order 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, Figure4(a) shows the logarithm of number of function evaluations against the maximum norm of error (log10yexact-yn∞), while Figure 4(b) depicts the logarithm of the maximum norm of error (log10yexact-yn∞) against the computational effort measured by the CPU time required by each method for Problem 3 for Tol=10-j,j=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.


## 5. 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.



---

*Source: 2459809-2019-05-20.xml*