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

1 Institute of Mathematical Sciences, University of Malaya, 50603 Kuala Lumpur, Malaysia 2 Department of Mathematics, Faculty of Mathematics and Computer Sciences, University of Kufa, Najaf, Iraq 3 Department of Mathematics and Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia 4Department of Physics, Faculty of Science, University of Malaya, 50603 Kuala Lumpur, Malaysia


Introduction
A special third-order differential equation (ODE) of the form   () =  (,  ()) ,  ( 0 ) = ,   ( 0 ) = ,   ( 0 ) = ,  ≥  0 (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 thirdorder 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.

Derivation of RKD Method
The general form of RKD method with  stages for solving initial value problem (1) can be written as where for  = 2, 3, . . ., .
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: 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 where the increment functions are 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 These expressions are best given in terms of elementary differentials and the Taylor series increment may be written as where for the scalar case the first few elementary differentials are Using the above terms, the increment functions Φ, Φ  , and Φ  for the RKD formula become The expressions for the local truncation errors in the  solution, the first derivative, and the second derivative are 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 4 Order 5 Order 6 Order conditions for   : Order 2 Order 3 Order 4 Order 5 Order 6 Order conditions for   : Order 1 Order 2 Order 3 Order 4 Order 5 Order 6 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: 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 ‖ (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 sixthorder conditions for   ,    , and    , respectively.
The RKD5 Method.Consider the following: 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: ) .
(30) Thus, the characteristic polynomial is Hence, the method is zero stable since the roots,  = 1, are less than or equal to one.

Problems Tested and Numerical Results
In this section, we will apply the new method to some thirdorder 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.

An Application to a Problem in Thin Film Flow
In this section, we will apply the proposed method to a wellknown problem in physics regarding the thin film flow of a liquid.Many problems have discussed this problem; see, for example, [16][17][18][19][20][21][22][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: 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 When the surface is prewetted by a thin film with thickness  > 0 (where  > 0 is very small), the function  is given by 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 with initial conditions 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: where we have taken  0 = 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.
For comparison purposes, we present efficiency curve where the common logarithm of the maximum global error 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.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.

𝑥
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.

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 threestage 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.

Table 4 :
Table comparing values of the numerical solution using a fifth-order Runge-Kutta method (DOPRI), fourth-order Runge-