Abstract
This article presents an improved spectralhomotopy analysis method (ISHAM) for solving nonlinear differential equations. The implementation of this new technique is shown by solving the FalknerSkan and magnetohydrodynamic boundary layer problems. The results obtained are compared to numerical solutions in the literature and MATLAB's bvp4c solver. The results show that the ISHAM converges faster and gives accurate results.
Keywords:
FalknerSkan flow; MHD flow; improved spectralhomotopy analysis methodIntroduction
Boundary layer flow problems have wide applications in fluid mechanics. In this article, we propose an improved spectralhomotopy analysis method (ISHAM) for solving general boundary layer problems. Three boundary layer problems are considered and solved in this study using the novel technique. The first problem considered is the classical twopoint nonlinear boundary value Blasius problem which models viscous fluid flow over a semiinfinite flat plate. Although solutions for this problem had been obtained as far back as 1908 by Blasius [1], the problem is still of great interest to many researchers as can be seen from the several recent studies [25].
The second problem considered in this article is the thirdorder nonlinear FalknerSkan equation. The FalknerSkan boundary layer equation has been studied by several researchers from as early as 1931 [6]. More recent studies of the solutions of the The FalknerSkan equation include those of Harries et al. [7], Pade [8] and Pantokratoras [9]. The third problem considered is magnetohydrodynamic (MHD) boundary layer flow. Such boundary layer problems arise in the study of the flow of electrically conducting fluids such as liquid metal. Owing to its many applications such as power generators, flow meters, and the cooling of reactors, MHD flow has been studied by many researchers, for example [10,11].
Owing to the nonlinearity of equations that describe most engineering and science phenomena, many authors traditionally resort to numerical methods such as finite difference methods [12], RungeKutta methods [13], finite element methods [14] and spectral methods [4] to solve the governing equations. However, in recent years, several analytical or semianalytical methods have been proposed and used to find solutions to most nonlinear equations. These methods include the Adomian decomposition method [1517], differential transform method [18], variational iteration method [19], homotopy analysis method (HAM) [2023], and the spectralhomotopy analysis (SHAM) (see Motsa et al. [24,25]) which sought to remove some of the perceived limitations of the HAM. More recently, successive linearization method [2628], has been used successfully to solve nonlinear equations that govern the flow of fluids in bounded domains.
In this article, boundary layer equations are solved using the ISHAM. The ISHAM is a modified version of the SHAM [24,25]. One strength of the SHAM is that it removes restrictions of the HAM such as the requirement for the solution to conform to the socalled rule of solution expression and the rule of coefficient ergodicity. Also, the SHAM inherits the strengths of the HAM, for example, it does not depend on the existence of a small parameter in the equation to be solved, it avoids discretization, and the solution obtained is in terms of an auxiliary parameter ħ which can conveniently be chosen to determine the convergence rate of the solution.
Mathematical formulation
We consider the general nonlinear thirdorder boundary value problem
subject to the boundary conditions
where c_{i}, b_{j }(i = 1, ..., 4 j = 1, 2, 3) are constants.
Equation 2.1 can be solved easily using methods such as the HAM and the SHAM. In each of these methods, an initial approximation f_{0}(η) is sought, which satisfies the boundary conditions. The speed of convergence of the method depends on whether f_{0}(η) is a good approximation of f (η) or not. The approach proposed here seeks to find an optimal initial approximation f_{0 }that would lead to faster convergence of the method to the true solution. We thus first seek to improve the initial approximation that is used later in the SHAM to solve the governing nonlinear equation.
We assume that the solution f(η) may be expanded as an infinite sum:
where f_{i}'s are unknown functions whose solutions are obtained using the SHAM at the ith iteration and f_{n}, (n ≥ 1) are known from previous iterations. The algorithm starts with the initial approximation f_{0}(η) which is chosen to satisfy the boundary conditions (2.2). An appropriate initial guess is
Substituting (2.3) in the governing equation (2.12.2) gives
subject to the boundary conditions
where the coefficient parameters a_{k,i1}, (k = 1, ..., 3) and r_{i1 }are defined as
Starting from the initial approximation (2.4), the subsequent solutions f_{i }(i ≥ 1) are obtained by recursively solving Equation 2.5 using the SHAM, [24,25]. To find the solutions of Equation 2.5, we begin by defining the following linear operator:
where q ∈ 0[1] is the embedding parameter, and F_{i}(η; q) is an unknown function.
The zerothorder deformation equation is given by
where ħ is the nonzero convergence controlling auxiliary parameter and is a nonlinear operator given by
Differentiating (2.10) m times with respect to q and then setting q = 0, and finally dividing the resulting equations by m! yield the mthorder deformation equations:
subject to the boundary conditions
where
The initial approximation f_{i,0 }that is used in the higherorder equations (2.12) is obtained on solving the linear part of Equation 2.5 which is given by
subject to the boundary conditions:
Since the coefficient parameters and the righthand side of Equation 2.15 for i = 1, 2, 3, ... are known (from previous iterations), the equation can easily be solved using numerical methods such as finite differences, finite elements, RungeKuttabased shooting methods or collocation methods. In this article, Equation 2.15 are solved using the Chebyshev spectral collocation method. The method (see, for example, [2931]), is based on the Chebyshev polynomials defined on the interval [1, 1] by
To implement the method, the physical region [0, ∞) is transformed into the region [1, 1] using the domain truncation technique whereby the problem is solved in the interval [0, L] instead of [0, ∞). This leads to the mapping
where L is the scaling parameter used to invoke the boundary condition at infinity. We use the popular GaussLobatto collocation points [29,31] to define the Chebyshev nodes in [1, 1], namely:
where N is the number of collocation points. The variable f_{i,0 }is approximated by the interpolating polynomial in terms of its values at each of the collocation points by employing the truncated Chebyshev series of the form:
where T_{k }is the kth Chebyshev polynomial. Derivatives of the variables at the collocation points may be represented by
where s is the order of differentiation and , with being the Chebyshev spectral differentiation matrix (see, for example [29,31]) whose entries are defined as
Substituting Equations 2.202.21 in 2.152.16 gives
subject to
where
In the above definitions, T stands for transpose and a_{k,i1}(k = 1, 2, 3) denotes a diagonal matrix of size (N + 1) × (N + 1). The boundary condition f_{i }(ξ_{N }) = 0 is implemented by deleting last row and last column of A_{i1}, and deleting the last rows of F_{i,0 }and R_{i1}. The derivative boundary conditions in (2.24) are then imposed on the resulting first row and last row of A_{i1 }and setting the first and last rows of F_{i,0 }and R_{i1 }to be zero. The solutions for f_{i.0}(ξ) are then obtained from soloving
In a similar manner, applying the Chebyshev spectral transformation on the higher order deformation equations (2.12)(2.13) gives
subject to the boundary conditions
where A_{i1 }and R_{i1}, are as defined in (2.25) and (2.27), respectively, and
To implement the boundary condition f_{i,m}(ξ_{N }) = 0, we delete the last rows of P_{i,m1 }and R_{i1 }and delete the last row and the last column of A_{i1 }in (2.29). The other boundary conditions in (2.30) are imposed on the first and the last rows of the modified A_{i1 }matrix on the left side of the equal sign in (2.29). The first and the last rows of the modified A_{i1 }matrix on the right side of the equal sign in (2.29) are then set to be zero. This results in the following recursive formula for m ≥ 1:
where Ã_{i1 }is the modified matrix A_{i1 }after incorporating the boundary conditions (2.30). Thus, starting from the initial approximation, which is obtained from (2.28), higherorder approximations f_{i,m}(ξ) for m ≥ 1, can be obtained through the recursive formula (2.33).
The solutions for f_{i }are then generated using the solutions for f_{i, m }as follows:
The [i, m] approximate solution for f (η) is then obtained by substituting f_{i }(obtained from 2.34) in equation 2.3.
Results and discussion
Table 1 shows the values of f" (0) at different orders [i, m] of the ISHAM approximation for the Blasius boundary layer flow when L = 30, ħ = 1 and N = 80. It is worth noting here that the numerical solution given by Howarth [32] is f" (0) = 0.332057, while the numerical result by the Matlab bvp4c routine is f" (0) = 0.33205734. Asaithambi [33] found this number correct to nine decimal positions as 0.332057336. It is evident that the ISHAM converges to the numerical result at orders [3,1] and [2,2]. Moreover, Table 1 shows that the ISHAM solution converges to the accurate solution of Howarth and the bvp4c result faster than the original SHAM results of which are those given in the first row of Table 1 (for the case when i = 1).
Table 1. Order [i, m] ISHAM approximate results for f" (0) of the Blasius boundary layer flow (Example 1) using L = 30, ħ = 1 and N = 80
In general, at order [i, m], i is the number of improvements of the initial approximation f_{0}(η) for f(η), and m is the number of improvements of the initial guess f_{q},_{0}(η); q = 1, 2, ..., i, for each application of the ISHAM. Table 2 gives a sense of the convergence rate of the ISHAM when compared with the numerical method for the Blasius problem at different values of η. In all the instances, convergence of the ISHAM is achieved at the second order.
Table 2. Comparison between the [m, m] ISHAM results and the bvp4c numerical results for the velocity pro le f' (η) at selected values of η for the Blasius boundary layer flow (Example 1) using L = 30, ħ = 1 and N = 200
Table 3 gives the values of f" (0) obtained used the ISHAM and the numerical method for various values of β for the FalknerSkan boundary layer problem. Full convergence is again achieved at order [2,2] for all the parameter values.
Table 3. Order [m, m] ISHAM approximate results for f" (0) of the FalknerSkan boundary layer flow (Example 2) using L = 30, ħ = 1 and N = 80
For the MHD boundary layer problem, Tables 4 and 5 illustrate the exact and approximate values of f' (η) and f" (0) at different values of η and the magnetic parameter M, respectively. The absolute errors in the approximations are also given. The tables show that the ISHAM converges rapidly with marginal or no errors after order [2,2].
Table 4. Order [m, m] ISHAM approximate results for the velocity profile f' (η) of the MHD boundary layer flow (Example 3) when M = 10 using L = 10, ħ = 1 and N = 200
Table 5. Order [m, m] ISHAM approximate results for f" (η) of the MHD boundary layer flow (Example 3) for different values of M using L = 10, ħ = 1 and N = 200
Conclusion
In this article, we have proposed an ISHAM for solving general nonlinear differential equations. This novel technique was compared against both numerical approximations and the MATLAB bvp4c routine for solving FalknerSkan and MHD boundary layer problems. The results demonstrate the relatively more rapid convergence of the ISHAM, and they show that the ISHAM is highly accurate.
Abbreviations
HAM: homotopy analysis method; ISHAM: improved spectralhomotopy analysis method; MHD: magnetohydrodynamic; SHAM: spectralhomotopy analysis.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
SSM developed the Matlab codes and generated the results. GTM and PS conceived of the study and formulated the problem. SS participated in the analysis of the results and manuscript coordination. All authors typed, read and approved the final manuscript.
Acknowledgements
The authors wish to acknowledge financial support from the University of Swaziland, University of KwaZuluNatal, University of Venda, and the National Research Foundation (NRF).
References

Blasius, H: Grenzschichten in Flussigkeiten mit kleiner Reibung. Z Math Phys. 56, 1–37 (1908)

Ahmad, F, AlBarakati, WH: An analytic solution of the Blasius problem. Commun Nonlinear Sci Numer Simul. 14, 1020–1024 (2009)

AlizadehPahlavan, A, BorjianBoroujeni, S: On the analytic solution of viscous fluid flow past a flat plate. Phys Lett A. 372, 3678–3682 (2008). Publisher Full Text

Parand, K, Shahini, M, Dehghan, M: Solution of a laminar boundary layer flow via a numerical method. Commun Nonlinear Sci Numer Simulat. 15, 360–367 (2010). Publisher Full Text

Yun, BI: Intuitive approach to the approximate analytical solution for the Blasius problem. Appl Math Comput. 208, 156–164 (2009). Publisher Full Text

Falkner, VM, Skan, SW: Some approximate solutions of the boundary layer equations. Philos Mag. 12, 865–896 (1931)

Harris, SD, Ingham, DB, Pop, I: Unsteady heat transfer in impulsive FalknerSkan flows: constant wall temperature case. Eur J Mech B. 21, 447–468 (2002). Publisher Full Text

Padé, O: On the solution of FalknerSkan equations. J Math Anal Appl. 285, 264–274 (2003). Publisher Full Text

Pantokratoras, A: The FalknerSkan flow with constant wall temperature and variable viscosity. Int J Thermal Sci. 45, 378–389 (2006). Publisher Full Text

Rashidi, MM: The modified differential transform method for solving MHD boundarylayer equations. Comput Phys Commun. 180, 2210–2217 (2009). Publisher Full Text

Parand, K, Rezai, AR, Ghaderi, SM: An approximate solution of the MHD FalknerSkan flow by Hermite functions pseudospectral method. Commun Nonlinear Sci Numer Simulat (2010)

Asaithambi, A: A secondorder finitedifference method for the FalknerSkan equation. Appl Math Comput. 156, 779–786 (2004). Publisher Full Text

Cortell, R: Numerical solutions of the classical Blasius flatplate problem. Appl Math Comput. 170, 706–710 (2005). Publisher Full Text

Asaithambi, A: Numerical solution of the FalknerSkan equation using piecewise linear functions. Appl Math Comput. 159, 267–273 (2004). Publisher Full Text

Elgazery, NS: Numerical solution for the FalknerSkan equation. Chaos Soliton Fract. 35, 738–746 (2008). Publisher Full Text

Wang, L: A new algorithm for solving classical Blasius equation. Appl Math Comput. 157, 1–9 (2004). Publisher Full Text

Abbasbandy, S: A numerical solution of Blasius equation by Adomian's decomposition method and comparison with homotopy perturbation method. Chaos Soliton Fract. 31, 257–260 (2007). Publisher Full Text

Kuo, B: Heat analysis for the FalknerSkan wedge flow by the differential transformation method. Int J Heat Mass Transfer. 48, 5036–5046 (2005). Publisher Full Text

Wazwaz, A: The variational iteration method for solving two forms of Blasius equation on a halfinfinite domain. Appl Math Comput. 188, 485–491 (2007). Publisher Full Text

Liao, SJ: Beyond Perturbation: Introduction to Homotopy Analysis Method. Chapman & Hall/CRC Press (2003)

Yao, B, Chen, J: A new analytical solution branch for the Blasius equation with a shrinking sheet. Appl Math Comput. 215, 1146–1153 (2009). Publisher Full Text

Yao, B, Chen, J: Series solution to the FalknerSkan equation with stretching boundary. Appl Math Comput. 215, 1146–1153 (2009). Publisher Full Text

Yao, B: Approximate analytical solution to the FalknerSkan wedge flow with the permeable wall of uniform suction. Commun Nonlinear Sci Numer Simulat. 14, 3320–3326 (2009). Publisher Full Text

Motsa, SS, Sibanda, P, Shateyi, S: A new spectralhomotopy analysis method for solving a nonlinear second order BVP. Commun Nonlinear Sci Numer Simulat. 15, 2293–2302 (2010). Publisher Full Text

Motsa, SS, Sibanda, P, Awad, FG, Shateyi, S: A new spectralhomotopy analysis method for the MHD JefferyHamel problem. Comput Fluids. 39, 1219–1225 (2010). Publisher Full Text

Makukula, Z, Motsa, SS, Sibanda, P: On a new solution for the viscoelastic squeezing flow between two parallel plates. J Adv Res Appl Math. 2(4), 31–38 (2010). Publisher Full Text

Makukula, ZG, Sibanda, P, Motsa, SS: A novel numerical technique for twodimensional laminar flow between two moving porous walls. Math Problems Eng. 15, Article ID 528956 (2010)

Shateyi, S, Motsa, SS: Variable viscosity on magnetohydrodynamic fluid flow and heat transfer over an unsteady stretching surface with Hall effect. Boundary Value Problems. 2010, 1–20 (2010)

Canuto, C, Hussaini, MY, Quarteroni, A, Zang, TA: Spectral Methods in Fluid Dynamics. SpringerVerlag, Berlin (1988)

Don, WS, Solomonoff, A: Accuracy and speed in computing the Chebyshev Collocation Derivative. SIAM J Sci Comput. 16(6), 1253–1268 (1995). Publisher Full Text

Howarth, L: On the solution of the laminar boundary layer equations. Proc R Soc Lond A. 164, 547–579 (1938). Publisher Full Text

Asaithambi, A: Solution of the FalkneSkan equation by recursive evaluation of Taylor coefficients. J Comput Appl Math. 176, 203–14 (2005). Publisher Full Text