# A Three-Stage Fifth-Order Runge-Kutta Method for Directly Solving Special Third-Order Differential Equation with Application to Thin Film Flow Problem

Authors: M. Mechee; N. Senu; F. Ismail; B. Nikouravan; Z. Siri
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

In this paper, a three-stage fifth-order Runge-Kutta method for the integration of
a special third-order ordinary differential equation (ODE) is constructed. The zero stability
of the method is proven. The numerical study of a third-order ODE arising
in thin film flow of viscous fluid in physics is discussed. The mathematical model of
thin film flow has been solved using a new method and numerical comparisons are
made when the same problem is reduced to a first-order system of equations which
are solved using the existing Runge-Kutta methods. Numerical results have clearly
shown the advantage and the efficiency of the new method.


---

## Body

## 1. Introduction

A special third-order differential equation (ODE) of the form(1)y′′′(x)=f(x,y(x)),y(x0)=α,y′(x0)=β,y′′(x0)=γ,x≥x0
which is not explicitly dependent on the first derivative y′(x) and the second derivative y′′(x) of the solution is frequently found in many physical problems such as electromagnetic waves, thin film flow, and gravity driven flows. The solution to (1) can be obtained by reducing it to an equivalent first-order system which is three times the dimension and can be solved using a standard Runge-Kutta method or a multistep method.Most researchers, scientists, and engineers solve problem (1) by converting the problem to a system of first-order equations. However, there are also studies on numerical methods which solve (1) directly. Such work can be seen in Awoyemi [1], Waeleh et al. [2], Zainuddin [3], and Jator [4]. Awoyemi and Idowu [5] and Jator [6] proposed a class of hybrid collocation methods for the direct solution of higher-order ordinary differential equations (ODEs). Samat and Ismail [7] developed an embedded hybrid method for solving special second-order ODEs. Waeleh et al. [2] developed a block multistep method which can directly solve general third-order equations; on the other hand, Ibrahim et al. [8] developed a multistep method that can directly solve stiff third-order differential equations. All of the methods discussed above are multistep methods; hence, they need starting values when used to solve ODEs such as (1). Senu et al. [9] derived the Runge-Kutta-Nyström method for solving second-order ODEs directly. Mechee et al. [10] constructed new Runge-Kutta methods for solving (1).In this paper, we are concerned with a one-step method, particularly the three-stage fifth-order Runge-Kutta method, for directly solving special third-order ODEs. Accordingly, we have developed a direct Runge-Kutta method (RKD) which can be directly used to solve (1). The advantage of the new method over multistep methods is that it initialises itself. The method produces yn+1, yn+1′, and yn+1′′ to approximate 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.


## 2. Derivation of RKD Method

The general form of RKD method withs stages for solving initial value problem (1) can be written as
(2)yn+1=yn+hyn′+h22yn′′+h3∑i=1s‍biki,yn+1′=yn′+hyn′′+h2∑i=1s‍bi′ki,yn+1′′=yn′′+h∑i=1s‍bi′′ki,
where
(3)k1=f(xn,yn),(4)ki=f(xn+cih,yn+hciyn′+h22ci2yn′′+h3∑j=1i-1‍aijkj)
for i=2,3,…,s.The parameters of the RKD method areci,aij,bi,bi′,bi′′ for i=1,2,…,s and j=1,2,…,s are assumed to be real. If aij=0 for i⩽j, it is an explicit method, and otherwise it is an implicit method. The RKD method can be expressed in the Butcher notation using the table of coefficients as follows:
(5)cAbTb′Tb′′TTo determine the coefficients of the method given by (2)–(4), the RKD method expression is expanded using Taylor’s series expansion. After some algebraic manipulations, this expansion is equated to the true solution that is given by Taylor’s series expansion. General order conditions for the RKD method can be found from the direct expansion of the local truncation error. This idea is based on the derivation of order conditions for the Runge-Kutta method introduced in Dormand [11]. The RKD formulae in (2) may be expressed as
(6)yn+1=yn+hΦ(xn,yn),yn+1′=yn′+hΦ′(xn,yn),yn+1′′=yn′′+hΦ′′(xn,yn),
where the increment functions are
(7)Φ(xn,yn)=yn′+h2yn′′+h2∑i=1s‍biki,Φ′(xn,yn)=yn′′+h∑i=1s‍bi′ki,Φ′′(xn,yn)=∑i=1s‍bi′′ki,
and ki is defined in formula (4). If Δ is the Taylor series increment function, then the local truncation errors of the solution, the first derivative, and the second derivative may be obtained by substituting the true solution y(x) of (1) into the RKD increment function. This gives
(8)tn+1=h[Φ-Δ],tn+1′=h[Φ′-Δ′],tn+1′′=h[Φ′′-Δ′′].These expressions are best given in terms of elementary differentials and the Taylor series increment may be written as(9)Δ=y′+12hy′′+16h2F1(3)+124h3F1(4)+O(h4),Δ′=y′′+12hF1(3)+16h2F1(4)+124h3F1(5)+O(h4),Δ′′=F1(3)+12hF1(4)+16h2F1(5)+O(h3),
where for the scalar case the first few elementary differentials are
(10)F1(3)=f,F1(4)=fx+fyy′,F1(5)=fxx+2fxyy′+fyy′′+fyy(y′)2.Using the above terms, the increment functionsΦ,Φ′, and Φ′′ for the RKD formula become
(11)∑i=1s‍biki=∑i=1s‍bif+∑i=1s‍bici(fx+fyy′)h+12∑i=1s‍bici2(fxx+2fxyy′+fyy′′+fyy(y′)2)h2+O(h3),∑i=1s‍bi′ki=∑i=1s‍bi′f+∑i=1s‍bi′ci(fx+fyy′)h+12∑i=1s‍bi′ci2(fxx+2fxyy′+fyy′′+fyy(y′)2)h2+O(h3),∑i=1s‍bi′′ki=∑i=1s‍bi′′f+∑i=1s‍bi′′ci(fx+fyy′)h+12∑i=1s‍bi′′ci2(fxx+2fxyy′+fyy′′+fyy(y′)2)h2+O(h3).The expressions for the local truncation errors in they solution, the first derivative, and the second derivative are
(12)tn+1=h3[∑i=1s‍biki-(16F1(3)+124F1(4)+⋯)],tn+1′=h2[∑i=1s‍bi′ki-(12F1(3)+16F1(4)+⋯)],tn+1′′=h[∑i=1s‍bi′′ki-(F1(3)+12F1(4)+16h2F1(5)+⋯)].Substituting (11) into (12) and expanding as a Taylor expansion using MAPLE 4 software as introduced by Gander and Gruntz [12], the error equations or the order conditions for s-stage six-order RKD method can be written as follows.Order conditions fory:Order 3(13)∑‍bi=16.Order 4(14)∑‍bici=124.Order 5(15)∑‍bici2=160.Order 6(16)∑‍bici3=1120,∑‍biaij=1720.Order conditions fory′:Order 2(17)∑‍bi′=12.Order 3(18)∑‍bi′ci=16.Order 4(19)∑‍bi′ci2=112.Order 5(20)∑‍bi′ci3=120,∑‍bi′aij=1120.Order 6(21)∑‍bi′ci4=130,∑‍bi′aijcj=1720,∑‍bi′ciaij=1180.Order conditions fory′′:Order 1(22)∑‍bi′′=1.Order 2(23)∑‍bi′′ci=12.Order 3(24)∑‍bi′′ci2=13.Order 4(25)∑‍bi′′ci3=14,∑‍bi′′aij=124.Order 5(26)∑‍bi′′ci4=15,∑‍bi′′aijcj=1120,∑‍bi′′ciaij=130.Order 6(27)∑‍bi′′ci2aij=136,∑‍bi′′aijcj2+∑‍bi′′ciaijcj=7720,∑‍bi′′ci5=16,∑‍bi′′aijcj2=1360,∑‍bi′′ciaijcj=1144,12∑‍bi′′aijcj2+∑‍bi′′ciaijcj=1120.All indices are from 1 tos. To obtain the fifth-order RKD method, the following simplifying assumption is used in order to reduce the number of equations to be solved:
(28)bi′=bi′′(1-ci),i=1,…,s.To derive the three-stage fifth-order RKD method, we use the algebraic conditions up to order five (13)–(15), (17)–(20), and (22)–(26). The resulting system of equations consists of 16 nonlinear equations with 13 unknown variables for which we need to solve. Consequently, there is no solution since the number of equations exceeds the number of unknowns to be solved for. To overcome this, the simplifying assumption (28) is imposed. This will reduce the number of equations to 11 with 11 unknowns making the system solvable. This system has no free parameter but is consistent and yields a unique solution. The coefficients of the method are given below (see The RKD5 method in (29)). The error norms for yn, yn′, and yn′′ are given by ∥τ(6)∥2=4.392052306×10-4,∥τ′(6)∥2=4.108422255×10-3,  ∥τ′′(6)∥2=4.925724188×10-3, respectively, where τ(6), τ′(6), and τ′′(6) are error terms of the sixth-order conditions for yn, yn′, and yn′′, respectively.The RKD5 Method. Consider the following:(29)0035+61027500+1961000035-610332500-5165000511250-11612500118118-648118+64819736-618736+6181949-63649+636Next, we will discuss the zero stability of the method which is one of the criteria for the method to be convergent. Zero stability is an important tool for proving the stability and convergence of linear multistep methods. The interested reader is referred to the textbooks by Lambert [13] and Butcher [14] in which zero stability is discussed. Zero stability has been discussed in Hairer et al. [15] to determine an upper bound on the order of convergence of linear multistep methods. In studying the zero stability of the RKD method, we can write method (2) as follows:
(30)(100010001)(yn+1hyn′h2yn+1′′)=(1112011001)(yn+1hyn′h2yn+1′′),p(ξ)=det[Iξ-A]=det(ξ-1-1-120ξ-1-100ξ-1).
Thus, the characteristic polynomial is
(31)p(ξ)=(ξ-1)3.
Hence, the method is zero stable since the roots, ξ=1, are less than or equal to one.


## 3. Problems Tested and Numerical Results

In this section, we will apply the new method to some third-order differential equation problems. The following explicit RK methods are selected for the numerical comparisons.(i)
RKD5: the three-stage fifth-order RKD method derived in this paper.(ii)
RK4: the four-stage fourth-order RK method as in Butcher [14].(iii)
DOPRI: the seven-stage fifth-order Runge-Kutta derived by Dormand [11].Problem 1.
Consider(32)d3y(x)dx3=38(1+x)5,y(0)=1,y′(0)=12,y′′(0)=-14.
Exact solution: y(x)=1+x.Problem 2.
Consider(33)d3y(x)dx3=y2(x)+cos2(x)-cos(x)-1,y(0)=0,y′(0)=1,y′′(0)=0.
Exact solution: y(x)=sin(x).


## 4. An Application to a Problem in Thin Film Flow

In this section, we will apply the proposed method to a well-known problem in physics regarding the thin film flow of a liquid. Many problems have discussed this problem; see, for example, [16–23]. In a survey paper, Tuck and Schwartz [16] discussed the flow of a thin film of viscous fluid over a solid surface. Tension and gravity as well as viscosity are taken into account. The problem was formulated using third-order ordinary differential equations (ODEs) as follows:
(34)d3ydx3=f(y)
for the film profile y(x) in a coordinate frame moving with the fluid. The form of f(y) varies according to the physical context. Different forms of the function f are studied in Tuck and Schwartz [16]. For drainage down a dry surface, the form of f(y) was given as
(35)d3ydx3=-1+1y2.When the surface is prewetted by a thin film with thicknessδ>0 (where δ>0 is very small), the function f is given by
(36)f(y)=-1+1+δ+δ2y2-δ+δ2y3.Problems concerning the flow of thin films of viscous fluid with a free surface in which surface tension effects play a role typically lead to third-order ordinary differential equations governing the shape of the free surface of the fluid,y=y(x). According to Tuck and Schwartz [16], one such equation is
(37)y′′′=y-k,x≥x0
with initial conditions
(38)y(x0)=α,y′(x0)=β,y′′(x0)=γ,
where α, β, and γ are constants, is of particular importance because it describes the dynamic balance between surface and viscous forces in a thin fluid layer in the absence (or neglect) of gravity.For comparison purposes, we will use Runge-Kutta methods which are fourth-order (RK4) and fifth-order (DOPRI) methods, respectively. To use Runge-Kutta methods we write (1) as a system of three first-order equations. Following Biazar et al. [17], we can write (37) as the following system:
(39)dy1dx=y2(x),dy2dx=y3(x),dy3dx=y1-k(x),
where
(40)y1(0)=1,y2(0)=1,y3(0)=1,
we have taken x0=0 and α=β=γ=1.Unfortunately, for generalk, (37) cannot be solved analytically. We can, however, use these reductions to determine an efficient way to solve (1) numerically. We focus on the cases k=2 and k=3. The results are displayed in Tables 1 and 2 for the case k=2 and Tables 3 and 4 for the case k=3.Table 1
Table comparing values of the numerical solution, a fifth-order Runge-Kutta method (DOPRI), fourth-order Runge-Kutta method (RK4), and our new RKD5 method atx∈[0,0.2,0.4,0.6,0.8,1.0] taking h=0.01 and k=2 with the initial conditions y(0)=y′(0)=y′′(0)=1.
x
Exact solution
DOPRI
RK4
RKD5
0.0
1.000000000
1.0000000000
1.0000000000
1.0000000000
0.2
1.221211030
1.2212100045
1.2212100046
1.2212100045
0.4
1.488834893
1.4888347799
1.4888347800
1.4888347799
0.6
1.807361404
1.8073613977
1.8073613978
1.8073613977
0.8
2.179819234
2.1798192339
2.1798192341
2.1798192339
1.0
2.608275822
2.6082748676
2.6082748678
2.6082748676Table 2
Table comparing values of the numerical solution, a fifth-order Runge-Kutta method (DOPRI), fourth-order Runge-Kutta method (RK4), and our new RKD5 method atx∈[0,0.2,0.4,0.6,0.8,1.0] taking h=0.1 and k=2 with the initial conditions y(0)=y′(0)=y′′(0)=1.
x
Exact solution
DOPRI
RK4
RKD5
0.0
1.000000000
1.0000000000
1.0000000000
1.0000000000
0.2
1.221211030
1.2212100041
1.2212105060
1.2212099592
0.4
1.488834893
1.4888347796
1.4888356990
1.4888344801
0.6
1.807361404
1.8073613979
1.8073626884
1.8073605532
0.8
2.179819234
2.1798192349
2.1798208831
2.1798175255
1.0
2.608275822
2.6082748696
2.6082768844
2.6082719667Table 3
Table comparing values of the numerical solution using a fifth-order Runge-Kutta method (DOPRI), fourth-order Runge-Kutta method (RK4), and our new RKD5 method atx∈[0,0.2,0.4,0.6,0.8,1.0] taking h=0.01 and k=3 with the initial conditions y(0)=y′(0)=y′′(0)=1.
x
RK4
DOPRI
RKD5
0.0
1.0000000000
1.0000000000
1.0000000000
0.2
1.2211551425
1.2211551424
1.2211551423
0.4
1.4881052844
1.4881052842
1.4881052838
0.6
1.8042625484
1.8042625481
1.8042625471
0.8
2.1715227984
2.1715227981
2.1715227960
1.0
2.5909582594
2.5909582591
2.5909582556Table 4
Table comparing values of the numerical solution using a fifth-order Runge-Kutta method (DOPRI), fourth-order Runge-Kutta method (RK4), and our new RKD5 method atx∈[0,0.2,0.4,0.6,0.8,1.0] taking h=0.1 and k=3 with the initial conditions y(0)=y′(0)=y′′(0)=1.
x
RK4
DOPRI
RKD5
0.0
1.0000000000
1.0000000000
1.0000000000
0.2
1.2211559590
1.2211551421
1.2211550887
0.4
1.4881067401
1.4881052848
1.4881049238
0.6
1.8042645823
1.8042625503
1.8042615558
0.8
2.1715254210
2.1715228023
2.1715208324
1.0
2.5909615178
2.5909582657
2.5909549758For comparison purposes, we present efficiency curve where the common logarithm of the maximum global error along the integration versus the function evaluations as shown in Figures1 and 2. From Figures 1 and 2, we observed that the method RKD5 performs better when integrating third-order ODE compared to RK4 and DOPRI methods.Figure 1
The efficiency curves of the RKD5 method and its comparisons for Problem1 with xend=10 and h=1/2i, i=1,…,5.Figure 2
The efficiency curves of the RKD5 method and its comparisons for Problem2 with xend=π and h=1/2i, i=1,…,5.From Tables1 and 2 we observe that the numerical results using RKD5 are correct to five decimal places. Applying the fourth- and fifth-order RK methods (RK4 and DOPRI) to (39) for k=2 also yields five-decimal-place accuracy. While from Tables 3 and 4 we observed the numerical results for the new method, RKD5 agrees to eight decimal places when compared with the fourth- and fifth-order RK methods (RK4 and DOPRI). This is consistent with the results displayed in Tables 1 and 2. In Figures 3 and 4, we plot the numerical solution, yi for k=2 and k=3, respectively, with h=0.01. Figure 5 shows that the new RKD5 method requires less function evaluations than the RK4 and DOPRI methods. This is because when problem (37) is solved using RK4 and DOPRI method, it need to be reduced to a system of first-order equations which is three times the dimension.Figure 3
Plot of the solutionyi for problem (37) for k=2, h=0.01.Figure 4
Plot of the solutionyi for problem (37) for k=3, h=0.01.Figure 5
Plot of graph for function evaluations against step size,h for k=3, h=1/10i, i=1,…,4.


## 5. Discussion and Conclusion

In this paper, we have derived the order conditions for a Runge-Kutta method which can be used to directly solve special third-order ordinary differential equations. A three-stage fifth-order RKD method has been derived and applied to the thin film flow problem. Numerical results show that the new method agrees very well with well-known existing methods in the literature and required less function evaluations. As such, this method is more cost effective in terms of computation time than other existing methods.



---

*Source: 795397-2013-06-06.xml*