Approximate Analytical Solution of a Third-Order IVP arising in Thin Film Flows driven by Surface Tension

In this paper, we present a way of applying the so-called He’s variational iteration method (VIM) to numerically solve the non linear autonomous third-order ordinary differential equation (ODE) y = y obtained by considering a traveling wave solution admitted by a lubrication equation modeling a twodimensional spreading of a thin viscous film on a inclined slope. Approximate analytical solution is derived and compared to the results obtained from the Adomian decomposition method (ADM) proposed in [20], to the exact analytical solution [7,8], to a fifth order Runge-Kutta method (DOPRI), a fourth order Runge-Kutta method (RK4), a three-stage fifth order Runge-Kutta method (RKD5) developed in [18]. A very good agreement and accuracy is observed. Comparisons are obtained using symbolic capabilities of Maple 18.0 package.


Introduction 117
2 Problem statement and mathematical formulation 119 3 Basic concepts of He's variational iteration method 120

Introduction
In fluid dynamics of viscous fluids, the exact analytical solutions of the flow problems are usually difficult to obtain since governing equations of motion, in general, are highly nonlinear partial differential equations.In some situations, by means of similarity transformations, the system of partial differential equations are reduced to that of ordinary differential equations, which, on few occasions, for a long time, before the advent computers, the researchers mainly directed their efforts at obtaining some forms of approximate solutions.One of the key issues of approximate solutions has always been the accuracy of the solutions.The accuracy, generally speaking, is measured in terms of the norm of the relative percent error Fabrizio Morlando δ where the relative error being the difference of the approximate solution from the exact solution divided by the exact solution.In the absence of an exact solution, (analytical or numerical) a heuristic approach consisting of the convergence of successive approximate solution.With the advent of computers, the approximate solutions in fluid dynamics have lost some of their importance as more and better numerical algorithms have been developed to solve the increasingly realistic, but more but more complicated problems numerically.Nevertheless, approximate analytical solutions still have their relevance because, firstly, they give the solution for each point within the domain of interest, unlike the numerical solutions, which are available, for a particular run, only for a set of discrete points in the domain, secondly, compared to a numerical solution, a nicely produced approximate solution, requiring a minimal effort and having a reasonable amount of accuracy, is always handy for an engineer, scientist or an applied mathematician, who can obtain a solution quickly, thereby gaining a valuable insight into the essential of the problem and, thirdly, even with most of the scientific packages, some initial guess is required for the solution, as the algorithm, in general, are nor globally convergent.In such situations, approximate solutions can provide an excellent starting guess, that can be readily refined to the exact numerical solution in a few iterations.Problems concerning the flow of thin films of viscous fluid with a free surface in which surface tension effects play a key role typically lead to third-order ordinary differential equations governing the shape of the free surface of the fluid y = y(x).Among these, as pointed out in [33], it is of particular interest the third order ODE obtained by investigating traveling-wave solutions or steady-state solutions of the lubrication equation subject to the initial conditions where ′ = d/dx.This equation is of particular importance, because it describes the thickness y = y(x) of a fluid layer draining down a vertical wall dynamic balanced between surface tension effects, represented analytically by the thirdorder derivative term, and viscous shearing forces, represented by the term y −2 , in the absence (or neglect) of gravity.Tanner, in [30], also derived this equation in order to investigate the motion of the contact line for a thin oil droplet spreading on a horizontal surface.We may consider equation (1.1) as an intermediate limit problem being the small-y limit of a fluid dry wall draining model where the extra constant term −1 represents the presence of gravity and the large-y limit of its generalization to wet wall draining model, Approximate Analytical Solution of a Third-Order IVP 119 where δ > 0 stands for a very small input parameter measuring wall wetness.Remark that when δ = 0 (1.4) reduces to (1.3).The interested reader is referred to [8,33] and the references therein for a review of thin film flow theory in which some third-order ODEs occurring.
In this work, we aim to apply the variational iteration method (abbr.VIM) in finding the approximate analytical solutions of a highly nonlinear differential equation arising in the thin film flow problems.The VIM is a powerful analytical technique introduced by J.H.He in [10,11,12,13] as a modification of a general Lagrange multiplier method.This method gives rapidly convergent successive approximations of the exact solutions if such solution exists.We propose here this method for solving numerically equation (1.1) with condition (1.2) and, then, we compare the results with some exact and numerical results well-known in literature.Notice that in [20] the problem is also solved by the Adomian decomposition method (abbr.ADM).
For comparison purpose, we show that difference between the VIM and ADM solutions is negligible.This comparison is benched-marked against numerical solutions.An important advantage of these methods over the numerical methods is that they provide series solutions in the form of functions of a single variable.These forms of solutions can be used to evaluate analytical expressions of various flow parameters of physical relevance.The convergence of the VIM and the ADM is systematically discussed in [32].
This paper is divided up as follows.In Section 2, we formulate the concerned mathematical model.In Section 3, we give a brief introduction of the He's variational iteration method.In Section 4, we apply the variational iteration method to give an approximate solution of the considered problem.In Section 5, the results obtained are compared with the numerical (Runge-Kutta type methods) results.Furthermore, we show that our approximated solution is in good agreement with the previous results presented as far in literature and the absolute percent error of approximation slightly improves.Concluding remarks are made in Section 6.

Problem statement and mathematical formulation
In the contest of standard lubrication theory, typically surface tension driven thin film flows are described by a fourth order nonlinear parabolic partial differential equation of the form modeling the two-dimensional spreading of a thin viscous film on a slope inclined at an angle α to the horizontal, where C is the inverse capillary number (i.e. the ratio of surface tension to viscous forces), B the Bond number (i.e. the ratio of gravity to viscous forces), δ the aspect ratio and h = h(t, x) the time evolution fluid film height, see ( [8], Appendix) for a derivation.For the flow of a thin film down a vertical wall, α = π/2 and (2.1) reduces to For steady situations, it is well-known that (2.2) may be integrated once and a third-order ordinary differential equation is obtained.Henceforth, we deal with traveling wave solutions, i.e. solutions of the form admitted by (2.2), where V is the wave velocity.The investigation of the traveling wave solution (2.3) can be interpreted as investigating steady-state solutions with respect to a moving frame where the frame moves with velocity V , see [30].Substituting (2.3) into (2.2) we obtain the fourth order ordinary differential equation where, in this case, ′ = d/dx.Integrating (2.4) once and setting the integration constant to zero, B = 0 and C = 3V we obtain (1.1) where the overhead bars are suppressed for convenience.Another derivation of (1.1) comes from considering traveling wave solutions of the lubrication equation: 3) into (2.5)we obtain the third-order ODE where c is a constant of integration.Finally, (2.6) is rescaled to (1.2), see [3].Duffy and Wilson in [7] use the solution of (1.1) obtained by Ford [8] (see also [22]) to plot solutions of (1.1) subject to the initial conditions (1.2).They, then, show that the critical initial condition y ′′ (0) = 1.28359871 leads to the boundary condition that is also obtained by Tanner [30] and Tuck and Schwartz [33].In this work, we choose y ′′ (0) = 1 to simplify our calculations.

Basic concepts of He's variational iteration method
For the purpose of illustration of the methodology to the proposed method using variational iteration method (VIM), we begin by considering a non-linear differential equation given in the following formal form where L = d m dt m , m ∈ N, is a linear operator, N a non linear operator and g is known analytical non-homogeneous source term, subjected to the initial conditions Approximate Analytical Solution of a Third-Order IVP 121 where c k is a real number.According to the variational iteration method, a correction functional can be constructed as where λ = λ(t) is a general Lagrangian multiplier, the subscript n denotes the nth order approximation and u 0 is an initial approximation which can be known according to the initial/boundary conditions.This can be considered as the main shortcoming of the algorithm.Because of the existence of nonlinear part in (3.1), it is not possible to exactly find the optimal value of Lagrange multiplier.Hence, it is necessary to consider a limitation on the non linear part causing this part to be ignored.Therefore, ũn is allocated to show the non linear part which has a special property.It is considered as a restricted variation [12,13], i.e. δũ n = 0. Firstly, we identify the Lagrange multiplier λ optimally via integration by parts.The successive approximations u n+1 (x) of the solution will be readily obtained upon using the obtained Lagrange multiplier and by using any selective function u 0 .After identifying the multiplier in (3.2), we have the following iteration algorithm Notice that after identifying the Lagrange multiplier we can construct the iteration formula instead of the iteration algorithm (3.3).We call (3.5) Variational Iteration Algorithm-II.This formulation is derived in [15].

Application of VIM
To solve the equation (1.1) with initial condition (1.2) by the VIM-II, the correction functional can be written as follows: By taking the variation on both sides of above equation with respect to y n and using the fact that δũ n = 0 we can deduce Via integration by parts, we find Integrating by parts (4.3) twice we obtain To find optimal value of λ, we impose the stationary condition δy n+1 (x) = 0 obtaining the following equations from which the Lagrangian multiplier can be identified as Next, by substituting (4.6) into (4.1), the desired iterative relation can be constructed as In order to start iteration using (4.7), y 0 (x) is needed.We represent it by Maclaurin series up the order three and, then, by substituting the given boundary conditions (1.2) into we conclude that From (4.7) and (4.8) with n = 0, we have the first VIM-iteration With the aid of Maple 18.0 package, we get and, clearly, yn for a fixed n ∈ N implementing an opportune do-loop.Here we point out that just the first approximation y1 is capable to be an efficient approximation of the exact solution as is the case of the Blasius equation, see [17].
Approximate Analytical Solution of a Third-Order IVP 123

Numerical results and comparisons
In this section, we compare the VIM solution with the numerical solution obtained by the following selected explicit Runge-Kutta methods: 1. DOPRI: the seven-stage fifth order Runge-Kutta method as in [6], 2. RK4: the four-stage fourth order Runge-Kutta method as in [4], 3. RKD5: the three-stage fifth order direct Runge-Kutta (RKD) method derived in [18].
To use Runge-Kutta methods, following [2], (1.1) need to be reduced to the following system of three first-order equations: In [18], the authors introduce a new method called RKD5 that requires less function evaluations than the RK4 and DOPRI methods, essentially because the reduced system 5.1 is three times the dimension.From [7,8,20] the exact solution of (1.1) is given in parametric form as where Ai(s) and Bi(s) denote the Airy functions and α, β and s 0 are constants to be determined.Imposing the initial conditions (1.2) we find that α = 0.676482, β = 1.60629, s 0 = −0.39685.
The plot for the exact solution is presented in Figure 1.We plot only the solution for y ≥ 0 since the case y < 0 is a drop with negative height and can be ignored because not physically interest.
From Table 1 we observe that the numerical results using DOPRI, RK4, RKD5 methods and VIM yield five-decimal-place accuracy and, further, Figure 4, obtained with GNU Octave 4.0.0,visualize the good agreement of those methods.In [20], the authors use the Adomian decomposition method (ADM) to obtain a power series approximation to the exact solution.Basically, the ADM provides an effective algorithm for approximating the nonlinear term y −2 in 5.1.Here, we consider two polynomials approximation from [20] Notice that as the number of iterations increases, the amount of calculations needed to obtain the new Adomian polynomial also increases and the process of obtaining   these polynomials becomes tedious and time consuming.Such drawback is eliminated through the application of He's variational iteration method, which does not require the additional computations of the Adomian polynomials.In Figure 2 we plot a comparison between y 6 , y 7 and the exact solution according to ( [20], Fig. 1, pp. .Notice that Theorem 1 in [20] holds also for our type of solutions and, then, as a byproduct, the analytical solution, the ADM solution and the VIM solution converge for x 2 < 1 and diverge for x 2 > 1.In Figure 3 we plot a comparison of first VIM iteration and the exact solution observing a very good agreement.
1.The natural first attempt to solve equation (1.1) by simply calling Maple's dsolve command, without imposing initial values conditions, just gives something like the an implicit formulation as the answer in terms of Airy functions, constants _C1, _C2, _C3 and using the RootOf command in terms of the global variable _Z.From here, applying the boundary conditions (1.2) seems quite difficult.
2. The exact solution (5.2) is plotted by using the function odeplot in the plots package with the ref ine = 1 option which tells odeplot to us all the stored points for the plot.The odeplot function plots solution curves obtained from the output of a call to dsolve is able to detect a singularity of the problem, in fact after run the command it returns the following warning message: cannot evaluate solution further left of −1.1175363, probably a singularity.We find, simply using fsolve command, that the considered solutions of (1.1), (1.2) have a zero at the point x as reported in Table 2.This points lie outside the domain of convergence x 2 < 1 and, thus, are considered singularities for the solutions.Notice that in [20] the authors, in order to include also singularity points the domain of convergence, truncate the obtained ADM solution at a suitable order to ensure that the truncated solution

Conclusions
In this paper, the VIM has been successfully applied to finding the solutions of (1.1) and (1.2).The obtained solutions are compared with numerical and ADM solutions.Those comparisons point out that the VIM may be considered as an efficient alternative method for solving a wide range of physical nonlinear problems, especially those arising in the area of fluid dynamics.We point out how the VIM can also be applied to any heat transfer and flow problem leading to a coupled nonlinear system of ODEs, see for instance [5,14,16,17,23,31,35] without any claim to completeness.The major advantage of the VIM, unlike the common numerical methods, is of providing analytical approximation or an approximated solution without linearization, perturbation, closure approximation, or discretization.Furthermore, unlike the ADM, VIM does not require the additional computations of the Adomian polynomials.

Equation ( 3 . 3 )
is called Variational Iteration Algorithm-I.Hence, according to iterations of this sequence, we can determine approximations of the exact solution as follow u(x) = lim n→∞ u n (x).(3.4)

Figure 3 :
Figure 3: Plot comparing the 1 st VIM iteration solution with the exact solution.

Table 1 :
Table comparing values of DOPRI, RK4 and RKD5 methods taking h = 0.01 as step size, the exact solution and the first VIM-iteration at x ∈ [0.0, 0.2, 0.4, 0.6, 0.8, 1.0].Analytical Solution of a Third-Order IVP127satisfies the so-called contact line condition.Based on the fact that (1.1) is autonomous it is possible to opportunely shift x → x + α to get a numerical solution in the interval x 2 < 1 subject to the initial conditions in α.

Table 2 :
Table comparing zero values for different classes of solutions.