Differential equations (DEs) are of great use in mode-ling different real life problems arising in science and engineering (Arora, 2019). Model equations for-med by using DEs get complicated and several times it becomes quite difficult to find its exact solution (Ahmed and Iqbal, 2020). However, to find the exact solution of a complicated model equation a practice is to simplify the model equation and then find the exact solution of the simplified equation, after then they obtained result is used to approximate the original equation (Islam, 2015). 

In this circumstance, the approximate result differs from the real one. To avoid such inconvenience researchers, find their faith in numerical techniques to find out the approximate solution of a complicated model equation. In the recent era as a result of the technological revolution, numerical techniques have become the most desirable technique in searching the solution to such a problem that cannot be solved analytically.  There are several numerical methods to solve DEs.  However, the Euler method is one of the basic methods for solving initial value problems of ordinary differential equations (ODEs). Problems arising in DEs are of two types based on the condition given at the endpoints namely initial value problems (IVPs) and boundary value problems (BVPs), and it was introduced in 1768 by British mathematician Leonhard Euler (Hossain et al., 2017). After then several numerical methods developed for solving DEs namely Higher-order Taylor methods, Runge-Kutta methods, Multistep methods, Extrapolation methods, shooting method, Finite difference methods, Finite element methods, Finite volume methods, and so on, which may find any standard textbook (Burden & Faires, 1985; Balagurusamy, 2006; Mathews, 1992; Gerald, 2004; Eaqub, 2006; Lambert, 1973; Sastry, 2012; Iqbal et al., 2020). 

However, the researcher paid a little more attention to Runge-Kutta (RK) method to solve ODEs due to its simplicity to apply and able to provide accuracy. The method was introduced around 1900 by German mathematicians Carl Runge and Martin Wilhelm Kutta. There are several types of Runge-Kutta method (Gowri et al., 2017; Ababneh et al., 2009; Ademiluyi & Babatola, 2001; Akanbi, 2011; Agam & Yahaya, 2014; Chauhan & Srivastava, 2018; Dormand & Prince, 1980; Paul et al., 2014; Ralston, 1962; Senthi-lkumar, 2009; Paul et al., 2016; Senthilkumar & Paul, 2012; Yaakub & Evans, 1999; Paul & Senthil Kumar, 2015; Fehlberg, 1964; Fehlberg, 1966). Among them, RK fourth-order method is more popular. Nowadays it is widely used for solving DEs (Arora, 2019). 

In this article, we have used RK fourth order and RK sixth order method to solve initial and boundary value DEs. Also, we have to use the shooting method to convert the boundary value DE to the initial value DE to employ the RK method. The rest of the paper is organized as follows: Section 2 deals with a brief discussion of the RK method, in section 3 our selected problem is introduced, application of RK fourth order and sixth order method to the selected problem is presented in section 4 and 5, respectively, section 6 consists of result and discussion and conclusion is placed in section 7. Finally, the paper ends with a list of references.

2. RUNGE-KUTTA METHOD

Runge-Kutta method is a universal method that is widely used for solving ODE. The Runge-Kutta techniques are designed in such a way to achieve greater accuracy in approximation where functional values are required to know at some certain point of the subinterval. A brief discussion of the RK 4th order and 6th order method is presented in the adjacent subsections.

Runge-Kutta 4th order method 

Runge-Kutta 4th order (RK4) method is the most powerful method to solve ODE. The classical RK method of 4th order requires four evaluations per step and it gives more accurate results. A widely used technique of classical RK fourth order method was used in this study whose derivation is also presented here.

The Runge¬-Kutta method finds approximate value y for a given x. Only first-order ODEs can be solved by using the Runge-Kutta 4th order method. The formula used to compute the next value y_(n+1) from the previous valuey_n given below. The value of   n are   0,1,2,3,…...,s(x-x_0)/ h. Here, h is step height and x_(n+1)=x_n+h

y(x+h)-y(x)≅ak_1+bk_2+ck_3+dk_4

k_1=hf(x,y)

k_2=hf(x+mh,y+mk_1)

k_3=hf(x+mh,y+mk_2)

k_4=hf(x+ph,y+pk_1)

k_1 is the increment based on the slope at the beginning of the interval, using y

k_2 is the increment based on the slope at the midpoint of the interval, using y+(hk_1)/2

k_3 is again the increment based on the slope at the midpoint, using y+(hk_2)/2

k_4 is the increment based on the slope at the end of the interval, using y+hk_3

To find the coefficients a,b,c,d,m,n,p from Runge-Kutta formulas, re-product the Taylor series in term of h. The last formula is not a polynomial approximation.

F_1=f_x+ff_y

F_2=f_xx+2ff_xy+f^2 f_yy

F_3=f_xxx+ff_xxy+3f^2 f_xyy+f^3 f_yyy

Then differentiating the equation(y=) ́f(x,y), we get

y^2=f_x+f_y y^=f_x+f_y f=F_1

y^3=f_xx+2ff_xy+f^2 f_yy+f_y (f_x+ff_y )=F_2+f_y F_1

y^4=f_xxx+3ff_xxy+3f^2 f_xyy+f^3 f_yyy+f_y (f_xx+2f_xy+f^2 f_yy )+3(f_x+ff_y )(f_xx+ff_yy )+〖f_y〗^2 (f_x+ff_y )

=F_3+f_y F_2+3F_1 (f_xy+ff_yy+〖f_y〗^2 F_1)

Which follows the Taylors series is to be written in the form

y(x+h)-y(x)=hf+1/2 h^2 F_1+1/6 h^3 〖(F〗_2+f_y F_1)+1/24 h^4 [F_3+f_y F_2+3(f_xx+ff_yy ) F_1+〖f_x〗^2 F_1 ]

Turning now to the various k values similar com-putations produce

k_1=hf

k_2=h[f+mhF_1+1/2 m^2 h^2 F_2+1/6 m^3 h^3 F_3+⋯⋯⋯]

k_3=[{f+nhF_1+1/2 h^2 (n^2 F_2+2mnf_y F_1 )+1/6 h^3 (n^2 F_3+3m^2 nf_y F_2+6mn^2 (f_xy+ff_yy ) F_1 )}+6mnp〖f_y〗^2 F_1+⋯⋯⋯]

k_4=[{f+phF_1+1/2 h^2 (p^2 F_2+2pnf_y F_1 )+1/6 h^3 (p^3 F_3+3n^2 pf_y F_2+6np^2 (f_xy+ff_yy)F_1)}+6mnp〖f_y〗^2 F_1+⋯⋯⋯]

Combining these according to the final Runge-Kutta formula,

y(x+h)-y(x)=(a+b+c+d)hf+(bm+cn+dp) h^2 F_1+1/2 (bm^2+cn^2+dp^2 ) h^3 F^2+1/6 (bm^3+cn^3+dp^3 ) h^4 F_3+(cmn+dnp) h^3 f_y F_1+1/2 (cm^2 n+dn^2 p) h^4 f_y F_2+(cmn^2+dnp^2 ) h^4 (f_xy+ff_yy ) F_1+dmnph^4 〖f_y〗^2 F_1+⋯⋯⋯

Comparison with the Taylor series suggests the eight conditions

a+b+c+d=1

bm+cn+dp=1/2

bm^2+cn^2+dp^2=1/3

bm^3+cn^3+dp^3=1/4

cmn+dnp=1/6

cmn^2+dnp^2=1/8

cm^2 n+dn^2 p=1/12

dmnp=1/3

The eight equations in seven unknowns are actually excessive. The classical solution set is 

m=n=1/2

p=1

a=d=1/6

b=c=1/3

Leading to the Runge-Kutta formulas

k_1=hf(x,y)

k_2=hf(x+1/2 h,y+1/2 k_1 )

k_3=hf(x+1/2 h,y+1/2 k_2 )

k_4=hf(x+h,y+k_3 )

From equation (1) we get

y(x+h)=y(x)+1/6(k_1+〖2k〗_2+2k_3+k_4)

y_(r+1) (x+h)=y_r (x)+1/6(k_1+〖2k〗_2+2k_3+k_4)

Where r=1,2,3,…….,n, which is the general 4th order RK method.

Runge-Kutta 6th order method

Runge-Kutta Sixth order (RK6) method is better for solving an ODE as well as a partial differential equation (PDE) with an IVP. Though in general, higher-order Runge-Kutta gives a better solution than the lower order ones (Hossain et al., 2017). But it is not true for all purposes and differs from case to case. However, it is a well-known result that the Runge-Kutta method of order ρ needs at least ρ stages to the process. It is true for ρ=1,2,3,4 but suddenly for order ρ=5 there exist a process with stages 6 due to (Kutta, 1901) but corrected by (Nyström, 1925). Again, in the study of (Huïa, 1956) there are 8 stages processed for ρ=6. To mitigate, these varieties and complexity Butcher (Butcher, 1963) provide a common form to all Runge-Kutta method by a matrix and two vectors and it is known as Butcher tableau. This method approximates the solution to an ODE of the form y=f(t,y) with s stages is

y_(n+1)=y_n+h∑_(i=1)^s▒〖b_i k_ni 〗

With k_i=hf(x_n+c_i h,y_n+∑_(j=1)^(i-1)▒〖a_ij k_j 〗),

Where, n is the nodes running index (n=0,1,2,⋯⋯), i is the index used to label stages (1≤i≤s), the nodes c_i=∑_(j=1)^(i-1)▒a_ij  and for the consistency the weights b_i are defined as ∑_(i=1)^s▒b_i =1.

And the Butcher tableau is

0

c_1 a_21

c_2 a_31 a_32

c_3 a_41 a_42 a_43

⋮

c_s a_s1 a_s2 a_s3 ⋯ a_(s,(s-1))

b_1 b_2 b_3 ⋯ b_(s-1) b_s

Derivation of RK6 method from Butcher tableau can find in any standard numerical analysis book. Due to the paucity of space, we are not presenting the derivation here. For details of derivation please see (Butcher, 1963). And the Butcher tableau for RK6 method is -

0

1/3 1/3

2/3 0 2/3

1/3 1/12 1/3 -1/12

1/2 -1/16 9/8 -3/16 -3/8

1/2 0 9/8 -3/8 -3/4 1/2

1 9/44 -9/11 63/44 18/11 0 -16/11

11/120 0 27/40 27/40 -4/15 -4/15 11/120

Leading to the RK6 formulas

y_(n+1)=y_n+h/120 (11k_1+81k_3+81k_4-32k_5-32k_6+11k_7 ),

Where,

k_1=f(x_n,y_n,z_n )

k_2=f(x_n+h/3,y_n+k_1/3,z_n+l_1/3)

k_3=f(x_n+2h/3,y_n+(2k_2)/3)

k_4=f(x_n+h/3,y_n+k_1/12+k_2/3-k_3/12)

k_5=f(x_n+h/2,y_n-k_1/16+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/8)

k_6=f(x_n+h/2,y_n+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/4+k_5/2)

k_7=f(x_n+h,y_n+(9k_1)/44-〖9k〗_2/11+〖63k〗_3/44+〖18k〗_4/11-(16k_6)/11)

3. Problem statement

There are two parts of differential equations, IVP and BVP on the basics of the given conditions indicated at the endpoints. In our present work, second order initial and boundary value ODEs are solved by RK4 and RK6 methods. Among many different formulas, one most common form of fourth ordered which derived and Butchers sixth ordered Runge-Kutta methods are used. The BVP problems are firstly converted to a system of IVP through the medium of shooting method.

The numerical solution of IVPs

Consider the following second order IVP of ODE:

█((d^2 y)/(dx^2 )=f(x,y,y^ ),                        (1))

With the initial condition y(x_0 )=y_0,y^ (x_0 )=α.

Now this second order IVP can be written in a system of first order differential equations by considering

dy/dx=y^=z.

Then the new problem will be,

dz/dx=f(x,y,z)  dy/dx=g(x,y,z)

Where, g(x,y,z)=z

With the initial condition y(x_0 )=y_0,z(x_0 )=z_0=α.

Then the RK4 method will use the following formula: 

z(x_(n+1) )=z_n+h/6 (k_1+〖2k〗_2+2k_3+k4),

y(x_(n+1) )=y_n+h/6 (l_1+〖2l〗_2+2l_3+l_4 ),

Where,

k_1=f(x_n,y_n,z_n )

l_1=f(x_n,y_n,z_n )

k_2=f(x_n+h/2,y_n+k_1/2,z_n+l_1/2)

l_2=g(x_n+h/2,y_n+k_1/2,z_n+l_1/2)

k_3=f(x_n+h/2,y_n+k_2/2,z_n+l_2/2)

l_3=g(x_n+h/2,y_n+k_2/2,z_n+l_2/2)

k_4=f(x_n+h/2,y_n+k_3,z_n+l_3 )

l_4=g(x_n+h/2,y_n+k_3,z_n+l_3 )

For n=0,1,2,………………… .

And the RK6 method will use the following formula

z_(n+1)=z_n+h/120 (11k_1+81k_3+81k_4-32k_5-32k_6+11k_7 ),

y_(n+1)=y_n+h/120 (11l_1+81l_3+81l_4-32l_5-32l_6+11l_7 ),

Where,

k_1=f(x_n,y_n,z_n )

l_1=f(x_n,y_n,z_n )

k_2=f(x_n+h/3,y_n+k_1/3,z_n+l_1/3)

l_2=g(x_n+h/3,y_n+k_1/3,z_n+l_1/3)

k_3=f(x_n+2h/3,y_n+(2k_2)/3,z_n+〖2l〗_2/3)

l_3=g(x_n+2h/3,y_n+〖2k〗_2/3,z_n+〖2l〗_2/3)

k_4=f(x_n+h/3,y_n+k_1/12+k_2/3-k_3/12,z_n+l_1/12+l_2/3-l_3/12)

l_4=g(x_n+h/3,y_n+k_1/12+k_2/3-k_3/12,z_n+l_1/12+l_2/3-l_3/12)

k_5=f(x_n+h/2,y_n-k_1/16+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/8,z_n-l_1/16+(9l_2)/8-〖3l〗_3/16-〖3l〗_4/8)

l_5=g(x_n+h/2,y_n-k_1/16+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/8,z_n-l_1/16+(9l_2)/8-〖3l〗_3/16-〖3l〗_4/8)

k_6=f(x_n+h/2,y_n+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/4+k_5/2,z_n+(9l_2)/8-〖3l〗_3/16-〖3l〗_4/4+l_5/2)

l_6=g(x_n+h/2,y_n+(9k_2)/8-〖3k〗_3/16-〖3k〗_4/4+k_5/2,z_n+(9l_2)/8-〖3l〗_3/16-〖3l〗_4/4+l_5/2)

k_7=f(x_n+h,y_n+(9k_1)/44-〖9k〗_2/11+〖63k〗_3/44+〖18k〗_4/11-(16k_6)/11,z_n+(9l_1)/44-〖9l〗_2/11+〖63l〗_3/44+〖18l〗_4/11-(16l_6)/11)

l_7=g(x_n+h,y_n+(9k_1)/44-〖9k〗_2/11+〖63k〗_3/44+〖18k〗_4/11-(16k_6)/11,z_n+(9l_1)/44-〖9l〗_2/11+〖63l〗_3/44+〖18l〗_4/11-(16l_6)/11)

For n=1,2,3,………………………

The numerical solution of BVPs

Consider the following second order boundary value ODEs:

█((d^2 y)/(dx^2 )=f(x,y,y^ ),a≤x≤b             (2))

With the boundary condition y(a)=α,y(b)=β.

The same formulas of RK4 and RK6 methods that are used in solving IVP are used along with the shooting method to solve this BVP. Mainly, in the Shooting method, BVP has been transferred under consideration by an IVP with unknown parameter t as of the form,

█((d^2 y)/(dx^2 )=f(x,y,y^ );  a≤x≤b      (3))

y(a)=α,y^ (a)=t . 

The parameter t is to be chosen in a manner to ensure that   〖lim〗┬(k→∞)⁡〖y(b,t_k )=y(b)=β〗, where  y(x,t_k) denotes the solution to the initial-value problem with   t=t_k. Thus RK4 and RK6 methods are used for solving the IVPs of Eq. (3) for the repeated value of  t until we carry out the desired level of tolerance and accuracy and the solution y(x) of the BVP achieved.

4. Application of Runge-Kutta 4th order and 6th order method to the initial value problem

In this section, we consider the following examples of IVP to test the efficiency and convergence of num-erical solution along with actual solution where Runge-Kutta fourth order and sixth order methods were applied separately. The tables are containing the approximated result, actual result, and computed absolute error; the outcomes are displayed in the graph also. For better comparison, three different step sizes are used as 0.1, 0.05, and 0.025, respectively. The computations programming language, associated with the examples, is performed by MATLAB.

Example 1. Consider the IVP

y^=(12 log⁡x)/x^2 -y^/x,1≤x≤2,

y(1)=y^ (1)=1.

The actual solution of this equation is y=1+log⁡x+2〖(log⁡x)〗^3.

Table 1 to 3 and Fig 1 to 3 show the computed result and maximum absolute error relative to the actual solution. Where the line curve represents the approximated curve and the circle is for actual results.

Table 1:  Computed result and absolute maximum error of IVP with step size 0.1

5. Application of Runge-Kutta 4th order and 6th order method to the boundary value problem

The following BVP was taken for our study and test the efficiency and convergence of Runge-Kutta fourth order and sixth order methods separately and compared with actual solution. The tables are containing the approximated result, actual result, and computed absolute error; also the outcomes are displayed in the graph. For better comparison, three different step sizes are used as 0.1, 0.05, and 0.025, respectively. The computations programming lang-uage, for the affiliated examples, are performed by MATLAB.

Example 2: Consider the BVP〖y^〗^=2y^-y+xe^x-x,0≤x≤2, y(0)=0, y(2)=-4and compare the result to the actual solution y(x)=1/6 x^3 e^x-5/3 xe^x+2e^x-x-2.

Table 4-6, and Fig 4-6 show the computed result and maximum absolute error with actual solution. As earlier the line curved represents the approximated curve and circle for the actual results.

Table 4:  Computed result and the absolute maximum error of BVP with step size 0.1

6. Comparison between the obtained results

The obtained result of our model examples of IVP and BVP is expressed through Table 1-3 and Table 4-6 and graphical representations are in Fig 1-3 and Fig 4-6, respectively. From the tables and figures for both examples and both methods, we observed that if step size leads to deterioration the estimated solution converges faster to the exact solution such that in the limit when step size tends to zero the errors go to zero. Also, we observed from Table 1-3 that in solving IVP the method RK4 gives a more accurate and better result than the method RK6 and it has been clarified with the adjacent graphs Fig 1-3. 

We also see from Table 4-6 and Fig 4-6 that in solving BVP the integrator RK4 gives better approximation and converges faster to the actual solution than RK6.

In this paper, we have applied Runge-Kutta fourth order and Butchers sixth ordered methods to solve IVP and BVP of ODEs and the results are found in good agreement with exact solutions. We assume three different step sizes for each problem to arrive at a more accurate result with competition between the exact and approximated solutions in tables and figures. From the above investigation of these methods, observed that the convergence of approximated solutions to exact solutions increases with decreasing step size in both cases and the convergence rate of RK4 is superior to RK6 in comparison to the exact solutions. Hence it is cleared from this study that to find the more accurate result the lower order RK4 method is appropriate than the higher order RK6 method.

We would like to express our gratitude to our dear students RoksanaYeasmin, Sumaiya Pervin, Sayantani Biswas for their effortless support in completing this paper. We would also like to thanks our colleagues for their advice and support.

The authors declare that there is no conflict of interest concerning the publication of this paper.