# Two Embedded Pairs of Runge-Kutta Type Methods for Direct Solution of Special Fourth-Order Ordinary Differential Equations

Authors: Kasim Hussain; Fudziah Ismail; Norazak Senu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We present two pairs of embedded Runge-Kutta type methods for direct solution of fourth-order ordinary differential equations (ODEs) of the formy(iv)=f(x,y) denoted as RKFD methods. The first pair, which we will call RKFD5(4), has orders 5 and 4, and the second one has orders 6 and 5 and we will call it RKFD6(5). The techniques used in the derivation of the methods are that the higher order methods are very precise and the lower order methods give the best error estimate. Based on these pairs, we have developed variable step codes and we have used them to solve a set of special fourth-order problems. Numerical results show the robustness and the efficiency of the new RKFD pairs as compared with the well-known embedded Runge-Kutta pairs in the scientific literature after reducing the problems into a system of first-order ordinary differential equations (ODEs) and solving them.


---

## Body

## 1. Introduction

This paper deals with embedded RKFD methods for directly solving special fourth-order ordinary differential equations (ODEs) of the form(1)yivx=fx,y,x≥x0with initial conditions(2)yx0=y0,y′x0=y0′,y′′x0=y0′′,y′′′x0=y0′′′,in which the first, second, and third derivatives do not appear explicitly. This type of problems can be found in various fields of applied science and engineering such as beam theory [1, 2], fluid dynamics [3], neural networks [4], and electric circuits [5]. Traditionally, the fourth-order ordinary differential equations are transformed to a first-order system of ordinary differential equations, so that standard numerical methods can be applied (see [6–11]). However, several researchers (see [1, 12, 13]) observed the drawback of this technique as it wastes a lot of computing time and human effort. Therefore, direct integration methods have attracted significant attention from several authors for solving higher order ODEs, because these direct methods demonstrated the features in accuracy and speed (see [14–23]). However, all the methods discussed above are multistep methods in nature. This paper primarily aims to construct a one-step method to solve special fourth-order ODEs directly; this new method is self-starting in nature.The general form of RKFD method withs-stage for solving special fourth-order ODEs (1) can be expressed as follows [24]:(3)yn+1=yn+hyn′+h22yn′′+h36yn′′′+h4∑i=1sbiki,yn+1′=yn′+hyn′′+h22yn′′′+h3∑i=1sbi′ki,yn+1′′=yn′′+hyn′′′+h2∑i=1sbi′′ki,yn+1′′′=yn′′′+h∑i=1sbi′′′ki,where(4)k1=fxn,yn,ki=fxn+cih,yn+hciyn′+h22ci2yn′′+h36ci3yn′′′+h4∑j=1saijkj,for i=2,3,…,s.The parametersbi, bi′, bi′′, bi′′′, aij, and ci of the RKFD method are to be determined for i=1,2,…,s and j=1,2,…,s and supposed to be real. The RKFD method is an explicit method if aij=0 for i≤j and is an implicit method if aij≠0 for some i such that i≤j.To determine the parameters of the RKFD method given in (3)-(4), the RKFD method expression is expanded using the Taylor series expansion. After doing some algebraic simplifications, this expansion is equated to the true solution that is given by the Taylor series expansion. The direct expansion of the truncation error is used to derive the order conditions for the RKFD method [25]. A good deal of algebraic and numerical calculations are required for the above operation which were carried out using algebra package Maple [26]. Algebraic order conditions for the RKFD method can be obtained from the direct expansion of the local truncation error.In this paper we will derive embedded Runge-Kutta pairs for direct integration of special fourth-order ODEs. Embedded pairs of RK type methods have a built-in local truncation error estimate; as a result, the step size can be controlled at virtually no extra cost, and hence an efficient variable step size code can be developed.In recent years, the construction of embedded Runge-Kutta method is an effective research area yielding continuous development to the existing codes. The present paper is primarily dedicated as an extra work in this research area. This technique involves two Runge-Kutta formulae of ordersr and v (r>v, usually r=v+1) (see [25, 27–33]). We are interested in deriving the effective embedded r(v) pairs of RKFD methods that provide a cheap error estimation for variable step size codes. They depend on the methods (c,A,b,b′,b′′,b′′′) of order r and (c,A,b^,b^′,b^′′,b^′′′) of order v. Butcher tableau of embedded RKFD pair can be written as follows:(5)
​
​
c
A
b
T
b
′
T
b
′
′
T
b
′
′
′
T
b
^
T
b
^
′
T
b
^
′
′
T
b
^
′
′
′
TThe method will compute yn+1, yn+1′, yn+1′′, and yn+1′′′ to approximate y(xn+1), y′(xn+1), y′′(xn+1), and y′′′(xn+1), where yn+1 is the computed solution and y(xn+1) is the exact solution.The remainder of this paper is organized as follows. In Section2, we present the order conditions of RKFD method as well as the basic concepts and notations which are used for embedded method. In Section 3, we present the construction of the new embedded RKFD pairs of orders 5(4) and 6(5), respectively. In Section 4, we carry out the numerical experiments to show the efficiency of the new embedded RKFD pairs when compared with the well-known Runge-Kutta pairs from the scientific literature. Conclusions of the paper are given in Section 5.


## 2. The Order Conditions of RKFD Method

The order conditions of RKFD method up to fifth order have been derived using Taylor series expansion by Hussain et al. [24]. We use the same approach to derive the order conditions up to seventh order and for convenience we will present the algebraic order conditions of the RKFD method given in [24] together with the seventh order method in this paper. Next, we give the order conditions up to order 7.The order conditions fory are as follows:Fourth order:(6)∑i=1sbi=124.Fifth order:(7)∑i=1sbici=1120.Sixth order:(8)∑i=1sbici2=1360.Seventh order:(9)∑i=1sbici3=1840.The order conditions fory′ are as follows:Third order:(10)∑i=1sbi′=16.Fourth order:(11)∑i=1sbi′ci=124.Fifth order:(12)∑i=1sbi′ci2=160.Sixth order:(13)∑i=1sbi′ci3=1120.Seventh order:(14)∑i=1sbi′ci4=1210,(15)∑i,j=1sbi′aij=15040.The order conditions fory′′ are as follows:Second order:(16)∑i=1sbi′′=12.Third order:(17)∑i=1sbi′′ci=16.Fourth order:(18)∑i=1sbi′′ci2=112.Fifth order:(19)∑i=1sbi′′ci3=120.Sixth order:(20)∑i=1sbi′′ci4=130,(21)∑i,j=1sbi′′aij=1720.Seventh order:(22)∑i=1sbi′′ci5=142,∑i=1sbi′′aijcj=15040,∑i=1sbi′′ciaij=11008.The order conditions fory′′′ are as follows:First order:(23)∑i=1sbi′′′=1.Second order:(24)∑i=1sbi′′′ci=12.Third order:(25)∑i=1sbi′′′ci2=13.Fourth order:(26)∑i=1sbi′′′ci3=14.Fifth order:(27)∑i=1sbi′′′ci4=15,(28)∑i,j=1sbi′′′aij=1120.Sixth order:(29)∑i=1sbi′′′ci5=16,(30)∑i,j=1sbi′′′aijcj=1720,∑i=1sbi′′′ciaij=1144.Seventh order:(31)∑i=1sbi′′′ci6=17,∑i=1sbi′′′cj2aij=1168,∑i=1sbi′′′aijcj2=12520,∑i=1sbi′′′ciaijcj=1840.The following strategies are utilized for developing efficient embedded pairs.(1)
The quantities ofτ(r+1)2 and τ^(v+1)2 should be as small as possible for higher and lower order RKFD method, respectively, where(32)τr+12=∑i=1n1τir+12+∑i=1n2τi′r+12+∑i=1n3τi′′r+12+∑i=1n4τi′′′r+12,τ^v+12=∑i=1n^1τ^iv+12+∑i=1n^2τ^i′v+12+∑i=1n^3τ^i′′v+12+∑i=1n^4τ^i′′′v+12,whereτi(r+1), τi′(r+1), τi′′(r+1), and τi′′′(r+1) are called the error terms for y, y′, y′′, and y′′′, respectively.(2)
The following quantities given in [34] should be as small as possible:
(i)
(33)Cv+2=τ^v+2-τv+22τ^v+12,C′v+2=τ^′v+2-τ′v+22τ^′v+12,C′′v+2=τ^′′v+2-τ′′v+22τ^′′v+12,C′′′v+2=τ^′′′v+2-τ′′′v+22τ^′′′v+12.(ii)
(34)Bv+2=τ^v+22τ^v+12,B′v+2=τ^′v+22τ^′v+12,B′′v+2=τ^′′v+22τ^′′v+12,B′′′v+2=τ^′′′v+22τ^′′′v+12,whereτ^(v+1), τ^′(v+1), τ^′′(v+1), and τ^′′′(v+1) are called the error coefficients for y^, y^′, y^′′, and y^′′′ of the embedded RKFD pairs, respectively.(3)
We defined the local error estimation at the pointtn+1 by the following formula:(35)EST=max⁡βn+1∞,βn+1′∞,βn+1′′∞,βn+1′′′∞,where(36)βn+1=y^n+1-yn+1,βn+1′=y^n+1′-yn+1′,βn+1′′=y^n+1′′-yn+1′′,βn+1′′′=y^n+1′′′-yn+1′′′,wherey, y′, y′′, and y′′′ and y^, y^′, y^′′, and y^′′′ are solutions using the higher order formula and the lower order formula, respectively.The local error estimation, EST, can be used to control the step size h by the standard formula as given in [35–38](37)hn+1=0.9hnTOLEST1/v+1,where 0.9 is a safety factor, the local error estimation at each step is represented by EST, and TOL is the maximum allowable local error which is the precision required.IfEST≤TOL, then the step is accepted and we applied the procedure of performing local extrapolation (or higher order mode) meaning that the more accurate approximation will be used to advance the integration. If EST>TOL, then the step is rejected and the step size h will be updated using formula (37).


## 3. Construction of Embedded Explicit RKFD Pairs

The construction of embedded explicit RKFD pairs will be discussed in this section. In specific, we will derive two embedded RKFD pairs of orders 5(4) and 6(5) with three and four stages per step, respectively.


### 3.1. The Derivation of Embedded RKFD5(4) Pair

This section will focus on the derivation of embedded RKFD5(4) pair with three stages. The authors in [24] derived three-stage fifth-order RKFD method and the solution is given as follows:(38)c2=35+610,c3=35-610,b1′′′=19,b2′′′=49-636,b3′′′=49+636,b1′′=19,b2′′=736-618,b3′′=736+618,b1′=118,b2′=118-648,b3′=118+648,b1=191080,b2=131080-1162160,b3=131080+1162160,a21=4059187793,a31=-1502532215,a32=1826569317.Now, based on the above solution for values of A and c, we derive a three-stage fourth-order embedded formula. Solving the algebraic conditions (6), (10)-(11), (16)–(18), and (23)–(26) simultaneously gives a solution for b^i in terms of b^1 and b^3 and the solution of b^i′ in terms of b^3′ which are given as follows:(39)b^2=-b^1-b^3+124,b^1′=112+1726+25-256b^3′,b^2′=112-1726-75-256b^3′while the values of b^i′′ and b^i′′′, i=1,2,3, are the same as the fifth-order method. Our objective now is to select the values of the free parameters b^1, b^3, and b^3′ and thus the values of C(6), B(6), C′(6), B′(6), τ^(5)2, and τ^′(5)2 are as small as possible. Using the above solution, we get C′(6)=1.200 is a constant value; however, B′(6) and τ^′(5)2 are functions in terms of b^3′, by plotting the graph of τ^′(5)2 versus b^3′ (see Figure 1). From the numerical experiment and for the accurate result, we choose b^3′=1/10 usingMinimize command in Maple Software, which gives the error norm τ^′(5)2=1.14566368×10-3 and B′(6)=0.4726195824.Figure 1
The graph ofτ′(5)2 versus b^3′.The quantitiesC(6), B(6), and τ^(5)2 are functions in terms of b^1, b^3. Choose b^1=1/100 and plot the graph of C(6), B(6), and τ^(5)2 against b^3 in the interval 0,0.08. From numerical experiment and for the optimized pair, therefore b^3=3/100 is chosen which gives B(6)=0.5887532657, C(6)=0.5887532657, and τ^(5)2=3.726445732×10-3. Consequently, we denote this pair as RKFD5(4) method and it can be written in Butcher tableau as follows:(40)
​
​
3
5
+
6
10
4059
187793
3
5
-
6
10
-
1502
532215
1826
569317
19
1080
13
1080
-
11
6
2160
13
1080
+
11
6
2160
1
18
1
18
-
6
48
1
18
+
6
48
1
9
7
36
-
6
18
7
36
+
6
18
1
9
4
9
-
6
36
4
9
+
6
36
1
100
1
600
3
100
37
300
-
47
6
1800
-
17
300
+
47
6
1800
1
10
1
9
7
36
-
6
18
7
36
+
6
18
1
9
4
9
-
6
36
4
9
+
6
36


### 3.2. The Derivation of Embedded RKFD6(5) Pair

In this section a four-stage embedded RKFD6(5) pair will be derived. For the sixth-order(r=6) method, the order conditions of RKFD method up to order six need to be solved to derive a four-stage sixth-order RKFD method, namely, RKFD6 method. To determine the parameters of RKFD6 method, we choose (10)–(13) from order conditions for y′, (16)–(20) from order conditions for y′′, and (23)–(27) and (29) from order conditions for y′′′. It involved 15 nonlinear equations with 15 unknowns and then was solved simultaneously, which results in a unique solution as follows:(41)b1′=124,b2′=116-548,b3′=116+548,b4′=0,b1′′=112,b2′′=524-524,b3′′=524+524,b4′′=0,b1′′′=112,b2′′′=512,b3′′′=512,b4′′′=112,c2=12+510,c3=12-510,c4=1.Substituting the values of the above c2, c3, and c4 into the order conditions for y we get(42)b2=172-75360+5b1,b3=172+75360-5b1,b4=-b1+172.Here we have one free parameter b1 which can be chosen by minimizing the error norm of the seventh order conditions for y according to Dormand et al. [34]. The error norms and the global error of the seventh order conditions are defined as follows:(43)τ72=∑i=1np+1τi72,τ′72=∑i=1np′+1τi′72,τ′′72=∑i=1np′′+1τi′′72,τ′′′72=∑i=1np′′′+1τi′′′72,τg72=∑i=1np+1τi72+∑i=1np′+1τi′72+∑i=1np′′+1τi′′72+∑i=1np′′′+1τi′′′72,where τ(7), τ′(7), τ′′(7), and τ′′′(7) are the local truncation errors norms for y, y′, y′′, and y′′′ of the RKFD6 method, respectively, and τg(7) is the global error. As a result, we find the error equations of y as(44)τ72=16300-17+1260b12.By minimizing τ(7)2 with respect to the free parameter b1, we obtain b1=17/1260 that is the optimal value and gives τ(7)2=0 and(45)b2=172-5168,b3=172+5168,b4=12520.To find the coefficients a21, a31, a32, a41, a42, and a43, substituting the values of b1, b2, b3, b4, b1′, b2′, b3′, b4′, b1′′, b2′′, b3′′, b4′′, b1′′′, b2′′′, b3′′′, and b4′′′ into (28), (30), and (31) of the order conditions of y′′′ and (15) of the order conditions of y′, the following system of nonlinear equations needs to be solved:(46)112a32+124a42+a43+5120a42-a43=1840,(47)18+524a32+140a42+a43+5120a42-a43=12520,(48)116-548a21+116+548a31+a32=15040,(49)512a21+a31+a32+112a41+a42+a43=1120,(50)524+55120a32+124+5120a42+124-5120a43=1720,(51)524+55120a21+524-55120a31+a32+112a41+a42+a43=1144.Solving (46)–(51) simultaneously, we get(52)a21=1168+1154200,a31=1420-5700,a32=1280-5840,a41=-1840-5168,a42=184-5140,a43=5168-115840.Consequently, the local truncation error norms and the global error of the seventh order conditions of RKFD6 method are calculated and given as follows:(53)τ72=0,τ′72=2.380952381×10-4,τ′′72=4.761904762×10-4,τ′′′72=4.761904762×10-4,τg72=7.142857143×10-4.Finally, all the parameters of sixth-order RKFD method with four-stage and denoted as RKFD6 are written in Butcher tableau as follows:(54)
​
​
1
2
+
5
10
1
168
+
11
5
4200
1
2
-
5
10
1
420
-
5
700
1
280
-
5
840
1
-
1
840
-
5
168
1
84
-
5
140
5
168
+
11
5
840
17
1260
1
72
-
5
168
1
72
+
5
168
1
2520
1
24
1
16
-
5
48
1
16
+
5
48
0
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12Now, depending on the above values of A and c of RKFD6 method together with the equations up to order five, that is, (6)-(7), (10)–(12), (16)–(19), and (23)–(28), solving the system simultaneously produces a solution of an embedded RKFD formula of order five which is given by(55)b^1=130-12-510b^2-12+510b^3,b^4=1120-12+510b^2-12-510b^3,b^1′=116+580-55b^3′,b^2′=18-b^3′,b^4′=-148-580+55b^3′while the values of b^i′′ and b^i′′′, i=1,2,3,4, are the same as the RKFD6 method. Now our goal is to determine the values of the free parameters b^2, b^3, and b^3′ such that the values of C(7), B(7), C′(7), B′(7), τ^′(6)2, and τ^(6)2 are as small as possible. The quantities C(7), B(7), and τ^(6)2 are functions depending on b^2 and b^3. Picking b^2=1/40 and plotting the graph of C(7), B(7), and τ^(6)2 against b^2 and from the numerical experiment, therefore we choose b^3=1/20 that is the optimal value which gives C(7)=1.507670351, B(7)=1.507670351, and the local truncation error τ^(6)2=9.444444448×10-3. The quantities C′(7), B′(7), and τ^′(6)2 are functions in terms of b^3′. Plotting the graph of these functions and for the optimized pair, we found that b^3′=1/10 is the optimal value usingMinimize command in Maple Software, which yields C′(7)=1.683161973, B′(7)=1.708476855, and τ^′(6)2=8.1256471×10-4 (see Figure 2).Figure 2
The graph ofτ′(6)2 versus b^3′.Finally, all the parameters of the embedded pair are represented in the Butcher tableau and denoted as RKFD6(5) method as follows:(56)
​
​
1
2
+
5
10
1
168
+
11
5
4200
1
2
-
5
10
1
420
-
5
700
1
280
-
5
840
1
-
1
840
-
5
168
1
84
-
5
140
5
168
+
11
5
840
17
1260
1
72
-
5
168
1
72
+
5
168
1
2520
1
24
1
16
-
5
48
1
16
+
5
48
0
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12
-
1
240
-
5
400
1
40
1
20
-
7
240
+
5
400
1
16
-
3
5
400
1
40
1
10
-
1
48
+
3
5
400
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12


## 3.1. The Derivation of Embedded RKFD5(4) Pair

This section will focus on the derivation of embedded RKFD5(4) pair with three stages. The authors in [24] derived three-stage fifth-order RKFD method and the solution is given as follows:(38)c2=35+610,c3=35-610,b1′′′=19,b2′′′=49-636,b3′′′=49+636,b1′′=19,b2′′=736-618,b3′′=736+618,b1′=118,b2′=118-648,b3′=118+648,b1=191080,b2=131080-1162160,b3=131080+1162160,a21=4059187793,a31=-1502532215,a32=1826569317.Now, based on the above solution for values of A and c, we derive a three-stage fourth-order embedded formula. Solving the algebraic conditions (6), (10)-(11), (16)–(18), and (23)–(26) simultaneously gives a solution for b^i in terms of b^1 and b^3 and the solution of b^i′ in terms of b^3′ which are given as follows:(39)b^2=-b^1-b^3+124,b^1′=112+1726+25-256b^3′,b^2′=112-1726-75-256b^3′while the values of b^i′′ and b^i′′′, i=1,2,3, are the same as the fifth-order method. Our objective now is to select the values of the free parameters b^1, b^3, and b^3′ and thus the values of C(6), B(6), C′(6), B′(6), τ^(5)2, and τ^′(5)2 are as small as possible. Using the above solution, we get C′(6)=1.200 is a constant value; however, B′(6) and τ^′(5)2 are functions in terms of b^3′, by plotting the graph of τ^′(5)2 versus b^3′ (see Figure 1). From the numerical experiment and for the accurate result, we choose b^3′=1/10 usingMinimize command in Maple Software, which gives the error norm τ^′(5)2=1.14566368×10-3 and B′(6)=0.4726195824.Figure 1
The graph ofτ′(5)2 versus b^3′.The quantitiesC(6), B(6), and τ^(5)2 are functions in terms of b^1, b^3. Choose b^1=1/100 and plot the graph of C(6), B(6), and τ^(5)2 against b^3 in the interval 0,0.08. From numerical experiment and for the optimized pair, therefore b^3=3/100 is chosen which gives B(6)=0.5887532657, C(6)=0.5887532657, and τ^(5)2=3.726445732×10-3. Consequently, we denote this pair as RKFD5(4) method and it can be written in Butcher tableau as follows:(40)
​
​
3
5
+
6
10
4059
187793
3
5
-
6
10
-
1502
532215
1826
569317
19
1080
13
1080
-
11
6
2160
13
1080
+
11
6
2160
1
18
1
18
-
6
48
1
18
+
6
48
1
9
7
36
-
6
18
7
36
+
6
18
1
9
4
9
-
6
36
4
9
+
6
36
1
100
1
600
3
100
37
300
-
47
6
1800
-
17
300
+
47
6
1800
1
10
1
9
7
36
-
6
18
7
36
+
6
18
1
9
4
9
-
6
36
4
9
+
6
36


## 3.2. The Derivation of Embedded RKFD6(5) Pair

In this section a four-stage embedded RKFD6(5) pair will be derived. For the sixth-order(r=6) method, the order conditions of RKFD method up to order six need to be solved to derive a four-stage sixth-order RKFD method, namely, RKFD6 method. To determine the parameters of RKFD6 method, we choose (10)–(13) from order conditions for y′, (16)–(20) from order conditions for y′′, and (23)–(27) and (29) from order conditions for y′′′. It involved 15 nonlinear equations with 15 unknowns and then was solved simultaneously, which results in a unique solution as follows:(41)b1′=124,b2′=116-548,b3′=116+548,b4′=0,b1′′=112,b2′′=524-524,b3′′=524+524,b4′′=0,b1′′′=112,b2′′′=512,b3′′′=512,b4′′′=112,c2=12+510,c3=12-510,c4=1.Substituting the values of the above c2, c3, and c4 into the order conditions for y we get(42)b2=172-75360+5b1,b3=172+75360-5b1,b4=-b1+172.Here we have one free parameter b1 which can be chosen by minimizing the error norm of the seventh order conditions for y according to Dormand et al. [34]. The error norms and the global error of the seventh order conditions are defined as follows:(43)τ72=∑i=1np+1τi72,τ′72=∑i=1np′+1τi′72,τ′′72=∑i=1np′′+1τi′′72,τ′′′72=∑i=1np′′′+1τi′′′72,τg72=∑i=1np+1τi72+∑i=1np′+1τi′72+∑i=1np′′+1τi′′72+∑i=1np′′′+1τi′′′72,where τ(7), τ′(7), τ′′(7), and τ′′′(7) are the local truncation errors norms for y, y′, y′′, and y′′′ of the RKFD6 method, respectively, and τg(7) is the global error. As a result, we find the error equations of y as(44)τ72=16300-17+1260b12.By minimizing τ(7)2 with respect to the free parameter b1, we obtain b1=17/1260 that is the optimal value and gives τ(7)2=0 and(45)b2=172-5168,b3=172+5168,b4=12520.To find the coefficients a21, a31, a32, a41, a42, and a43, substituting the values of b1, b2, b3, b4, b1′, b2′, b3′, b4′, b1′′, b2′′, b3′′, b4′′, b1′′′, b2′′′, b3′′′, and b4′′′ into (28), (30), and (31) of the order conditions of y′′′ and (15) of the order conditions of y′, the following system of nonlinear equations needs to be solved:(46)112a32+124a42+a43+5120a42-a43=1840,(47)18+524a32+140a42+a43+5120a42-a43=12520,(48)116-548a21+116+548a31+a32=15040,(49)512a21+a31+a32+112a41+a42+a43=1120,(50)524+55120a32+124+5120a42+124-5120a43=1720,(51)524+55120a21+524-55120a31+a32+112a41+a42+a43=1144.Solving (46)–(51) simultaneously, we get(52)a21=1168+1154200,a31=1420-5700,a32=1280-5840,a41=-1840-5168,a42=184-5140,a43=5168-115840.Consequently, the local truncation error norms and the global error of the seventh order conditions of RKFD6 method are calculated and given as follows:(53)τ72=0,τ′72=2.380952381×10-4,τ′′72=4.761904762×10-4,τ′′′72=4.761904762×10-4,τg72=7.142857143×10-4.Finally, all the parameters of sixth-order RKFD method with four-stage and denoted as RKFD6 are written in Butcher tableau as follows:(54)
​
​
1
2
+
5
10
1
168
+
11
5
4200
1
2
-
5
10
1
420
-
5
700
1
280
-
5
840
1
-
1
840
-
5
168
1
84
-
5
140
5
168
+
11
5
840
17
1260
1
72
-
5
168
1
72
+
5
168
1
2520
1
24
1
16
-
5
48
1
16
+
5
48
0
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12Now, depending on the above values of A and c of RKFD6 method together with the equations up to order five, that is, (6)-(7), (10)–(12), (16)–(19), and (23)–(28), solving the system simultaneously produces a solution of an embedded RKFD formula of order five which is given by(55)b^1=130-12-510b^2-12+510b^3,b^4=1120-12+510b^2-12-510b^3,b^1′=116+580-55b^3′,b^2′=18-b^3′,b^4′=-148-580+55b^3′while the values of b^i′′ and b^i′′′, i=1,2,3,4, are the same as the RKFD6 method. Now our goal is to determine the values of the free parameters b^2, b^3, and b^3′ such that the values of C(7), B(7), C′(7), B′(7), τ^′(6)2, and τ^(6)2 are as small as possible. The quantities C(7), B(7), and τ^(6)2 are functions depending on b^2 and b^3. Picking b^2=1/40 and plotting the graph of C(7), B(7), and τ^(6)2 against b^2 and from the numerical experiment, therefore we choose b^3=1/20 that is the optimal value which gives C(7)=1.507670351, B(7)=1.507670351, and the local truncation error τ^(6)2=9.444444448×10-3. The quantities C′(7), B′(7), and τ^′(6)2 are functions in terms of b^3′. Plotting the graph of these functions and for the optimized pair, we found that b^3′=1/10 is the optimal value usingMinimize command in Maple Software, which yields C′(7)=1.683161973, B′(7)=1.708476855, and τ^′(6)2=8.1256471×10-4 (see Figure 2).Figure 2
The graph ofτ′(6)2 versus b^3′.Finally, all the parameters of the embedded pair are represented in the Butcher tableau and denoted as RKFD6(5) method as follows:(56)
​
​
1
2
+
5
10
1
168
+
11
5
4200
1
2
-
5
10
1
420
-
5
700
1
280
-
5
840
1
-
1
840
-
5
168
1
84
-
5
140
5
168
+
11
5
840
17
1260
1
72
-
5
168
1
72
+
5
168
1
2520
1
24
1
16
-
5
48
1
16
+
5
48
0
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12
-
1
240
-
5
400
1
40
1
20
-
7
240
+
5
400
1
16
-
3
5
400
1
40
1
10
-
1
48
+
3
5
400
1
12
5
24
-
5
24
5
24
+
5
24
0
1
12
5
12
5
12
1
12


## 4. Numerical Experiments

In this section, we present some problems which involve special fourth-order ODEs of the formy(iv)(x)=f(x,y(x)) in order to test the performance of the two new embedded pairs derived in Section 3. The numerical results obtained are compared with the well-known embedded Runge-Kutta methods which are chosen from the scientific literature after transforming the same problems to a system of first-order differential equations and solving them. The methods chosen in the numerical experiments are as follows:(i)
RKFD6(5): the new embedded Runge-Kutta pair of orders 6(5) derived in Section3 in this paper.(ii)
RKFD5(4): the new embedded Runge-Kutta pair of orders 5(4) derived in Section3 in this paper.(iii)
RK6(5)V: the embedded Runge-Kutta pair of orders 6(5) derived by Verner as given in [39].(iv)
RK5(4)D: the embedded Runge-Kutta pair of orders 5(4) with FSAL property derived by Dormand and Prince [28].Problem 1.
We consider that the homogeneous linear problem is as follows:(57)yiv=-4y,y0=0,y′0=1,y′′0=2,y′′′0=2.The exact solution is given by y(x)=exsin⁡(x). The problem is integrated in the interval 0,5.Problem 2.
We consider that the homogeneous nonlinear problem is as follows:(58)yiv=-1516y7,y0=1,y′0=12,y′′0=-14,y′′′0=38.The exact solution is given by y(x)=1+x. The problem is integrated in the interval 0,10.Problem 3.
We consider that the inhomogeneous nonlinear problem is as follows:(59)yiv=y2+cos2⁡x+sin⁡x-1,y0=0,y′0=1,y′′0=0,y′′′0=-1.The exact solution is given by y(x)=sin⁡(x). The problem is integrated in the interval 0,8.Problem 4.
We consider that the linear system is as follows:(60)yiv=e3xu,y0=1,y′0=-1,y′′0=1,y′′′0=-1,ziv=16e-xy,z0=1,z′0=-2,z′′0=4,z′′′0=-8,wiv=81e-xz,w0=1,w′0=-3,w′′0=9,w′′′0=-27,uiv=256e-xw,u0=1,u′0=-4,u′′0=16,u′′′0=-64.The exact solution is given by(61)y=e-x,z=e-2x,w=e-3x,u=e-4x.The problem is integrated in the interval 0,3.In the numerical experiments, we have computed for each method and problem the maximum global error and the number of function evaluations used in the integration. In Figures3–6, we plot these values in double logarithmic scale. From Table 1 and Figures 3–6, we observed that the new RKFD6(5) and RKFD5(4) methods are more efficient to solve directly fourth-order differential equations compared with RK6(5)V and RK5(4)D methods in terms of function evaluations and the total number of steps. This is due to the fact that when using RK6(5)V and RK5(4)D methods, the fourth-order ODEs need to be transformed into a first-order system and hence the dimension will increase four times.Table 1
The comparison results of the total number of successful steps and failure steps between RKFD6(5), RKFD5(4), RK6(5)V, and RK5(4)D methods for all problems.
Problems
Methods
Successful steps
Failure steps
1
RKFD6(5)
440
4
RKFD5(4)
705
0
RK6(5)V
461
9
RK5(4)D
768
0
2
RKFD6(5)
151
0
RKFD5(4)
200
0
RK6(5)V
208
0
RK5(4)D
253
0
3
RKFD6(5)
330
4
RKFD5(4)
466
0
RK6(5)V
348
8
RK5(4)D
497
0
4
RKFD6(5)
264
5
RKFD5(4)
340
6
RK6(5)V
339
5
RK5(4)D
484
5Figure 3
The efficiency curves for Problem1.Figure 4
The efficiency curves for Problem2.Figure 5
The efficiency curves for Problem3.Figure 6
The efficiency curves for Problem4.


## 5. Conclusion

We have constructed the two pairs of embedded RKFD methods for directly solving special fourth-order ODEs of the formy(iv)(x)=f(x,y(x)) using variable step size codes in this paper. These methods are denoted as RKFD5(4) and RKFD6(5), respectively. Based on the methods, the variables step size codes are developed and used to solve special fourth-order ODEs. Numerical results show that the two new embedded RKFD pairs are more effective as compared with the existing embedded Runge-Kutta pairs of the same order in the literature. From the numerical results we conclude that the two new embedded RKFD methods are computationally more efficient in solving special fourth-order ODEs and outperformed the well-known embedded Runge-Kutta methods.



---

*Source: 196595-2015-12-17.xml*