Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 

92 

 This work is licensed under a Creative Commons Attribution 4.0 International License. 

 

  

 

 

The Implementations Special Third-Order Ordinary Differential 

Equations (ODE) for 5th-order 3rd-stage Diagonally Implicit Type Runge-

Kutta Method (DITRKM) 

 

 

Department of Mathematics, Facultyof Computer Sciences and Mathematics,Tikrit University, Iraq. 

                                                                                                             
 

 

Abstract 

    The derivation of 5th order diagonal implicit type Runge Kutta methods (DITRKM5) for 

solving 3rd special order ordinary differential equations (ODEs) is introduced in the present 

study. The DITRKM5 techniques are the name of the approach. This approach has three 

equivalent non-zero diagonal elements. To investigate the current study, a variety of tests for 

five various initial value problems (IVPs) with different step sizes h were implemented. Then, 

a comparison was made with the methods indicated in the other literature of the implicit RK 

techniques. The numerical techniques are elucidated as the qualification regarding the 

efficiency and number of function evaluations compared with another literature of the implicit 

RK approaches from the result of the computations. In addition, the stability polynomial for 

DITRK method is derived and analyzed.  

 

Keywords: Numerical Methods, Ordinary Differential Equations, Diagonal Implicit Type 

Runge Kutta Methods, Initial Value Problems, Stability Polynomial Analysis. 

 

1. Introduction 
Third-order ODEs are used in neural network engineering and applied sciences, the dynamics 

of fluid flow, the ship's motion, and electric circuits, among other fields [1-6]. Consider the 

numerical method for solving the special "initial value problems" (IVPs) for order three as the 

following form  

𝑦′′′(𝑥) = 𝑓(𝑥, 𝑦(𝑥))                                       (2.1) 

                           𝑦(𝑥0) = 𝛼,  𝑦
′(𝑥0) = 𝛽     and    𝑦

′′(𝑥0) = 𝛾 

Ibn Al Haitham Journal for Pure and Applied Science 

Journal homepage: http://jih.uobaghdad.edu.iq/index.php/j/index 

 

Doi: 10.30526/35.1.2803 

Article history: Received 29 August, 2021, Accepted 26, September, 2021, Published in January 2022. 

 

Firas A. Fawzi 

firasadil01@gmail.com 

 

Mustafa H. Jumaa 

mustafa.hassan551991@gmail.com 

 

https://creativecommons.org/licenses/by/4.0/
mailto:firasadil01@gmail.com
mailto:mustafa.hassan551991@gmail.com


Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

93 

 

   The implicit methods are important because they can reach high orders of accuracy at the 

equivalent number of stages, which can be represented as an advantage that leads to the more 

accurate than the explicit approaches. This manufactures it easier to exist the solution to the 

difficulties of the problems. 

So, the implicit RK techniques play an important role for denomination the physical and 

mathematical problems, like a differential algebraic equation.  

    In addition, diagonal implicit RK (DIRK) techniques are also pointed to as semi-implicit 

approaches or semi explicit RK techniques since they obtained at minimum one value does 

not zero for the lower of the triangular diagonal matrices. Therefore, to solve Eq. (2.1), two 

general strategies can be employed. The elementary way is to transfer the Eq. (2.1) into a 

problem with first-order then apply any pattern of the RK approach to it.  

      As a result, numerous implicit RK approaches, such as Ismail et al. [7] and others, have 

been developed. The second option is to use the RK Type method to directly solve Eq. (2.1). 

For second-order systems, several scholars provided an efficient implicit RK technique (see 

[8-13]). Ghawadri et al. [14], constructed a solution to the ill-posed issue for a beam with an 

elastically base using special fourth-order ODEs. Moreover, [15-17] developed a solution of 

the special 3rd order for the ODEs directly by RK technique. Finally, Senu [18] and Fawzi et 

al. [19] constructed the embedded the RK technique to solve 3rd order for the ODEs.  

   A significant objective for current research is to show how particular third-order ODEs are 

solving via the DIRECT method. Additionally, while solving eq. (2.1) numerically, the 

algebraic order of the technique used must be taken into account, as this is the most important 

factor in achieving high accuracy. 

    Section 2.2 demonstrates the basic idea of construction and derivation of the DITRK system 

for addressing Initial Value Problems (IVPs). The DITRK technique's order criteria are 

outlined in Section 2.3. Section 2.4 describes the 3rd stage 5th order (DITRKM5) methods. In 

Section 2.5, the analyses of the stability polynomial for the DITRK method are presented. In 

Section 2.6, mentions the DITRK approach with five IVPs. In Section 2.7, the validation of 

the current approach compared with those in the other literatures of the implicit RK techniques.  

2. The Methodology of DITRK Techniques 
For solving IVPs in eq. (2.1), the prevalent formula of the implicit RK approach for the 𝑚 

stage can be expressed as follows: 

𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑦𝑛
′ +

ℎ2

2
𝑦𝑛

′′ + ℎ3 ∑ 𝑑𝑖
𝑚
𝑖=1 𝑘𝑖                                  (2.2) 

𝑦𝑛+1
′ = 𝑦𝑛

′ + ℎ 𝑦𝑛
′′ + ℎ2 ∑  𝑏𝑖

𝑚
𝑖=1 𝑘𝑖                                              (2.3) 

𝑦𝑛+1
′′ = 𝑦𝑛

′′ + ℎ ∑  𝑔𝑖
𝑚
𝑖=1 𝑘𝑖                                                           (2.4) 

and 

𝑘1 = 𝑓(𝑥𝑛, 𝑦𝑛)                                                                              (2.5) 

𝑘𝑖 = 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑦𝑛 + ℎ 𝑐𝑖 𝑦𝑛
′ +

ℎ2

2
𝑐𝑖

2 𝑦𝑛
′′ + ℎ3 ∑ 𝑎𝑖𝑗

𝑖−1
𝑗=1 𝑘𝑗 )          (2.6) 

    where 𝑖 = 2,3, … , 𝑚 . 



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

94 

 

     The parameters of diagonal implicit RK type (DITRK) methods are presumed as 

 𝑐𝑖 , 𝑎𝑖𝑗 , 𝑑𝑖 ,  𝑏𝑖 ,  𝑔𝑖  where 𝑖, 𝑗 = 1, 2, 3 … , 𝑠 are real numbers and 𝑚 is referred to stage digit for 

the approach. This scheme is known as diagonal implicit when 𝑎𝑖𝑗 ≠ 0 for 𝑗 > 𝑖. The last 

denomination includes the single DITRK techniques that 𝐴 indicate that the lower the 

triangular diagonal matric of 𝐴 have same values with 𝑎𝑖𝑗 ≠ 0 where 𝑖 = 𝑗 at the diagonal. 

The DITRK approach proposed from the work of Butcher, as illustrated in Table 2.1 [20]. 

 

Table 2.1: Butcher form DITRK method. 

 

 

 

 

 

 

 

3. Order Conditions of the DITRK Technique 
     According to Mechee et al. [17], the orders of algebraic criteria for RKD approached over 

order 6 are as follow: 

Order conditions of 𝑦: 

order 3                ∑ 𝑑𝑖 =
1

6
                                                                        (2.7)                                                                       

order 4        ∑ 𝑑𝑖 𝑐𝑖 =
1

24
                                                                   (2.8) 

order 5        ∑ 𝑑𝑖 𝑐𝑖
2 =

1

60
                                                                  (2.9) 

order 6          ∑ 𝑑𝑖 𝑐𝑖
3 =

1

120
      and    ∑ 𝑑𝑖 𝑎𝑖,𝑗 = 1/720.                  (2.10) 

Order conditions of 𝒚′: 

order 2              ∑ 𝑏𝑖 =
1

2
                                                                        (2.11)                                                                   

order 3            ∑ 𝑏𝑖 𝑐𝑖 =
1

6
                                                                    (2.12) 

order 4                  ∑ 𝑏𝑖 𝑐𝑖
2 =

1

12
                                                                 (2.13) 

order 5              ∑ 𝑏𝑖 𝑐𝑖
3 =

1

20
      and  ∑ 𝑏𝑖 𝑎𝑖,𝑗 =

1

120
                            (2.14) 

order 6         ∑ 𝑏𝑖 𝑐𝑖
4 =

1

30
 , ∑ 𝑏𝑖 𝑎𝑖,𝑗 𝑐𝑗 =

1

720
  and   ∑ 𝑏𝑖 𝑐𝑖  𝑎𝑖,𝑗 =

1

180
  .        (2.15) 

Order conditions of 𝒚′′: 

order 1                  ∑ 𝑔𝑖 = 1                                                                      (2.16)   



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

95 

 

order 2                         ∑ 𝑔𝑖 𝑐𝑖 =
1

2
                                                                   (2.17)                                                           

order 3                 ∑ 𝑔𝑖 𝑐𝑖
2 =

1

3
                                                                   (2.18)    

order 4               ∑ 𝑔𝑖 𝑐𝑖
3 =

1

4
    and  ∑ 𝑔𝑖 𝑎𝑖,𝑗 =

1

24
                                  (2.19) 

order 5            ∑ 𝑔𝑖 𝑐𝑖
4 =

1

5
 ,   ∑ 𝑔𝑖 𝑎𝑖,𝑗 𝑐𝑗 =

1

120
 and   ∑ 𝑔𝑖 𝑐𝑖  𝑎𝑖,𝑗 =

1

30
             (2.20) 

order 6           ∑ 𝑔𝑖 𝑐𝑖
2𝑎𝑖,𝑗 =

1

36
 , ∑ 𝑔𝑖 𝑎𝑖,𝑗  𝑐𝑗

2 + ∑ 𝑔𝑖 𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 =
7

720
 ,  

    ∑ 𝑔𝑖 𝑐𝑖
5 =

1

6
 ,  ∑ 𝑔𝑖 𝑎𝑖,𝑗  𝑐𝑗

2 = 
1

360
  , ∑ 𝑔𝑖 𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 =

1

144
   

and     

                                    
1

2
∑ 𝑔𝑖 𝑎𝑖,𝑗  𝑐𝑗

2 + ∑ 𝑔𝑖 𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 =
1

120
                               (2.21) 

 

4. Formation of the 3rd stage 5th order (DITRKM5) Method 

    We implement a diagonal implicit type Runge–Kutta approach using order conditions 

derivations as demonstrated in section 2.3, which is developed according to Mechee work 

[17]. For the 𝑝 order DITRK approach, the local truncation error is defined as follows: 

for order five: 

𝐿𝑇𝐸order 5 = [(∑ 𝑑𝑖 𝑐𝑖
3 −

1

120
)2 + (∑ 𝑑𝑖 𝑎𝑖,𝑗 −

1

720
)2 + (∑ 𝑏𝑖 𝑐𝑖

4 −
1

30
)2 + (∑ 𝑏𝑖 𝑎𝑖,𝑗 𝑐𝑗 −

1

720
)2 +

(∑ 𝑏𝑖 𝑎𝑖,𝑗 𝑐𝑖 −
1

180
)2 + (∑ 𝑔𝑖 𝑎𝑖,𝑗 𝑐𝑖

2 −
1

36
)2 + (∑ 𝑔𝑖  𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 −

1

144
)2 + (∑ 𝑔𝑖 𝑐𝑖

5 −
1

6
)2 +

(∑ 𝑔𝑖 𝑎𝑖,𝑗 𝑐𝑗
2 −

1

360
)2 + (∑ 𝑔𝑖  𝑎𝑖,𝑗 𝑐𝑗

2 + ∑ 𝑔𝑖  𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 −
1

720
)2 + (

1

2
∑ 𝑔𝑖  𝑎𝑖,𝑗 𝑐𝑗

2 + ∑ 𝑔𝑖  𝑐𝑖 𝑎𝑖,𝑗 𝑐𝑗 −

1

120
)2] 

1

2                                                                                                              (2.22) 

the error of local truncation terms for 𝑦, 𝑦′and 𝑦′′. The fifth-order three-stage of the DITRKM 

method in the present study can be computed by employing algebraic order conditions over 5. 

The System of the result includes 16 nonlinear equations with 16 unknown variables, assuming  

𝑎1,1 = 𝑎2,2   and      𝑎2,2 = 𝑎3,3                                                     (2.23) 

Thus, the calculations of the system products the set of solutions in terms of the parameters 

𝑎1,1 , 𝑎2,2  and 𝑐1  as follows: 

𝑎2,1 = 0, 𝑎3,1 =
3

20
−

3

10
RootOf(10z2 − 10z + 1), 𝑎3,2 = 0, 𝑎3,3 =

1

120
RootOf(10z2 − 10z + 1),

𝑏1 =
−5

18
RootOf(10z2 − 10z + 1) +

5

18
, 𝑏2 =

2

9
, 𝑏3 =  

5

18
RootOf(10z2 − 10z + 1), 𝑐1 =

RootOf(10z2 − 10z + 1) , 𝑐2 =
1

2
 , 𝑐3 = −RootOf(10z

2 − 10z + 1) + 1, d1 =
1

8
−

5

36
RootOf(10z2 − 10z + 1), d2 =

1

18
, d3 =

5

36
 RootOf(10z2 − 10z + 1) −

1

72
,  𝑔1 =

5

18
, 𝑔2 =

4

9
 ,    𝑔3 =  

5

18
                                    (2.24) 



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

96 

 

 Finally, coefficients of the DITRKM method for 5-order 3 stages indicated by DITRKM5 can 

be read as shown in Table 2.2. 

 

Table 2.2: third-stage fifth-order DITRK Method (DITRKM5). 

 

 

 

 

 

 

 

 

5. The Stability Polynomial of DITRK Method 

       In order to study the stability polynomial of the DITRK method, the following equation 

is suggested : 

               𝑦′′′ = −𝛾 3𝑦                                                                                        (2.25) 

By substituting the DITRK method from eqs. (2.2) - (2.4) to exam eq. (2.25) with stage 𝑚 =

3  are yield 

𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑦𝑛
′ +

ℎ2

2
𝑦𝑛

′′ + ℎ3 ∑ 𝑑𝑖
3
𝑖=2 𝑘𝑖                                     (2.26) 

𝑦𝑛+1
′ = 𝑦𝑛

′ + ℎ 𝑦𝑛
′′ + ℎ2 ∑  𝑏𝑖

3
𝑖=2 𝑘𝑖                                                  (2.27) 

𝑦𝑛+1
′′ = 𝑦𝑛

′′ + ℎ ∑  𝑔𝑖
3
𝑖=2 𝑘𝑖                                                               (2.28) 

and 

              

𝑘1 = 𝑓(𝑥𝑛, 𝑦𝑛 )                                                                              
 

𝑘𝑖 = 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑦𝑛 + ℎ 𝑐𝑖 𝑦𝑛
′ +

ℎ2

2
𝑐𝑖

2 𝑦𝑛
′′ + ℎ3 ∑ 𝑎𝑖𝑗

𝑖−1
𝑗=1 𝑘𝑗 ) 

              (2.29)                        

where 𝑖 = 2,3, … , 𝑚 and 

            𝑌𝑖 = 𝑦𝑛 + ℎ 𝑐𝑖 𝑦𝑛
′ +

ℎ2

2
𝑐𝑖

2 𝑦𝑛
′′ + ℎ3 ∑ 𝑎𝑖𝑗

𝑖−1
𝑗=1 (−𝛾

3𝑌𝑗 )                            (2.30) 

where    

          𝑌1 = 𝑦𝑛                                                                                                    (2.31) 

         𝑌𝑖 = 𝑦𝑛 + ℎ 𝑐𝑖 𝑦𝑛
′ +

ℎ2

2
𝑐𝑖

2 𝑦𝑛
′′ + ℎ3 ∑ 𝑎𝑖𝑗

𝑖−1
𝑗=1 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑌𝑖 )                     (2.32) 

more simplification  

       𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑦𝑛
′ +

ℎ2

2
𝑦𝑛

′′ + ℎ3 ∑ 𝑑𝑖
𝑚
𝑖=2 (−𝛾

3)𝑌𝑗                                       (2.33) 

       𝑦𝑛+1
′ = 𝑦𝑛

′ + ℎ 𝑦𝑛
′′ + ℎ2 ∑  𝑏𝑖

𝑚
𝑖=2 (−𝛾

3)𝑌𝑗                                                   (2.34) 

       𝑦𝑛+1
′′ = 𝑦𝑛

′′ + ℎ ∑  𝑔𝑖
𝑚
𝑖=2 (−𝛾

3)𝑌𝑗                                                                 (2.35) 

From eq. (2.25), the above equations can be written as follow 

       𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑦𝑛
′ +

ℎ2

2
𝑦𝑛

′′ + ℎ3 ∑ 𝑑𝑖
𝑚
𝑖=2 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑌𝑖 )                           (2.36) 



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

97 

 

       𝑦𝑛+1
′ = 𝑦𝑛

′ + ℎ 𝑦𝑛
′′ + ℎ2 ∑  𝑏𝑖

𝑚
𝑖=2 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑌𝑖 )                                       (2.37) 

       𝑦𝑛+1
′′ = 𝑦𝑛

′′ + ℎ ∑  𝑔𝑖
𝑚
𝑖=2 𝑓(𝑥𝑛 + 𝑐𝑖 ℎ, 𝑌𝑖 )                                                     (2.38) 

From eq. (2.30), multiply eq. (2.34) by ℎ and eq. (2.35) by ℎ2,yield 

       ℎ𝑦𝑛+1
′ = ℎ𝑦𝑛

′ + ℎ2 𝑦𝑛
′′ + ℎ3 ∑  𝑏𝑖

𝑚
𝑖=2 (−𝛾

3)𝑌𝑗                                           (2.39) 

       ℎ2𝑦𝑛+1
′′ = ℎ2𝑦𝑛

′′ + ℎ3 ∑  𝑔𝑖
𝑚
𝑖=2 (−𝛾

3)𝑌𝑗                                                   (2.40) 

Then 

[

𝑦𝑛+1
ℎ𝑦𝑛+1

′

ℎ2𝑦𝑛+1
′′

] = [
1 1 1

2

0 1 1
0 0 1

] [

𝑦𝑛   

ℎ𝑦𝑛
′

ℎ2𝑦𝑛
′′

] + (−𝛾 3ℎ3) [
𝑑1 𝑑2 𝑑3
 𝑏1  𝑏2  𝑏3
 𝑔1  𝑔2  𝑔3

] [

𝑦𝑛   

ℎ𝑦𝑛
′

ℎ2𝑦𝑛
′′

]          (2.41) 

Also, the matrix format of eq. (2.30) can be defined as 

[

𝑌1
𝑌2
⋮

𝑌𝑚

] = [

1
1
⋮
1

0
𝑐2
⋮

𝑐𝑚

0
𝑐2

2

⋮
𝑐𝑚

2

] [

𝑦𝑛   

ℎ𝑦𝑛
′

⋮
ℎ2𝑦𝑛

′′

] + (−𝛾3ℎ3) [

0 0 ⋯ 0
 𝑎21 0 ⋯ 0

⋮
 𝑎𝑚1

⋮
 𝑎𝑚2

⋱ ⋮
…  𝑎𝑚−1

] [

𝑌1
𝑌2
⋮

𝑌𝑚

]       (2.42) 

where 𝑚 = 3. Therefore, 

[

𝑦𝑛+1
ℎ𝑦𝑛+1

′

ℎ2𝑦𝑛+1
′′

] = 𝑓(𝐻) [

𝑦𝑛   

ℎ𝑦𝑛
′

ℎ2𝑦𝑛
′′

] ,      𝐻 = (−𝛾 3ℎ3)                                                      (2.43)  

and 

𝑓(𝐻) = [

1 + 𝐻𝑑𝑇 𝑃−1𝐸1 1 + 𝐻𝑑
𝑇 𝑃−1𝐸2 1 + 𝐻𝑑

𝑇 𝑃−1𝐸3
𝐻𝑏𝑇 𝑃−1𝐸1 1 + 𝐻𝑏

𝑇 𝑃−1𝐸2 1 + 𝐻𝑏
𝑇 𝑃−1𝐸3

𝐻𝑔𝑇 𝑃−1𝐸1 𝐻𝑔
𝑇 𝑃−1𝐸2 1 + 𝐻𝑏

𝑇 𝑃−1𝐸3

]                    (2.44)    

where 

𝐸1 = [
1
1
1

] , 𝐸2 = [
0
𝑐2
𝑐3

] , 𝐸3 = [

0
𝑐2

2

𝑐𝑚
2

]     and 𝑃−1 = (1 − 𝐻𝐴)−1                         (2.45)    

𝐴 = [

0 0 ⋯ 0
 𝑎21 0 ⋯ 0

⋮
 𝑎𝑚1

⋮
 𝑎𝑚2

⋱ ⋮
…  𝑎𝑚−1

]    ,   𝐵 = [

𝑑1 𝑑2 𝑑3
 𝑏1  𝑏2  𝑏3
 𝑔1  𝑔2  𝑔3

]                            (2.46) 

Thus, the stability polynomial of the DITRK method can be written as 

    ∅(𝜑, 𝐻) = |𝜑𝐼 − 𝑓(𝐻)|                                                                             (2.47)  

Where 𝑓(𝐻) is given value, the characteristic equation is defined as follow, 

      ∅(𝜑, 𝐻) = 𝑃0(𝐻)𝜑
3 + 𝑃1(𝐻)𝜑

2 + 𝑃2(𝐻)𝜑 + 𝑃3(𝐻)                            (2.48) 

6. Test of Problems 

The approaches that demonstrated in section 2.3 tested with 5 various problems in this part. 

The numerical results of the suggested approaches compared with those of other RK 

techniques at equivalent order which are already available. The numerical experiments were 

conducted using the following methods: 

(1) DITRKM5: 3rd stage 5th order DITRK approach computed in the present work. 

(2) Radau I: 3rd stage 5th order RK technique presented in [20].  



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

98 

 

(3) Radau IA: 3rd stage 5th order RK method studied in [21].  

(4) Radau II: 3rd stage 5th order RK scheme tested in [20].  

(5) Radau IIA: 3rd stage 5th order RK approach noted in [21].  

Problem (1): Consider a nonhomogeneous linear ODE given in [23] 

𝑦′′′(𝑥) = 𝑦(𝑥) +  cos(𝑥), with  𝑦(0) = 0, 𝑦′(0) = 0,  𝑦′′(0) = 1 where  𝑥 ∈ [0,1], 

and analytic solution  𝑦(𝑥) =
(e𝑥  −cos(𝑥)−sin (𝑥))

2
. 

Problem (2): Consider the nonhomogeneous nonlinear ODE 

𝑦′′′(𝑥) = (𝑦(𝑥))2 + cos2(𝑥) − cos(𝑥) − 1, with     𝑦(0) = 0, 𝑦′(0) = 1, 𝑦′′(0) = 1 where  

0 ≤ 𝑥 ≤ 2,    the exact solution   𝑦(𝑥) = sin (𝑥). 

Problem (3): The nonhomogeneous nonlinear ODEs is considered as 

          𝑦1
′′′(𝑥) = 𝑦2(𝑥), with  𝑦1(0) = 1, 𝑦1

′ (0) = 0, 𝑦1
′′(0) = 1,   

𝑦2
′′′(𝑥) = −𝑦1(𝑥) − 2 𝑦2(𝑥) + 2 𝑦3(𝑥) with  𝑦2(0) = 0, 𝑦2

′ (0) = 1, 𝑦2
′′(0) = 0, 

         𝑦3
′′′(𝑥) = 𝑦1(𝑥) + 𝑦2(𝑥) with  𝑦3(0) = 1, 𝑦3

′ (0) = 1, 𝑦3
′′(0) = 1,With analytic 

 solution 𝑦1(𝑥) = cosh (𝑥) ,  𝑦2(𝑥) = sinh (𝑥) and  𝑦3(𝑥) = e
𝑥   where  0 ≤ 𝑥 ≤ 1. 

7. Numerical Results 

     Figure (2.1) shows the efficiency of the DITRKM methods created by charting of decimal 

logarithm for the highest "global error" versus logarithm of function estimate.  

 

 

 

 

 

 
 

                    problem 1                                                   problem 2                                             problem 3                                                               

Figure 2.1: Accuracy curve for DITRKM5, Radau I, Radau IA, Radau II and Radau IIA with h = 0.1, 0.05 

0.025, 0.00125, 0.00625 for the problem 1, problem 2 and problem 3. 

When compared the current study with another implicit RK approach for equivalent order, the 

DITRKM5 method requires fewer "function evaluations". The digit of equations increased 

three times with the problems turned to a system of 1st order ODEs. In the comparison, the 

existing implicit the RK approach with the equivalent order, the "global error" and digit of 

"function estimate" contain the smallest maximum for the DITRKM5 method at each iteration, 

as shown in Figure (2.1) that obtained from Table (2.1). As shown in Figure (2.1), the fifth-

order three stage results DITRK method (DITRKM5) produces more accurate findings than 

the other results in the literature (Radau I, Radau IA, Radau II, and Radau IIA). In this work, 

the logarithm of "maximum global error" is known as a logarithm function for "function 

1.8 2 2.2 2.4 2.6 2.8 3 3.2
-7

-6

-5

-4

-3

-2

-1

Log (Number of function call)

L
o
g
 (

M
a
x
 G

lo
b
a
l 
E

rr
o
rs

) 

 

 

DITRKM5

Radau I

Radau IA

Radau II

Radau IIA

2 2.5 3 3.5
-6

-5.5

-5

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

Log (Number of function call)

L
o
g
 (

M
a
x
 G

lo
b
a
l 
E

rr
o
rs

) 

 

 

DITRKM5

Radau I

Radau IA

Radau II

Radau IIA

1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6
-7

-6

-5

-4

-3

-2

-1

0

Log (Number of function call)

L
o
g
 (

M
a
x
 G

lo
b
a
l 
E

rr
o
rs

) 

 

 

DITRKM5

Radau I

Radau IA

Radau II

Radau IIA



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

99 

 

evaluation" with different step size ℎ = 0.1, 0.05, 0.025, 0.0125,0.00625  for five test 

problems. 

Table 2.1: Comparisons of number of function call and maximum global error for DITRKM5, Radau I, Radau 

IA, Radau II and Radau IIA Methods with ℎ = 0.1, 0.05 0.025, 0.00125, 0.00625 for the problem 1, problem 2 

and problem 3. 

Problem 3 Problem 2 Problem 1 Method 
Step size 

(h) 

Maximum 

Global 

Error 

Maximum 

Global 

Error 

Maximum 

Global 

Error 

No. of 

Function 

Call 

Maximum 

Global Error 

No. of 

Function 

Call 

  

9.74073E-05 88 0.000277124 160 9.16352E-05 88 DITRKM5 

0.1 

0.05488292 99 0.02489423 180 0.02986368 99 Radau I 

0.06335586 99 0.02867864 180 0.03449261 99 Radau IA 

0.05488292 99 0.02647426 180 0.02878467 99 Radau II 

0.06335586 99 0.0298746 180 0.03336877 99 Radau IIA 

1.8304E-05 160 6.8347E-05 320 1.74282E-05 160 DITRKM5 

0.05 

0.02293305 180 0.01292113 360 0.01273117 180 Radau I 

0.02656216 180 0.01494992 360 0.01475188 180 Radau IA 

0.02293305 180 0.01360768 360 0.01229666 180 Radau II 

0.02656216 180 0.015539 360 0.0143079 180 Radau IIA 

4.67234E-06 320 1.79521E-05 648 4.45635E-06 320 DITRKM5 

0.025 

0.0115538 360 0.006847777 729 0.006421128 360 Radau I 

0.01340439 360 0.007940234 729 0.007451075 360 Radau IA 

0.0115538 360 0.007152807 729 0.006195875 360 Radau II 

0.01340439 360 0.008219884 729 0.007223443 360 Radau IIA 

1.22981E-06 648 4.35276E-06 1288 1.1728E-06 648 DITRKM5 

0.0125 

0.005945227 729 0.003387183 1449 0.003298104 729 Radau I 

0.006903284 729 0.003931975 1449 0.003829934 729 Radau IA 

0.005945227 729 0.003535766 1449 0.003179202 729 Radau II 

0.006903284 729 0.004074293 1449 0.003710422 729 Radau IIA 

3.02667E-07 1288 2.88938E-06 2560 2.88938E-07 1288 DITRKM5 

0.00625 

 

0.002941504 1449 0.001667861 2880 0.00163419 1449 Radau I 

0.003416977 1449 0.001937206 2880 0.001898433 1449 Radau IA 

0.002941504 1449 0.001741931 2880 0.001575303 1449 Radau II 

0.003416977 1449 0.001741931 2880 0.001839395 1449 Radau IIA 

 

8. Conclusion 

     The test of the fifth-order three-stage DITRKM5 methods for the integration of ODEs 

obtained by testing the minimized error norm and studied the digit of the evaluations of 

function that described in this paper. The digit of "function evaluations" and the maximum 

error of the supposed technique is lower than those in implicit of the RK approaches, as 

indicated by numerical results in all Figure (2.1) and Table (2.1), and the introduced 



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

100 

 

techniques are most accurate by solving the directly specific ODEs of the order of three. The 

stability polynomial for DITRK method is analyzed and presented.  

 

References 

1. Wu, X.; Wang, Y.; Price, W. Multiple resonances, responses, and parametric instabilities in 

offshore structures. J. Ship Res. 1988, 32, 285–296. 

2.  Cortell, R. Application of the fourth-order Runge–Kutta method for the solution of high-

order general initial value problems. Comput. Struct. 1993, 49, 897–900. 

3. Malek, A.; Beidokhti, R. Numerical solution for high order differential equations using a 

hybrid neural network—Optimization method. Appl. Math. Comput. 2006, 183, 260–271. 

4.  Boutayeb, A.; Chetouani, A. A mini-review of numerical methods for high-order problems. 

Int. J. Comput. Math. 2007, 84, 563–579. 

5.  Alomari, A.; Anakira, N.R.; Bataineh, A.S.; Hashim, I. Approximate solution of nonlinear 

system of bvp arising in fluid flow problem. Math. Probl. Eng. 2013, 136043. 

6.  Twizell, E. A family of numerical methods for the solution of high-order general initial 

value problems. Comput. Methods Appl. Mech. Eng. 1988, 67, 15–25. 

7.  Ismail, F.; Jawias, N.I.C.; Suleiman, M.; Jaafar, A. Solving linear ordinary differential 

equations using singly diagonally implicit Runge–Kutta fifth order five-stage method. 

WSEAS Trans. Math. 2009, 8, 393–402. 

8.  Ismail, F. Diagonally implicit Runge–Kutta Nystrom general method order five for solving 

second order ivps. WSEAS Trans. Math. 2010, 9, 550–560. 

10.  Senu, N.; Suleiman, M.; Ismail, F.; Arifin, N.M. New 4 (3) pairs diagonally implicit 

Runge–Kutta-Nyström method for periodic IVPs. Discret. Dyn. Nat. Soc. 2012, 324989. 

11.  Senu, N.; Suleiman, M.; Ismail, F.; Othman, M. A singly diagonally implicit Runge–

Kutta-Nyström method for solving oscillatory problems. IAENG Int. J. Appl. Math. 2011, 41, 

155–161. 

12.  Senu, N.; Suleiman, M.; Ismail, F.; Othman, M. A fourth-order diagonally implicit Runge–

Kutta-Nyström method with dispersion of high order. In Proceedings of the 4th International 

Conference on Applied Mathematics, Simulation, Modelling (ASM’10), Corfu Island, Greece, 

22–25 July 2010; 78–82. 

13.  Attili, B.S.; Furati, K.; Syam, M.I. An efficient implicit Runge–Kutta method for second 

order systems. Appl. Math. Comput. 2006, 178, 229–238. 

14.  Wen, Z.; Zhu, T.; Xiao, A. Two classes of three-stage diagonally-implicit Runge–Kutta 

methods with an explicit stage for stiff oscillatory problems. Math. Appl. 2011, 1, 15. 

15.  Ghawadri, N., Senu, N., Adel Fawzi, F., Ismail, F., & Ibrahim, Z. B. Diagonally implicit 

Runge–Kutta type method for directly solving special fourth-order ordinary differential 

equations with Ill-Posed problem of a beam on elastic foundation. Algorithms, 2019 12(1), 

10.  



Ibn Al-Haitham Jour. for Pure & Appl. Sci. 53 (1)2022 
 

101 

 

16.  Mechee, M., Senu, N., Ismail, F., Nikouravan, B., & Siri, Z. (2013). A three-stage fifth-

order Runge-Kutta method for directly solving special third-order differential equation with 

application to thin film flow problem. Mathematical Problems in Engineering, 2013.  

17.  F A Fawzi, N Senu, F Ismail and Z A Majid A New Integrator of Runge-Kutta Type for 

Directly Solving General Third-order ODEs with Application to Thin Film Flow Problem 

Applied Mathematics & Information Sciences, 2018. 12 (4), 775-784. 

18.  F A Fawzi and N Senu and F Ismail An efficient of direct integrator of Runge-Kutta type 

method for solving y''' = f(x,y,y') with application to thin film flow problem international 

journal of pure and applied mathematics, 2018,  120, 27-50. 

19.  N Senu, M Suleiman, F Ismail. An embedded explicit Runge-Kutta-Nyström method for 

solving oscillatory problems, Physica Scripta, 2009,  80 , 015005. 

20. Fawzi, F. A., Abdullah, Z. M., & Hussein, N. K. Two Embedded Pairs for Solve Directly 

Third-Order Ordinary Differential Equation by Using Runge-Kutta Type Method (RKTGD), 

Journal of Physics: Conference Series (Vol. 1879, No. 2, p. 022123). IOP Publishing, 2021. 

21. Butcher, John Charles, and Nicolette Goodwin. Numerical methods for ordinary 

differential equations. New York: Wiley, 2008. 2. 

22. Dormand, J. R. Numerical methods for differential equations: a computational approach 

CRC press, 1996.3. 

23. Hussain, K. A., Ismail, F., Senu, N., Rabiei, F., & Ibrahim, R. Integration for special third-

order ordinary differential equations using improved Runge-Kutta direct method, MJS, 2015 

34(2), 172-179.