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

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.

1. Introduction

A special third-order differential equation (ODE) of the form 𝑦^′′′(𝑥)=𝑓𝑥,𝑦(𝑥),𝑦𝑥₀𝑦=𝛼,_′𝑥₀=𝛽,𝑦^′′𝑥₀=𝛾,𝑥≥𝑥₀ (1) which is not explicitly dependent on the first derivative 𝑦^′(𝑥) and the second derivative 𝑦^′′(𝑥) 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 𝑦_𝑛+1, 𝑦^′_𝑛+1, and 𝑦^′′_𝑛+1 to approximate 𝑦(𝑥_𝑛+1), 𝑦^′(𝑥_𝑛+1), and 𝑦^′′(𝑥_𝑛+1), where 𝑦_𝑛+1 is the computed solution and 𝑦(𝑥_𝑛+1) is the exact solution.

2. Derivation of RKD Method

The general form of RKD method with 𝑠 stages for solving initial value problem (1) can be written as 𝑦_𝑛+1=𝑦_𝑛+ℎ𝑦^′_𝑛+ℎ²2𝑦_𝑛^′′+ℎ^3𝑠_𝑖=1𝑏_𝑖𝑘_𝑖,𝑦^′_𝑛+1=𝑦^′_𝑛+ℎ𝑦_𝑛^′′+ℎ^2𝑠_𝑖=1𝑏^′_𝑖𝑘_𝑖,𝑦^′′_𝑛+1=𝑦_𝑛^′′+ℎ^𝑠_𝑖=1𝑏_𝑖^′′𝑘_𝑖, (2) where 𝑘₁=𝑓𝑥_𝑛,𝑦_𝑛,𝑘_𝑖=𝑓𝑥_𝑛+𝑐_𝑖ℎ,𝑦_𝑛+ℎ𝑐_𝑖𝑦^′_𝑛+ℎ²2𝑐²_𝑖𝑦_𝑛^′′+ℎ^3𝑖−1_𝑗=1𝑎_𝑖𝑗𝑘_𝑗 (3)(4) for 𝑖=2,3,…,𝑠.

The parameters of the RKD method are 𝑐_𝑖,𝑎_𝑖𝑗,𝑏_𝑖,𝑏^′_𝑖,𝑏_𝑖^′′ for 𝑖=1,2,…,𝑠 and 𝑗=1,2,…,𝑠 are assumed to be real. If 𝑎_𝑖𝑗=0 for 𝑖⩽𝑗, 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)

To 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 𝑦_𝑛+1=𝑦_𝑛+ℎΦ𝑥_𝑛,𝑦_𝑛,𝑦^′_𝑛+1=𝑦^′_𝑛+ℎΦ^′𝑥_𝑛,𝑦_𝑛,𝑦^′′_𝑛+1=𝑦_𝑛^′′+ℎΦ^′′𝑥_𝑛,𝑦_𝑛, (6) where the increment functions are Φ𝑥_𝑛,𝑦_𝑛=𝑦^′_𝑛+ℎ2𝑦_𝑛^′′+ℎ^2𝑠_𝑖=1𝑏_𝑖𝑘_𝑖,Φ^′𝑥_𝑛,𝑦_𝑛=𝑦_𝑛^′′+ℎ^𝑠_𝑖=1𝑏^′_𝑖𝑘_𝑖,Φ^′′𝑥_𝑛,𝑦_𝑛=^𝑠_𝑖=1𝑏_𝑖^′′𝑘_𝑖, (7) and 𝑘_𝑖 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 𝑦(𝑥) of (1) into the RKD increment function. This gives 𝑡_𝑛+1=ℎ[Φ−Δ],𝑡^′_𝑛+1=ℎΦ^′−Δ^′,𝑡^′′_𝑛+1=ℎΦ^′′−Δ^′′. (8)

These expressions are best given in terms of elementary differentials and the Taylor series increment may be written as Δ=𝑦^′+12ℎ𝑦^′′+16ℎ²𝐹₁⁽³⁾+1ℎ24³𝐹₁⁽⁴⁾+𝑂ℎ⁴,Δ^′=𝑦^′′+12ℎ𝐹₁⁽³⁾+16ℎ²𝐹₁⁽⁴⁾+1ℎ24³𝐹₁⁽⁵⁾+𝑂ℎ⁴,Δ^′′=𝐹₁⁽³⁾+12ℎ𝐹₁⁽⁴⁾+16ℎ²𝐹₁⁽⁵⁾+𝑂ℎ³, (9) where for the scalar case the first few elementary differentials are 𝐹₁⁽³⁾𝐹=𝑓,₁₍₄₎=𝑓_𝑥+𝑓_𝑦𝑦^′,𝐹₁⁽⁵⁾=𝑓_𝑥𝑥+2𝑓_𝑥𝑦𝑦^′+𝑓_𝑦𝑦^′′+𝑓_𝑦𝑦𝑦^′2. (10)

Using the above terms, the increment functions Φ,Φ^′, and Φ^′′ for the RKD formula become ^𝑠_𝑖=1𝑏_𝑖𝑘_𝑖=^𝑠_𝑖=1𝑏_𝑖𝑓+^𝑠_𝑖=1𝑏_𝑖𝑐_𝑖𝑓_𝑥+𝑓_𝑦𝑦^′ℎ+12^𝑠_𝑖=1𝑏_𝑖𝑐²_𝑖𝑓_𝑥𝑥+2𝑓_𝑥𝑦𝑦^′+𝑓_𝑦𝑦^′′+𝑓_𝑦𝑦𝑦^′2ℎ²+𝑂ℎ³,_𝑠𝑖=1𝑏^′_𝑖𝑘_𝑖=^𝑠_𝑖=1𝑏^′_𝑖𝑓+^𝑠_𝑖=1𝑏^′_𝑖𝑐_𝑖𝑓_𝑥+𝑓_𝑦𝑦^′ℎ+12^𝑠_𝑖=1𝑏^′_𝑖𝑐²_𝑖𝑓_𝑥𝑥+2𝑓_𝑥𝑦𝑦^′+𝑓_𝑦𝑦^′′+𝑓_𝑦𝑦𝑦^′2ℎ²+𝑂ℎ³,_𝑠𝑖=1𝑏_𝑖^′′𝑘_𝑖=^𝑠_𝑖=1𝑏_𝑖^′′𝑓+^𝑠_𝑖=1𝑏_𝑖^′′𝑐_𝑖𝑓_𝑥+𝑓_𝑦𝑦^′ℎ+12^𝑠_𝑖=1𝑏_𝑖^′′𝑐²_𝑖𝑓_𝑥𝑥+2𝑓_𝑥𝑦𝑦^′+𝑓_𝑦𝑦^′′+𝑓_𝑦𝑦𝑦^′2ℎ²+𝑂ℎ³. (11)

The expressions for the local truncation errors in the 𝑦 solution, the first derivative, and the second derivative are 𝑡_𝑛+1=ℎ^3𝑠_𝑖=1𝑏_𝑖𝑘_𝑖−16𝐹₁⁽³⁾+1𝐹24¹⁽⁴⁾,𝑡+⋯_′𝑛+1=ℎ^2𝑠_𝑖=1𝑏^′_𝑖𝑘_𝑖−12𝐹₁⁽³⁾+16𝐹₁⁽⁴⁾,𝑡+⋯_′′𝑛+1=ℎ^𝑠_𝑖=1𝑏_𝑖^′′𝑘_𝑖−𝐹₁⁽³⁾+12𝐹₁⁽⁴⁾+16ℎ²𝐹₁⁽⁵⁾+⋯. (12)

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 𝑠-stage six-order RKD method can be written as follows.

Order conditions for 𝑦:?Order 3 𝑏_𝑖=16. (13)?Order 4 𝑏_𝑖𝑐_𝑖=1.24 (14)?Order 5 𝑏_𝑖𝑐²_𝑖=1.60 (15)?Order 6 𝑏_𝑖𝑐³_𝑖=1120,𝑏_𝑖𝑎_𝑖𝑗=1.720 (16)

Order conditions for 𝑦^′:?Order 2 𝑏^′_𝑖=12. (17)?Order 3 𝑏^′_𝑖𝑐_𝑖=16. (18)?Order 4 𝑏^′_𝑖𝑐²_𝑖=1.12 (19)?Order 5 𝑏^′_𝑖𝑐³_𝑖=120,𝑏^′_𝑖𝑎_𝑖𝑗=1.120 (20)?Order 6 𝑏^′_𝑖𝑐⁴_𝑖=130,𝑏^′_𝑖𝑎_𝑖𝑗𝑐_𝑗=1,720𝑏^′_𝑖𝑐_𝑖𝑎_𝑖𝑗=1.180 (21)

Order conditions for 𝑦^′′:?Order 1 𝑏_𝑖^′′=1. (22)?Order 2 𝑏_𝑖^′′𝑐_𝑖=12. (23)?Order 3 𝑏_𝑖^′′𝑐²_𝑖=13. (24)?Order 4 𝑏_𝑖^′′𝑐³_𝑖=14,𝑏_𝑖^′′𝑎_𝑖𝑗=1.24 (25)?Order 5 𝑏_𝑖^′′𝑐⁴_𝑖=15,𝑏_𝑖^′′𝑎_𝑖𝑗𝑐_𝑗=1120,𝑏_𝑖^′′𝑐_𝑖𝑎_𝑖𝑗=1.30 (26)?Order 6 𝑏_𝑖^′′𝑐²_𝑖𝑎_𝑖𝑗=136,𝑏_𝑖^′′𝑎_𝑖𝑗𝑐²_𝑗+𝑏_𝑖^′′𝑐_𝑖𝑎_𝑖𝑗𝑐_𝑗=7,720𝑏_𝑖^′′𝑐⁵_𝑖=16,𝑏_𝑖^′′𝑎_𝑖𝑗𝑐²_𝑗=1,360𝑏_𝑖^′′𝑐_𝑖𝑎_𝑖𝑗𝑐_𝑗=1,11442𝑏_𝑖^′′𝑎_𝑖𝑗𝑐²_𝑗+𝑏_𝑖^′′𝑐_𝑖𝑎_𝑖𝑗𝑐_𝑗=1.120 (27)

All indices are from 1 to 𝑠. To obtain the fifth-order RKD method, the following simplifying assumption is used in order to reduce the number of equations to be solved: 𝑏^′_𝑖=𝑏_𝑖^′′1−𝑐_𝑖,𝑖=1,…,𝑠. (28)

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 𝑦_𝑛,?𝑦^′_𝑛, and 𝑦_𝑛^′′ are given by ‖𝜏⁽⁶⁾‖₂=4.392052306×10⁻⁴,‖𝜏^′(6)‖₂=4.108422255×10⁻³,??‖𝜏^′′(6)‖₂=4.925724188×10⁻³, respectively, where 𝜏⁽⁶⁾,?𝜏^′(6), and 𝜏^′′(6) are error terms of the sixth-order conditions for 𝑦_𝑛,?𝑦^′_𝑛, and 𝑦_𝑛^′′, respectively.

The RKD5 Method. Consider the following:0035+√61027+√5001960310005−√61033−√2500516500051−√1250116012501118−√186148+√18648197−√366718+√366181949−√64369+√636 (29)

Next, 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: 100𝑦010001^𝑛+1ℎ𝑦^′_𝑛ℎ²𝑦^′′_𝑛+1=1112𝑦011001^𝑛+1ℎ𝑦^′_𝑛ℎ²𝑦^′′_𝑛+1,𝑝(𝜉)1=det𝐼𝜉−𝐴=det𝜉−1−1−20𝜉−1−1.00𝜉−1 (30) Thus, the characteristic polynomial is 𝑝(𝜉)=(𝜉−1)³. (31) 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 𝑑³𝑦(𝑥)𝑑𝑥³=38(1+𝑥)⁵,𝑦(0)=1,𝑦^′(0)=12,𝑦^′′(0)1=−4. (32) Exact solution: 𝑦(𝑥)=√1+𝑥.

Problem 2. Consider 𝑑³𝑦(𝑥)𝑑𝑥³=𝑦²(𝑥)+cos²(𝑥)(𝑥)(0)𝑦−cos−1,𝑦=0,_′(0)=1,𝑦^′′(0)=0. (33) Exact solution: 𝑦(𝑥)=sin(𝑥).

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) for the film profile 𝑦(𝑥) in a coordinate frame moving with the fluid. The form of 𝑓(𝑦) varies according to the physical context. Different forms of the function 𝑓 are studied in Tuck and Schwartz [16]. For drainage down a dry surface, the form of 𝑓(𝑦) was given as 𝑑³𝑦𝑑𝑥³1=−1+𝑦². (35)

When the surface is prewetted by a thin film with thickness 𝛿>0 (where 𝛿>0 is very small), the function 𝑓 is given by 𝑓𝑦=−1+1+𝛿+𝛿²𝑦²−𝛿+𝛿²𝑦³. (36)

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, 𝑦=𝑦(𝑥). According to Tuck and Schwartz [16], one such equation is 𝑦^′′′=𝑦^−𝑘,𝑥≥𝑥₀ (37) with initial conditions 𝑦𝑥₀=𝛼,𝑦^′𝑥₀=𝛽,𝑦^′′𝑥₀=𝛾, (38) 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) where 𝑦₁(0)=1,𝑦₂(0)=1,𝑦₃(0)=1, (40) we have taken 𝑥₀=0 and 𝛼=𝛽=𝛾=1.

Unfortunately, for general 𝑘, (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 𝑘=2 and 𝑘=3. The results are displayed in Tables 1 and 2 for the case 𝑘=2 and Tables 3 and 4 for the case 𝑘=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 at 𝑥∈[0,0.2,0.4,0.6,0.8,1.0] taking ℎ=0.01 and 𝑘=2 with the initial conditions 𝑦(0)=𝑦^′(0)=𝑦^′′(0)=1.

Table 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 at 𝑥∈[0,0.2,0.4,0.6,0.8,1.0] taking ℎ=0.1 and 𝑘=2 with the initial conditions 𝑦(0)=𝑦^′(0)=𝑦^′′(0)=1.

Table 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 at 𝑥∈[0,0.2,0.4,0.6,0.8,1.0] taking ℎ=0.01 and 𝑘=3 with the initial conditions 𝑦(0)=𝑦^′(0)=𝑦^′′(0)=1.

Table 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 at 𝑥∈[0,0.2,0.4,0.6,0.8,1.0] taking ℎ=0.1 and 𝑘=3 with the initial conditions 𝑦(0)=𝑦^′(0)=𝑦^′′(0)=1.

For 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 Figures 1 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 Problem 1 with 𝑥end=10 and ℎ=1/2𝑖, 𝑖=1,…,5.

Figure 2 

The efficiency curves of the RKD5 method and its comparisons for Problem 2 with 𝑥end=𝜋 and ℎ=1/2𝑖, 𝑖=1,…,5.

From Tables 1 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 𝑘=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, 𝑦_𝑖 for 𝑘=2 and 𝑘=3, respectively, with ℎ=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 solution 𝑦𝑖 for problem (37) for 𝑘=2, ℎ=0.01.

Figure 4 

Plot of the solution 𝑦𝑖 for problem (37) for 𝑘=3, ℎ=0.01.

Figure 5 

Plot of graph for function evaluations against step size, ℎ for 𝑘=3, ℎ=1/10𝑖, 𝑖=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.