An improved spectral homotopy analysis method ... - Semantic Scholar
Numer Algor (2012) 60:463–481 DOI 10.1007/s11075-011-9523-0 ORIGINAL PAPER
An improved spectral homotopy analysis method for MHD flow in a semi-porous channel Sandile Sydney Motsa · Stanford Shateyi · Gerald T. Marewo · Precious Sibanda
Received: 1 April 2011 / Accepted: 28 November 2011 / Published online: 20 December 2011 © Springer Science+Business Media, LLC 2011
Abstract In this paper we report on a novel method for solving systems of highly nonlinear differential equations by blending two recent semi-numerical techniques; the spectral homotopy analysis method and the successive linearisation method. The hybrid method converges rapidly and is an enhancement of the utility of the original spectral homotopy analysis method (Motsa et al., Commun Nonlinear Sci Numer Simul 15:2293–2302, 2010; Computer & Fluids 39:1219–1225, 2010) and an improvement on other recent semianalytical techniques. We illustrate the application of the method by solving a system of nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. A comparison with the numerical solution confirms the validity and accuracy of the technique and shows that the method converges rapidly and gives very accurate results.
S. S. Motsa · P. Sibanda School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, Republic of South Africa S. S. Motsa e-mail: [email protected] P. Sibanda e-mail: [email protected] S. Shateyi (B) Department of Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, Republic of South Africa e-mail: [email protected] G. T. Marewo Department of Mathematics, University of Swaziland, Private Bag 4, Kwaluseni, Swaziland e-mail: [email protected]
464
Numer Algor (2012) 60:463–481
Keywords Improved spectral homotopy analysis method · Coupled nonlinear equations · Laminar two-dimensional flow · Electrically conducting incompressible viscous fluid
1 Introduction Most engineering and physical systems are described by complex systems of nonlinear equations whose solution pose major challenges. The majority of these equations have no closed-form solutions, and although there already exist a plethora of numerical methods for finding approximate solutions to nonlinear equations, it has always been an interesting problem and a challenge to devise new and more efficient algorithms, particularly analytical or semianalytical techniques for finding solutions of nonlinear equations. Analytic techniques such as perturbation methods have been extensively used for many years, for example, the homotopy continuation method has been used since the 1930’s [53]. Traditionally, almost all perturbation methods depend on the ‘small parameter assumption’ that greatly constrained the application of such methods for the vast majority of nonlinear problems which do not involve such a small parameter [46]. Such methods are therefore valid for only nonlinear problems [34]. However new and more efficient techniques are constantly being developed and in the last three decades the delta-expansion method [26], the Adomian decomposition method [6, 7], the homotopy analysis method (HAM) [25, 28, 29, 31], the homotopy perturbation method (HPM) [23, 24], the variational iteration method (VIM) [21, 22] and the homotopy decomposition method [15] are some of the methods that have been developed and used extensively to solve problems of nonlinear mechanics. The Adomian technique provides a very simple but efficient numerical technique for solving large classes of nonlinear equations. The method gives the solution as an infinite series that should generally converge to an accurate solution. The main weakness of the method however is that it is often difficult to prove the convergence of the series solution [15]. The modified Adomian decomposition method developed by Abbasbandy [1] is fairly efficient in solving nonlinear problems although the formula is prone to diverge if a poor choice of the initial guess is made. A method combining the Adomian method and the homotopy continuation method was developed by Wu [53] to address the problem of divergence of Abbasbandy’s iteration. An alternative approach that combines the Newton-Raphson and the homotopy continuation method, the so-called Newton-homotopy continuation method was proposed by Wu [54]. In his paper some robust rules for the choice of the auxiliary homotopy function to avoid the problem of divergence of many traditional numerical methods were proposed. Other attempts at modification of the Adomian method can be found, for example, in [4, 5, 10, 11, 16]. The homotopy analysis method (HAM) was developed by Liao [25, 28] and has been used to solve nonlinear differential equations in different fields
Numer Algor (2012) 60:463–481
465
of science and engineering. Unlike the Adomian method and other nonperturbation techniques that cannot guarantee the convergence of the solution series, the HAM provides for a simple mechanism to control and adjust the convergence rate and the convergence region of the solution series by means of an embedded convergence controlling parameter [15, 31]. The homotopy analysis method provides certain rules for the construction of the solution such as the rule of solution expression and the rule of coefficient ergodicity. The former is necessary for the concept of homotopy while the rule of coefficient ergodicity is based on the completeness of base functions [30, 31]. Various attempts have been made to improve the performance of the HAM and the homotopy decomposition method is one such modification that sought to combine the Adomian and homotopy methods [15]. Further modifications of the HAM have been made by Bataineh et al. [12, 13] to obtain the socalled modified homotopy analysis method (MHAM). An optimal homotopyanalysis method (OHAM) has been used in a recent study by Liao [32]. The homotopy perturbation method (HPM) is a coupling of the traditional perturbation method and homotopy in topology and was proposed by JiHuan He [23, 24]. It has been extensively used to solve a wide class of functional equations including integro-differential equations, see for example [55]. The HPM is convenient, efficient and gives reasonable accuracy in solving nonlinear problems and integro-differential equations. However, as shown in [34, 35], the HPM, just like the VIM, is prone to give divergent approximations. The HPM was thus modified in [20] and the modified algorithm used in [8, 46] to test its validity and accuracy. Sajid and Hayat [48] explicitly proved that for the nonlinear problems where is different from −1, one cannot get convergent results by using HPM. They proved that the HPM results can be obtained as a special case of the HAM when = 1. Similar conclusions were independently made by Abbasbandy [2], as well as by Liang and Jeffrey [34]. It should be noted that HPM and VIM are not the same though they can sometimes produce the same results [34]. The idea behind the variational iteration method is the construction of a correction function by a general Lagrange multiplier which can be identified optimally via the variational theory Ganji et al. [19] while homotopy perturbation method uses an embedding parameter and was developed by combining the standard homotopy and perturbation methods. In light of the limitations of some current analytical and semi-analytical approaches, it is important that new, accurate and robust methods are developed. One possible approach is to combine analytical and numerical techniques to produce hybrid methods that do not suffer from the limitations of analytical approaches. One such approach is the spectral homotopy analysis method (SHAM) which combines the Chebyshev spectral collocation method and the homotopy analysis method. This method (see Motsa et al. [41]) seeks to address some of the limitations of the standard homotopy analysis method. This method has the advantage that it is more flexible than the standard HAM in that it allows for a wider choice of linear operators, admits a much wider
466
Numer Algor (2012) 60:463–481
range of admissible values and converges much faster than the standard homotopy analysis method, [36, 37, 41–43]. Abbasbandy and Shivanian [3] presented a procedure based on pseudo-spectral method to calculate multiple (dual) solutions of a model of mixed convection in a porous medium with boundary conditions on a semi-infinite interval which admit dual solutions. In this paper we present a new and improved spectral homotopy analysis method (ISHAM) for solving highly nonlinear and coupled systems of ordinary differential equations. We use this method to solve a system of highly nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. The improved spectral homotopy analysis method is a hybrid method that blends the spectral homotopy analysis method and the successive linearisation method (SLM). The SLM is (see [9, 38, 39, 44, 45, 50]) a recent method that has been successfully utilized in solving several nonlinear problems. In the present study we show that the ISHAM is an enhancement of the original SHAM with better accuracy and rapid convergence.
2 Governing equations In this section we solve a test problem, the coupled nonlinear equations that governs the laminar two-dimensional flow of an electrically conducting incompressible viscous fluid between two parallel plates subject to a transverse magnetic field. We demonstrate how the hybrid spectral-homotopy analysis method can be applied to find solutions of the governing coupled second- and fourth-order nonlinear equations. The problem considered is described more fully in [17, 47, 52, 56]. The governing dimensionless equations are; ∂u ∂v + = 0, ∂x ∂y
(2.1)
2 ∂u v ∂u ∂ 2u Ha2 2 2 ∂ Py 2∂ u +v = −ε + ε u, + − u ∂x ∂y ∂x hq ∂ x2 ∂ y2 Re 2 ∂ Py ∂v v ∂v ∂ 2u 2∂ u +v =− + ε u + 2 , ∂x ∂y ∂x hq ∂ x2 ∂y
(2.2) (2.3)
where h is the distance between the two plates, Ha is the Hartmann number, Re is the Reynolds number, q is the normal velocity at the porous wall and ε is a small aspect ratio. Using the Berman similarity transformation v = −V(y),
u = u0 U(y) + x
dV , dy
it can easily be shown, see [56], that equations (2.1)–(2.2) reduce to the coupled system U − Ha2 U + Re(VU − U V ) = 0, V
IV
− Ha V + Re(VV − V V ) = 0, 2
(2.4) (2.5)
Numer Algor (2012) 60:463–481
467
subject to the boundary conditions U = 1, V = 0, V = 0, y = 0,
(2.6)
U = 0, V = 1, V = 0, y = 1.
(2.7)
Equations (2.4)–(2.5) were solved recently by Ziabakhsh and Domairry [56] using the homotopy analysis method. In this paper the equations are used as a test problem to show the capability and prove the validity and accuracy of the ISHAM. The results are compared to the numerical solutions reported in [56] and the results of the standard SHAM and SLM.
3 Improved spectral homotopy analysis method In this section we describe the use of the improved spectral homotopy analysis method (ISHAM) in the governing equations (2.4)–(2.7). To apply the ISHAM we begin with the SLM idea, where it is assumed that the solutions U(y) and V(y) can be expanded as U(y) = U i (y) +
i−1
U n (y),
V(y) = Vi (y) +
n=0
i−1
Vn (y), i = 1, 2, 3, . . . (3.8)
n=0
where U i and Vi are unknown functions whose solutions are obtained using the SHAM approach [36, 37, 41–43] at the ith iteration and U n , Vn (1 ≤ n ≤ i − 1) are known from previous iterations. The algorithm starts with the initial approximations U 0 (y) and V0 (y) which are chosen to satisfy the boundary conditions (2.6)–(2.7). We choose the following initial guesses U 0 (y) = 1 − y,
V0 (y) = 3y2 − 2y3 .
(3.9)
which satisfy the boundary conditions (2.6)–(2.7). Substituting (3.8) in the governing equations (2.4)–(2.7) gives (1) U i + d1,i−1 U i + d2,i−1 U i + d3,i−1 Vi + d4,i−1 Vi + Re(Vi U i − U i Vi ) = si−1 ,
(3.10) (2) , Viiv + d1,i−1 Vi + d2,i−1 Vi + e1,i−1 Vi + e2,i−1 Vi + Re(Vi Vi − Vi Vi ) = si−1
(3.11) subject to boundary conditions U i (0) = Vi (0) = Vi (0) = U i (1) = Vi (1) = Vi (1) = 0,
(3.12)
468
Numer Algor (2012) 60:463–481
where d1,i−1 = Re
i−1
Vn ,
d2,i−1 = −Ha2 − Re
n=0
d3,i−1 = −Re
Un,
d4,i−1 = Re
n=0
(1) si−1
=−
U n
− Ha
2
i−1
U n + Re
=−
i−1
Vn ,
e2,i−1 = Re
Vniv
i−1 n=0
− Ha
n=0
2
i−1
Vn
n=0
n=0
(2) si−1
U n ,
n=0
i−1
i−1
i−1 n=0
n=0
e1,i−1 = −Re
Vn ,
n=0
i−1
i−1
i−1
Vn
+ Re
n=0
i−1
U n
−
n=0
i−1
Un
i−1
n=0
Vn
.
n=0
Vn i−1 n=0
Vn
i−1 n=0
Vn
−
i−1
Vn
n=0
i−1
Vn
n=0
The successive linearisation method (see [9, 38, 39, 44, 45, 50]) is based on iteratively solving the linear components of equations (3.10) and (3.11) which are given as (1) U i + d1,i−1 U i + d2,i−1 U i + d3,i−1 Vi + d4,i−1 Vi = si−1 ,
(3.13)
(2) Viiv + d1,i−1 Vi + d2,i−1 Vi + e1,i−1 Vi + e2,i−1 Vi = si−1 ,
(3.14)
subject to boundary conditions 3.12. Equations (3.13)–(3.14) can be inverted using any numerical method of solutions. In some cases, it may be possible to solve the equations using analytical means. In the case of the ISHAM, an attempt is made to solve the full set of equations (3.10)–(3.11) without linearizing. Starting from the initial approximations U 0 (y), V0 (y) the subsequent solutions U i , Vi (i ≥ 1) are obtained by recursively solving equations (3.10)–(3.12) using the SHAM [36, 37, 41–43]. To find the SHAM solutions of (3.10)–(3.12) we begin by defining the following linear operators L1 [U˜ i (y; q), V˜ i (y; q)] =
∂ 2 U˜ i ∂ y2
˜
˜
+ d1,i−1 ∂∂Uy i + d2,i−1 U˜ i + d3,i−1 ∂∂Vy i + d4,i−1 V˜ i , (3.15)
L2 [V˜ i (y; q)] =
∂ 4 V˜ i ∂ y4
3 ˜ 2 ˜ ˜ + d1,i−1 ∂∂ yV3 i + d2,i−1 ∂∂ yV2 i + e1,i−1 ∂∂Vy i + e2,i−1 V˜ i ,
(3.16) where q ∈ [0, 1] is an embedding parameter and U˜ i (y; q), V˜ i (y; q) are unknown functions. If we also define auxiliary nonlinear operators N1 [U˜ i (y; q), V˜ i (y; q)] = L1 [U˜ i (y; q), V˜ i (y; q)] + Re(V˜ i U˜ i − U˜ i V˜ i ),
(3.17)
Numer Algor (2012) 60:463–481
469
N2 [V˜ i (y; q)] = L2 [V˜ i (y; q)] + Re(V˜ i V˜ i − V˜ i V˜ i ),
(3.18)
then the zero order deformation equations for this problem are
(1 − q)L1 [U˜ i (y; q), V˜ i (y; q)] − [U i,0 (y), Vi,0 (y)]
(1) = q N1 [U˜ i (y; q), V˜ i (y; q)] − si−1 ,
(3.19)
(2) , (1 − q)L2 [V˜ i (y; q) − Vi,0 (y)] = q N2 [V˜ i (y; q)] − si−1
(3.20)
where U i,0 (y), Vi,0 (y) are initial approximations for the unknown functions U(y), V(y) respectively and = 0 is the converging controlling auxiliary parameter. The corresponding m-th order deformation equations are L1 [U i,m , Vi,m ] − χm [U i,m−1 , Vi,m−1 ] = Sm,1 ,
(3.21)
L2 [Vi,m − χm Vi,m−1 ] = Sm,2 ,
(3.22)
subject to boundary conditions (0) = U i,m (1) = Vi,m (1) = Vi,m (1) = 0. U i (0) = Vi,m (0) = Vi,m
where
χm =
0, 1,
⎛ Sm,1 = L1 [U i,m−1 , Vi,m−1 ] + Re ⎝
m−1
m≤1 m>1
−(1 − and
(3.24)
Vi, jU i,m−1− j−
j=0
m−1
⎞ ⎠ U i, j Vi,m−1− j
j=0
(1) , χm )si−1
⎛
Sm,2 = L2 [Vi,m−1 ] + Re ⎝
(3.23)
(3.25)
m−1 j=0
Vi, j Vi,m−1− j−
m−1
⎞ (2) ⎠ Vi, j Vi,m−1− j − (1 − χm )si−1 ,
j=0
(3.26) The unknown functions U i,0 (y), Vi,0 (y) are obtained as solutions to the boundary value problem (1) + d1,i−1 U i,0 + d2,i−1 U i,0 + d3,i−1 Vi,0 + d4,i−1 Vi,0 = si−1 , U i,0
(3.27)
(2) iv Vi,0 + d1,i−1 Vi,0 + d2,i−1 Vi,0 + e1,i−1 Vi,0 + e2,i−1 Vi,0 = si−1 ,
(3.28)
470
Numer Algor (2012) 60:463–481
subject to boundary conditions U i,0 (0) = Vi,0 (0) = Vi,0 (0) = U i,0 (1) = Vi,0 (1) = Vi,0 (1) = 0.
(3.29)
We use the Chebyshev pseudospectral method to solve equations (3.27)– (3.29). Before using the Chebyshev pseudospectral method to solve this problem it is convenient to first use the mapping η = 2y − 1; 0 ≤ y ≤ 1 to transform the problem domain from [0, 1] into the domain [−1, 1] on which the Chebyshev pseudospectral method is implemented. The unknown functions U i,0 (η) and Vi,0 (η) are approximated as truncated series of Chebyshev polynomials of the form U i,0 (η) =
N
U i,0 (ηk )Tk (η j) and Vi,0 (η) =
k=0
N
Vi,0 (ηk )Tk (η j),
(3.30)
k=0
where j = 0, 1, . . . , N, Tk is the kth Chebyshev polynomial, and η0 , η1 , . . . , η N are Gauss-Lobatto collocation points (see [14]) defined by η j = cos
πj , N
j = 0, 1, . . . , N.
(3.31)
Derivatives of the functions U i,0 (y) and Vi,0 (η) at the collocation points are represented as ds U i,0 s = 2 Dsjk U i,0 (ηk ), dys N
k=0
ds Vi,0 s = 2 Dsjk Vi,0 (ηk ), dys N
(3.32)
k=0
where s is the order of differentiation and D is the Chebyshev spectral differentiation matrix (see for example [14, 51]). Substituting equations (3.30)– (3.32) in (3.27)–(3.29) yields
B11 B12 O B22
(1) Ui,0 S = i−1 (2) Vi,0 Si−1
(3.33)
subject to boundary conditions
U i,0 (η N ) = Vi,0 (η N ) =
N
D Nk Vi,0 (η N ) = U i,0 (η0 )
k=0
= Vi,0 (η0 ) =
N k=0
D0k Vi,0 (η0 ) = 0,
(3.34)
Numer Algor (2012) 60:463–481
471
where B11 = 24 D2 + 2d1,i−1 D + d2,i−1 , B12 = 2d3,i−1 D + d4,i−1 , B22 = 24 D4 + 23 d1,i−1 D3 + 22 d2,i−1 D2 + 2e1,i−1 D + e2,i−1 , Ui,0 = [U i,0 (η0 ), U i,0 (η1 ), . . . , U i,0 (η N )]T , Vi,0 = [Vi,0 (η0 ), Vi,0 (η1 ), . . . , Vi,0 (η N )]T , (1) (1) (1) (1) Si−1 = [Si−1 (η0 ), Si−1 (η1 ), . . . , Si−1 (η N )]T , (2) (2) (2) (2) Si−1 = [Si−1 (η0 ), Si−1 (η1 ), . . . , Si−1 (η N )]T , (2) dr,i−1 = diag(dr,i−1 (η1 ), dr,i−1 (η2 ), . . . , dr,i−1 (η N )); r = 1, 2, 3, 4, (2) er,i−1 = diag(er,i−1 (η1 ), er,i−1 (η2 ), . . . , er,i−1 (η N )); r = 1, 2,
and O is an (N + 1) × (N + 1) zero matrix. Once boundary conditions (3.34) are imposed on (3.33), the linear system is solved for Ui,0 , Vi,0 to reveal U i,0 (y), Vi,0 (y). Applying the Chebyshev spectral method in a similar manner to equations (3.21)–(3.23) gives B11 B12 Ui,m B11 B12 Ui,m−1 = (χm + ) O B22 Vi,m O B22 Vi,m−1 (1) S P (3.35) + i,m−1 − (1 − χm ) i−1 (2) Qi,m−1 Si−1 subject to boundary conditions U i,m (η N ) = Vi,m (η N ) =
N
D Nk Vi,m (η N ) = U i,m (η0 )
k=0
= Vi,m (η0 ) =
N
D0k Vi,m (η0 ) = 0.
(3.36)
k=0
where Ui,m = [U i,m (η0 ), U i,m (η1 ), . . . , U i,m (η N )]T , Vi,m = [Vi,m (η0 ), Vi,m (η1 ), . . . , Vi,m (η N )]T , ⎛ ⎞ m−1 m−1 Pi,m−1 = Re ⎝ Vi, jDUi,m−1− j − Ui, jDVi,m−1− j⎠ , ⎛ Qi,m−1 = Re ⎝
j=0 m−1 j=0
j=0
Vi, jD3 Vi,m−1− j −
m−1 j=0
⎞ DVi, jD2 Vi,m−1− j⎠ .
472
Numer Algor (2012) 60:463–481
Once boundary conditions are imposed, the solution to equations (3.35)–(3.36) is obtained as ˜ i,m−1 + B−1 [Ci,m−1 − (1 − χm )Si,m−1 ]; m = 1, 2, . . . , Zi,m = (χm + )B−1 BZ (3.37) ˜ are coefficient matrices in (3.35) modified to incorporate where B and B boundary conditions, and
Zi,m
(1) Ui,m Pi,m−1 S = , Ci,m−1 = , Si−1 = i−1 (2) Vi,m Qi,m−1 Si−1
(3.38)
Thus, starting from the initial approximation, which is obtained from (3.33), higher order approximations U i,m (y), Vi,m (η), for m ≥ 1, can be obtained through the recursive formula (3.37). The solutions for U i (y), Vi (y) are then generated using the solutions for U i,m , Vi,m as follows U i = U i,0 + U i,1 + U i,2 + U i,3 + · · · + U i,m ,
(3.39)
Vi = Vi,0 + Vi,1 + Vi,2 + Vi,3 + · · · + Vi,m .
(3.40)
The [i, m] approximate solutions for U(y) and V(y) are then obtained by substituting U i and Vi (obtained from (3.39)–(3.40)) in (3.8). We remark that the recursive scheme (3.37) reduces to the original spectral homotopy analysis method [36, 37, 41–43] when i = 1 and the successive linearisation method [9, 38, 39, 44, 45, 50] when m = 0.
4 Results and discussion The governing equations (2.4)–(2.7) were solved using the ISHAM approach described in the previous section. In this section we present the results showing the effect of the governing parameters on the important properties of the flow. To check the accuracy of the proposed ISHAM, comparison was made with numerical solutions obtained using the MATLAB routine bvp4c [27, 49], which is an adaptive Lobatto quadrature scheme. In computing all the results presented in this work the number of collocation points used was N = 60. This value N was found to be appropriate for all the physical parameters used in this section and gave results that matched the results of the MATLAB routine bvp4c when the specified relative tolerance, RelTol, and the specified absolute tolerance, AbsTol were both less than e−9 . We note that the iterative scheme (3.37) that generates the approximate solutions U(y), V(y) depends on the convergence controlling auxiliary parameter . Hence, a careful selection of must be made. In the case of the original homotopy analysis method (HAM), it was pointed out by Liao [28],
Numer Algor (2012) 60:463–481
473 14.041 m=2 m=4 m = 10 m = 16
14.0405 m=2 m=4 m = 10 m = 16
14.04 14.0395 14.039 14.0385 14.038
0
14.0375
0
Fig. 1 −curve for U (0) and V (0) when Re = 10, and Ha = 0
that the accuracy and convergence of the HAM series solution depends on the careful selection of the auxiliary parameter . It was suggested in [28] that the valid values of where the series converges correspond to the horizontal segment of each -curve, which is a plot of some property of the independent variable (or it’s derivatives) at a certain point in it’s governing domain. In this work, the optimal curve was identified as the value of at the critical point of the second order −curves generated using U (0) and V (0). Using numerical experimentation with values taken along the horizontal segment of the curves, it was found that the critical values of the second order curve resulted in the most convergent series. An sample curve is shown in Fig. 1. Tables 1 and 2 give a comparison of the SHAM (order [1, m] of the ISHAM), the homotopy analysis method (HAM) results reported in [56], and numerical results obtained using the Matlab in-built boundary value problem solver bvp4c. In particular, Tables 1 and 2 give a snap shot of the mean velocity values at different locations within the channel. We note that these
Table 1 Comparison of the numerical results against the SHAM (order [1, m] on the ISHAM) approximate solutions for U(y) and V(y) when Ha = 0 and Re = 1 with N = 60 and = −1 y U(y) 0.00 0.25 0.50 0.75 1.00 V(y) 0.00 0.25 0.50 0.75 1.00
1st order
2nd order
3rd order
bvp4c
Ref.[56]
1.000000000 0.661815551 0.380389261 0.164599897 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815298 0.380388853 0.164599643 0.000000000
0.000000000 0.164565728 0.514940157 0.852203849 1.000000000
0.000000000 0.164565723 0.514940146 0.852203844 1.000000000
0.000000000 0.164565723 0.514940146 0.852203844 1.000000000
0.000000000 0.164565723 0.514940146 0.852203845 1.000000000
0.000000000 0.164565724 0.514940140 0.852203842 1.000000000
474
Numer Algor (2012) 60:463–481
Table 2 Comparison of the numerical results, the SHAM (order [1, m] on the ISHAM) solutions for U(y) and V(y) when Re = 1, Ha = 1 and = −1 y U(y) 0.00 0.25 0.50 0.75 1.00 V(y) 0.00 0.25 0.50 0.75 1.00
1st order
2nd order
3rd order
bvp4c
Ref. [56]
1.000000000 0.624376612 0.342902417 0.144030095 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376291 0.342901931 0.144029842 0.000000000
0.000000000 0.165179106 0.514478084 0.851140074 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179101 0.514478095 0.851140064 1.000000000
results rapidly converge to the bvp4c results, with full convergence achieved at order 2. This is a good indication of the improved accuracy of the new solution approach. It is also worth mentioning here that, in their paper, Ziabakhsh and Domairry [56] did not disclose the order of approximation used to generate the HAM results (although we estimate it cannot be anything less than the 10th order) nor do they reveal the value of used. Reproducing their results for purposes of comparing with the present approximations is therefore impractical. Nonetheless, their approximations do not yield the same level of accuracy as the current method. Tables 3 and 4 demonstrate the effect of the variation of Reynolds number on U (0) and V (0), respectively, when the Hartman number is fixed. The tables give a comparison of the SHAM ([1, m]), ISHAM ([i, m] i > 1) and Table 3 ISHAM results for U (0) when Ha = 1 using various values of Re Re
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
10
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.91 −0.91 −0.91 −0.91 −0.91 (SLM →) −0.68 −0.68 −0.68 −0.68 −0.68 (SLM →) −0.55 −0.55 −0.55 −0.55 −0.55
−3.67855093 −3.57357902 −3.57381695 −3.57381721 −3.57381721 −3.57381721 −9.09456068 −7.49394063 −7.48966570 −7.48908565 −7.48908410 −7.48908409 −14.45969153 −10.51500928 −10.46272106 −10.45596459 −10.45586600 −10.45586393
−3.57397714 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.53495863 −7.48908437 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.68732625 −10.45589567 −10.45586397 −10.45586388 −10.45586388 −10.45586388
−3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.48910428 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.45635538 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388
−3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388
50
100
Numer Algor (2012) 60:463–481
475
Table 4 ISHAM results for V (0) when Ha = 1 using various values of Re Re
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
10
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.91 −0.91 −0.91 −0.91 −0.91 (SLM →) −0.68 −0.68 −0.68 −0.68 −0.68 (SLM →) −0.55 −0.55 −0.55 −0.55 −0.55
10.01969746 10.06518626 10.06486048 10.06485955 10.06485955 10.06485955 19.10301189 19.20789938 19.18030717 19.17886834 19.17885842 19.17885835 26.33729449 26.34830439 26.24922129 26.23538456 26.23501735 26.23500654
10.06478254 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.16893291 19.17885931 19.17885836 19.17885836 19.17885836 19.17885836 26.20696262 26.23503458 26.23500621 26.23500616 26.23500616 26.23500616
10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.17885554 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 26.23496043 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616
10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616
50
100
SLM ([i, 0]) generated results. The results presented in the tables were verified against the bvp4c numerical results. It can be seen that for large Reynolds numbers more terms in the SHAM solution series ar e required to achieve a match with the numerical results than in the case of the SLM. This shows that the SLM converges much faster than the SHAM. Essentially, the ISHAM results are a combination between and i steps of the SLM and m steps of the SHAM. The computing times required to generate the results are of the order of fractions of one second for all three methods. Hence, meaningful Table 5 ISHAM results for U (0) when Re = 10 using various values of Ha Ha 0
5
10
Order (m) 0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.90 −0.90 −0.90 −0.90 −0.90 (SLM →) −0.95 −0.95 −0.95 −0.95 −0.95 (SLM →) −0.97 −0.97 −0.97 −0.97 −0.97
[1,m](SHAM) −3.56623375 −3.46024971 −3.46050840 −3.46050933 −3.46050933 −3.46050933 −5.90262607 −5.82294691 −5.82300329 −5.82300332 −5.82300332 −5.82300332 −10.39354751 −10.33060703 −10.33062021 −10.33062025 −10.33062025 −10.33062025
[2,m] −3.46067933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82305801 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33063143 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
[3,m] −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
[4,m] −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
476
Numer Algor (2012) 60:463–481
Table 6 ISHAM results for V (0) when Re = 10 using various values of Ha Ha
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
0
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.90 −0.90 −0.90 −0.90 −0.90 (SLM →) −0.95 −0.95 −0.95 −0.95 −0.95 (SLM →) −0.97 −0.97 −0.97 −0.97 −0.97
9.98779654 10.03712703 10.03677402 10.03677181 10.03677181 10.03677181 10.94299138 10.94069274 10.94062918 10.94062917 10.94062917 10.94062917 14.04090221 14.03847281 14.03848032 14.03848033 14.03848033 14.03848033
10.03668989 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94059825 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03847155 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
5
10
comparison between the three methods is done using the total number of iteration steps required to achieve full convergence to the numerical solutions. For the SHAM and SLM approaches the total number of iteration steps is equivalent to their respective orders m and i. Since for i iterations there are m SHAM steps of the ISHAM, it follows that the total number of iterations of the ISHAM is i × m. Tables 3 and 4 show that the SLM converges much faster than the SHAM for all values of Re. Comparing the SLM and ISHAM results, it can be noted that the SLM converges faster than the ISHAM. For example, when Re = 100 we note from Table 3 that at the ISHAM order [2, 4] (total of 8 iterations),
1
1
0.9
Ha = 0 Ha = 5 Ha = 20
0.8
0.9
Ha = 0 Ha = 5 Ha = 20
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2 0.2 0.1 0.1
0 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 2 The comparison between the numerical results and [2,2] order ISHAM solution for U(y), V(y) and V (y), for different values of Ha, when Re = 1, N = 60, = −1
Numer Algor (2012) 60:463–481
477 1
1 Re = 1 Re = 10 Re = 20
0.9
Re = 1 Re = 10 Re = 20
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 3 The comparison between the numerical results and [2,2] order ISHAM solution for U(y), V(y) and V (y), for different values of Re, when Ha = 1, N = 60, = −1
U (0) is −10.45586397 compared to the converged solution of −10.45586388 which is obtained after only 4 iterations of the SLM. The same observation can be made by considering Table 4 for V (0). By the same argument, it can be seen that the ISHAM converges more rapidly than then SHAM approach. Tables 5 and 6 demonstrate the effect of the variation of Hartman number on U (0) and V (0), respectively, when the Reynolds number is fixed. Increasing the Hartmann number has no noticeable effect on the convergence rate of the SHAM. Again, it can be seen from Tables 5 and 6 that the SLM converges much faster than ISHAM which in turn converges much faster than the SHAM. Figures 1, 2, 3 and 4 give a comparison between the [2,2] order ISHAM approximations and the numerical results for different Hartmann numbers and increasing Reynolds numbers. The parameter values used are the same as those used by Ziabakhsh and Domairry [56] in order to facilitate a comparison
1.6
1.6 Ha = 0 Ha = 5 Ha = 20
1.4
Re = 1 Re = 10 Re = 20
1.4
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
Fig. 4 Comparison between the numerical results and [2,2] order ISHAM solution for V (y) when Re = 1, Ha = 1 and N = 60, = −1
478
Numer Algor (2012) 60:463–481
of the results. We remark here that reasonable accuracy and matching of the numerical results and the ISHAM approximations is achieved at the lower orders of approximation of the SHAM, SLM and ISHAM. When Ha = 0 and Re = 1 the results presented here match exactly with the HAM results given by Ziabakhsh and Domairry [56]. As expected, increasing the Hartmann number has the effect of suppressing the flow of the fluid due to a resistive Lorentz force now slowing down the motion of the fluid.
5 Conclusion In this study, a new numerical perturbation scheme for solving complex nonlinear boundary value problems is presented. The method is based on a blend between the spectral homotopy analysis method (SHAM) and the successive linearisation method (SLM). The implementation of the improved spectral homotopy analysis method (ISHAM), was demonstrated by solving a system of nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. Validation of the results obtained in this study was done against numerical results obtained using the bvp4c Matlab routine. The numerical experimentation conducted in this study suggests that the proposed method converges much more rapidly to the numerical when compared to the original SHAM. However, comparison of the convergence rates between the ISHAM and the SLM for this particular set of equations showed that the SLM convergence rate much more rapidly than the ISHAM. This preliminary study shows that the ISHAM is an enhancement and improvement of the accuracy and utility of the SHAM. Consequently, the proposed method would be expected to perform better, with more rapid convergence and better accuracy than the standard HAM and other recent semi-analytical methods when solving nonlinear boundary value problems that arise in science and engineering. Finally, it is worth noting that although in this study we have used the socalled -curves to find the value of the convergence-control parameter , it might be more worthwhile to explore the use of the so-called optimal HAM, Liao [32]. Acknowledgements The authors acknowledge the financial support from the National Research Foundation (NRF). We wish to thank the anonymous reviewers for their helpful comments and suggestions.
References 1. Abbasbandy, S.: Improving Newton-Raphson method for nonlinear equations by modified Adomian decomposition method. Appl. Math. Comput. 145, 887–893 (2003) 2. Abbasbandy, S.: The application of homotopy analysis method to solve a generalized HirotaSatsuma coupled KdV equation. Phys. Lett. A 361, 478–483 (2007)
Numer Algor (2012) 60:463–481
479
3. Abbasbandy, S., Shivanian, E.: Multiple solutions of mixed convection in a porous medium on semi-infinte interval using pseudo-spectral collocation method. Commun. Nonlinear Sci. Numer. Simul. 16, 2745–2752 (2011) 4. Abbaoui, K., Cherruault, Y.: Convergence of Adomian’s method applied to non-linear equations. Math. Comput. Model. 20, 69–73 (1994) 5. Abdulaziz, O., Hashim, I., Saif, A.: Series solutions of time-fractional PDEs by homotopy analysis method, differential equations and nonlinear mechanics, vol. 2008, Article ID 686512. doi:10.1155/2008/686512 (2008) 6. Adomian, G.: Nonlinear stochastic differential equations. J Math. Anal. Appl. 55, 441–452 (1976) 7. Adomian, G.: A review of the decomposition method and some recent results for nonlinear equations. Comput. Math. Appl. 21, 101–127 (1991) 8. Rafiq, A., Javeria, A.: New iterative methods for solving nonlinear equations by using the modified homotopy perturbation method. Acta Univ. Apulensis 18, 129–137 (2009) 9. Awad, F.G., Sibanda, P., Motsa, S.S., Makinde, O.D.: Convection from an inverted cone in a porous medium with cross-diffusion effects. Comput. Math. Appl. 61, 1431–1441 (2011) 10. Babolian, E., Biazar, J.: Solution of non-linear equations by modified Adomian decomposition method. Appl. Math. Comput. 132, 167–172 (2002) 11. Babolian, E., Biazar, J., Vahidi, A.R.: On the decomposition method for system of linear equations and system of linear Volterra integral equations. Appl. Math. Comput. 147, 19–27 (2004) 12. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: Modified homotopy analysis method for solving systems of second-order BVPs. Commun. Nonlinear Sci. Numer. Simul. 14, 430–442 (2009) 13. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: On a new reliable modification of homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 14, 409–423 (2009) 14. Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods in Fluid Dynamics. Springer, Berlin (1988) 15. Chen, X., L, Zheng, Zhang, X.: Convergence of the homotopy decomposition method for solving nonlinear equations. Adv. Dyn. Syst. Appl. 2(1), 59–64 (2007) 16. Cherruault, Y.: Convergence of Adomian’s method. Kybernetes 18, 31–38 (1989) 17. Desseaux, A.: Influence of a magnetic field over a laminar viscous flow in a semi-porous channel. Int. J. Eng. Sci. 37, 1781–1794 (1999) 18. Don, W.S., Solomonoff, A.: Accuracy and speed in computing the Chebyshev Collocation Derivative. SIAM J. Sci. Comput. 16(6), 1253–1268 (1995) 19. Ganji, D.D., Nourollahi, M., Rostamian, M.: A comparison of variational iteration method with Adomian’s decompostion method in some highly nonlinear equations. Int. J. Sci. Technol. 2(2), 179–188 (2007) 20. Golbabaia, A., Javidi, M.: Application of homotopy perturbation method for solving eighthorder boundary value problems. Appl. Math. Comput. 191, 334–346 (2007) 21. He, J.H.: Variational iteration method: a kind of non-linear analytical technique: some examples. Int. J. Non-Linear Mech. 34(4), 699–708 (1999) 22. He, J.H.: Variational iteration method for autonomous ordinary differential systems. Appl. Math. Comput. 114(2–3), 115–123 (2000) 23. He, J.H.: Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 178, 257– 262 (1999) 24. He, J.H.: A coupling method of a homotopy technique and a perturbation technique for nonlinear problems. Int. J. Non-linear Mech. 35, 37–43 (2000) 25. Liao, S.J.: The proposed homotopy analysis technique for the solution of nonlinear problems. PhD thesis, Shanghai Jiao Tong University (1992) 26. Karmishin, A.M., Zhukov, A.I., Kolsov, V.G.: Methods of Dynamics Calculation and Testing for Thin-walled Structures. Mashinostroyenie, Moscow (1990) 27. Kierzenka, J., Shampine, L.: A BVP solver based on residual control and the Matlab PSE. ACM Trans. Math. Softw. 27, 299–316 (2001) 28. Liao, S.J.: Beyond Perturbation: Introduction to Homotopy Analysis Method. Chapman & Hall/CRC Press (2003) 29. Liao, S.J.: Comparison between the homotopy analysis method and the homotopy perturbation method. Appl. Math. Comput. 169, 1186–1194 (2005)
480
Numer Algor (2012) 60:463–481
30. Liao, S.J., Tan, Y.: A general approach to obtain series solutions of nonlinear differential equations. Stud. Appl. Math. 119, 297–354 (2007) 31. Liao, S.J.: Notes on the homotopy analysis method: Some definitions and theores. Commun. Nonlinear Sci. Numer. Simul. 14, 983–997 (2009) 32. Liao, S.J.: An optimal homotopy-analysis approach for strongly nonlinear differential equations. Commun. Nonlinear Sci. Numer. Simul. 15, 2003–2016 (2010) 33. Li, J.L.: Adomian’s decomposition method and homotopy perturbation method in solving nonlinear equations. J. Comput. Appl. Math. 228, 168–173 (2009) 34. Liang, S., Jeffrey, D.J.: Comparison of homotopy analysis method and homotopy perturbation method through an evolution equation. Commun. Nonlinear Sci. Numer. Simul. 14, 4057–4064 (2009) 35. Liang, S., Jeffrey, D.J.: An efficient analytical approach for solving fourth order boundary value problems. Comput. Phys. Commun. 180, 2034–2040 (2009) 36. Makukula, Z., Sibanda, P., Motsa, S.S.: A note on the solution of the Von Kármán equations using series and Chebyshev spectral methods. BVP 2010, Article ID 471793, 1–17 (2010). doi:10.1155/2010/471793 37. Makukula, Z., Sibanda, P., Motsa, S.S.: A novel numerical technique for two-dimensional Laminar flow between two moving porous walls. Math. Probl. Eng. 2010, Article ID 528956, 1–15 (2010). doi:10.1155/2010/528956 38. Makukula, Z.G., Sibanda, P., Motsa, S.S.: On new solutions for heat transfer in a visco-elastic fluid between parallel plates. M3AS 4(4), 221–230 (2010) 39. Makukula, Z., Motsa, S.S., Sibanda, P.: On a new solution for the viscoelastic squeezing flow between two parallel plates. JARAM 2, 31–38 (2010) 40. Chun-Mei, C., Gao, F.: A few numerical methods for solving nonlinear equations. Int. Math. Forum 3(29), 1437–1443 (2008) 41. Motsa, S.S., Sibanda, P., Shateyi, S.: A new spectral-homotopy analysis method for solving a nonlinear second order BVP. Commun. Nonlinear Sci. Numer. Simul. 15, 2293–2302 (2010) 42. Motsa, S.S., Sibanda, P., Awad, F.G., Shateyi, S.: A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem. Comp. Fluid. 39, 1219–1225 (2010) 43. Motsa, S.S., Sibanda, P.: On the solution of MHD flow over a nonlinear stretching sheet by an efficient semi-analytical technique. Int. J. Numer. Methods Fluids doi:10.1002/fld.2541 44. Motsa, S.S., Shateyi, S.: A new approach for the solution of three-dimensional magnetohydrodynamic rotating flow over a shrinking sheet. Math. Probl. Eng. 2010, Article ID 586340, 1–15 (2010). doi:10.1155/2010/586340 45. Motsa, S.S., Shateyi, S.: Successive linearisation solution of free convection non-Darcy flow with heat and mass transfer. In: El-Amin, M. (ed.) Advanced Topics in Mass Transfer, pp. 425– 438. InTech Open Access Publishers (2011) 46. Mohyud-Din, S.T., Noor, M.A.: Homotopy perturbation method and pad approximants for solving Flierl-Petviashivili equation. Appl. Appl. Math. 3(2), 224–234 (2008) 47. Osterle, J.F., J Young, F.: Natural convection between heated vertical plates in a horizontal magnetic field. J. Fluid Mech. 11, 512–518 (1961) 48. Sajida, M., Hayat, T.: Comparison of HAM and HPM methods in nonlinear heat conduction and convection equations. Nonlinear Anal.: Real World Appl. 9, 2296–2301 (2008) 49. Shampine, L.F., Gladwell, I., Thompson, S.: Solving ODEs with MATLAB. Cambridge University Press, Cambridge (2003) 50. Shateyi, S., Motsa, S.S.: Variable viscosity on magnetohydrodynamic fluid flow and heat transfer over an unsteady stretching surface with Hall effect. BVP 2010, Article ID 257568, 1–20 (2010). doi:10.1155/2010/257568 51. Trefethen, L.N.: Spectral Methods in MATLAB. SIAM (2000) 52. Umavathi, J.C.: A note on magnetoconvection in a vertical enclosure. Int. J. Non-Linear Mech. 31(3), 371–376 (1996) 53. Wu, T.M.: A new formula of solving nonlinear equations by Adomian and homotopy methods. Appl. Math. Comput. 172, 903–907 (2006)
Numer Algor (2012) 60:463–481
481
54. Wu, T.M.: A study of convergence on the Newton-homotopy continuation method. Appl. Math. Comput. 168, 1169–1174 (2005) 55. Yildirim, A.: Solution of BVPs for fourth-order integro-differential equations by using homotopy perturbation method. Comput. Math. Appl. 56, 3175–3180 (2008) 56. Ziabakhsh, Z., Domairry, G.: Solution of the laminar viscous flow in a semi-porous channel in the presence of a uniform magnetic field by using the homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 14, 1284–1294 (2009)
An improved spectral homotopy analysis method for MHD flow in a semi-porous channel Sandile Sydney Motsa · Stanford Shateyi · Gerald T. Marewo · Precious Sibanda
Received: 1 April 2011 / Accepted: 28 November 2011 / Published online: 20 December 2011 © Springer Science+Business Media, LLC 2011
Abstract In this paper we report on a novel method for solving systems of highly nonlinear differential equations by blending two recent semi-numerical techniques; the spectral homotopy analysis method and the successive linearisation method. The hybrid method converges rapidly and is an enhancement of the utility of the original spectral homotopy analysis method (Motsa et al., Commun Nonlinear Sci Numer Simul 15:2293–2302, 2010; Computer & Fluids 39:1219–1225, 2010) and an improvement on other recent semianalytical techniques. We illustrate the application of the method by solving a system of nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. A comparison with the numerical solution confirms the validity and accuracy of the technique and shows that the method converges rapidly and gives very accurate results.
S. S. Motsa · P. Sibanda School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, Republic of South Africa S. S. Motsa e-mail: [email protected] P. Sibanda e-mail: [email protected] S. Shateyi (B) Department of Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, Republic of South Africa e-mail: [email protected] G. T. Marewo Department of Mathematics, University of Swaziland, Private Bag 4, Kwaluseni, Swaziland e-mail: [email protected]
464
Numer Algor (2012) 60:463–481
Keywords Improved spectral homotopy analysis method · Coupled nonlinear equations · Laminar two-dimensional flow · Electrically conducting incompressible viscous fluid
1 Introduction Most engineering and physical systems are described by complex systems of nonlinear equations whose solution pose major challenges. The majority of these equations have no closed-form solutions, and although there already exist a plethora of numerical methods for finding approximate solutions to nonlinear equations, it has always been an interesting problem and a challenge to devise new and more efficient algorithms, particularly analytical or semianalytical techniques for finding solutions of nonlinear equations. Analytic techniques such as perturbation methods have been extensively used for many years, for example, the homotopy continuation method has been used since the 1930’s [53]. Traditionally, almost all perturbation methods depend on the ‘small parameter assumption’ that greatly constrained the application of such methods for the vast majority of nonlinear problems which do not involve such a small parameter [46]. Such methods are therefore valid for only nonlinear problems [34]. However new and more efficient techniques are constantly being developed and in the last three decades the delta-expansion method [26], the Adomian decomposition method [6, 7], the homotopy analysis method (HAM) [25, 28, 29, 31], the homotopy perturbation method (HPM) [23, 24], the variational iteration method (VIM) [21, 22] and the homotopy decomposition method [15] are some of the methods that have been developed and used extensively to solve problems of nonlinear mechanics. The Adomian technique provides a very simple but efficient numerical technique for solving large classes of nonlinear equations. The method gives the solution as an infinite series that should generally converge to an accurate solution. The main weakness of the method however is that it is often difficult to prove the convergence of the series solution [15]. The modified Adomian decomposition method developed by Abbasbandy [1] is fairly efficient in solving nonlinear problems although the formula is prone to diverge if a poor choice of the initial guess is made. A method combining the Adomian method and the homotopy continuation method was developed by Wu [53] to address the problem of divergence of Abbasbandy’s iteration. An alternative approach that combines the Newton-Raphson and the homotopy continuation method, the so-called Newton-homotopy continuation method was proposed by Wu [54]. In his paper some robust rules for the choice of the auxiliary homotopy function to avoid the problem of divergence of many traditional numerical methods were proposed. Other attempts at modification of the Adomian method can be found, for example, in [4, 5, 10, 11, 16]. The homotopy analysis method (HAM) was developed by Liao [25, 28] and has been used to solve nonlinear differential equations in different fields
Numer Algor (2012) 60:463–481
465
of science and engineering. Unlike the Adomian method and other nonperturbation techniques that cannot guarantee the convergence of the solution series, the HAM provides for a simple mechanism to control and adjust the convergence rate and the convergence region of the solution series by means of an embedded convergence controlling parameter [15, 31]. The homotopy analysis method provides certain rules for the construction of the solution such as the rule of solution expression and the rule of coefficient ergodicity. The former is necessary for the concept of homotopy while the rule of coefficient ergodicity is based on the completeness of base functions [30, 31]. Various attempts have been made to improve the performance of the HAM and the homotopy decomposition method is one such modification that sought to combine the Adomian and homotopy methods [15]. Further modifications of the HAM have been made by Bataineh et al. [12, 13] to obtain the socalled modified homotopy analysis method (MHAM). An optimal homotopyanalysis method (OHAM) has been used in a recent study by Liao [32]. The homotopy perturbation method (HPM) is a coupling of the traditional perturbation method and homotopy in topology and was proposed by JiHuan He [23, 24]. It has been extensively used to solve a wide class of functional equations including integro-differential equations, see for example [55]. The HPM is convenient, efficient and gives reasonable accuracy in solving nonlinear problems and integro-differential equations. However, as shown in [34, 35], the HPM, just like the VIM, is prone to give divergent approximations. The HPM was thus modified in [20] and the modified algorithm used in [8, 46] to test its validity and accuracy. Sajid and Hayat [48] explicitly proved that for the nonlinear problems where is different from −1, one cannot get convergent results by using HPM. They proved that the HPM results can be obtained as a special case of the HAM when = 1. Similar conclusions were independently made by Abbasbandy [2], as well as by Liang and Jeffrey [34]. It should be noted that HPM and VIM are not the same though they can sometimes produce the same results [34]. The idea behind the variational iteration method is the construction of a correction function by a general Lagrange multiplier which can be identified optimally via the variational theory Ganji et al. [19] while homotopy perturbation method uses an embedding parameter and was developed by combining the standard homotopy and perturbation methods. In light of the limitations of some current analytical and semi-analytical approaches, it is important that new, accurate and robust methods are developed. One possible approach is to combine analytical and numerical techniques to produce hybrid methods that do not suffer from the limitations of analytical approaches. One such approach is the spectral homotopy analysis method (SHAM) which combines the Chebyshev spectral collocation method and the homotopy analysis method. This method (see Motsa et al. [41]) seeks to address some of the limitations of the standard homotopy analysis method. This method has the advantage that it is more flexible than the standard HAM in that it allows for a wider choice of linear operators, admits a much wider
466
Numer Algor (2012) 60:463–481
range of admissible values and converges much faster than the standard homotopy analysis method, [36, 37, 41–43]. Abbasbandy and Shivanian [3] presented a procedure based on pseudo-spectral method to calculate multiple (dual) solutions of a model of mixed convection in a porous medium with boundary conditions on a semi-infinite interval which admit dual solutions. In this paper we present a new and improved spectral homotopy analysis method (ISHAM) for solving highly nonlinear and coupled systems of ordinary differential equations. We use this method to solve a system of highly nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. The improved spectral homotopy analysis method is a hybrid method that blends the spectral homotopy analysis method and the successive linearisation method (SLM). The SLM is (see [9, 38, 39, 44, 45, 50]) a recent method that has been successfully utilized in solving several nonlinear problems. In the present study we show that the ISHAM is an enhancement of the original SHAM with better accuracy and rapid convergence.
2 Governing equations In this section we solve a test problem, the coupled nonlinear equations that governs the laminar two-dimensional flow of an electrically conducting incompressible viscous fluid between two parallel plates subject to a transverse magnetic field. We demonstrate how the hybrid spectral-homotopy analysis method can be applied to find solutions of the governing coupled second- and fourth-order nonlinear equations. The problem considered is described more fully in [17, 47, 52, 56]. The governing dimensionless equations are; ∂u ∂v + = 0, ∂x ∂y
(2.1)
2 ∂u v ∂u ∂ 2u Ha2 2 2 ∂ Py 2∂ u +v = −ε + ε u, + − u ∂x ∂y ∂x hq ∂ x2 ∂ y2 Re 2 ∂ Py ∂v v ∂v ∂ 2u 2∂ u +v =− + ε u + 2 , ∂x ∂y ∂x hq ∂ x2 ∂y
(2.2) (2.3)
where h is the distance between the two plates, Ha is the Hartmann number, Re is the Reynolds number, q is the normal velocity at the porous wall and ε is a small aspect ratio. Using the Berman similarity transformation v = −V(y),
u = u0 U(y) + x
dV , dy
it can easily be shown, see [56], that equations (2.1)–(2.2) reduce to the coupled system U − Ha2 U + Re(VU − U V ) = 0, V
IV
− Ha V + Re(VV − V V ) = 0, 2
(2.4) (2.5)
Numer Algor (2012) 60:463–481
467
subject to the boundary conditions U = 1, V = 0, V = 0, y = 0,
(2.6)
U = 0, V = 1, V = 0, y = 1.
(2.7)
Equations (2.4)–(2.5) were solved recently by Ziabakhsh and Domairry [56] using the homotopy analysis method. In this paper the equations are used as a test problem to show the capability and prove the validity and accuracy of the ISHAM. The results are compared to the numerical solutions reported in [56] and the results of the standard SHAM and SLM.
3 Improved spectral homotopy analysis method In this section we describe the use of the improved spectral homotopy analysis method (ISHAM) in the governing equations (2.4)–(2.7). To apply the ISHAM we begin with the SLM idea, where it is assumed that the solutions U(y) and V(y) can be expanded as U(y) = U i (y) +
i−1
U n (y),
V(y) = Vi (y) +
n=0
i−1
Vn (y), i = 1, 2, 3, . . . (3.8)
n=0
where U i and Vi are unknown functions whose solutions are obtained using the SHAM approach [36, 37, 41–43] at the ith iteration and U n , Vn (1 ≤ n ≤ i − 1) are known from previous iterations. The algorithm starts with the initial approximations U 0 (y) and V0 (y) which are chosen to satisfy the boundary conditions (2.6)–(2.7). We choose the following initial guesses U 0 (y) = 1 − y,
V0 (y) = 3y2 − 2y3 .
(3.9)
which satisfy the boundary conditions (2.6)–(2.7). Substituting (3.8) in the governing equations (2.4)–(2.7) gives (1) U i + d1,i−1 U i + d2,i−1 U i + d3,i−1 Vi + d4,i−1 Vi + Re(Vi U i − U i Vi ) = si−1 ,
(3.10) (2) , Viiv + d1,i−1 Vi + d2,i−1 Vi + e1,i−1 Vi + e2,i−1 Vi + Re(Vi Vi − Vi Vi ) = si−1
(3.11) subject to boundary conditions U i (0) = Vi (0) = Vi (0) = U i (1) = Vi (1) = Vi (1) = 0,
(3.12)
468
Numer Algor (2012) 60:463–481
where d1,i−1 = Re
i−1
Vn ,
d2,i−1 = −Ha2 − Re
n=0
d3,i−1 = −Re
Un,
d4,i−1 = Re
n=0
(1) si−1
=−
U n
− Ha
2
i−1
U n + Re
=−
i−1
Vn ,
e2,i−1 = Re
Vniv
i−1 n=0
− Ha
n=0
2
i−1
Vn
n=0
n=0
(2) si−1
U n ,
n=0
i−1
i−1
i−1 n=0
n=0
e1,i−1 = −Re
Vn ,
n=0
i−1
i−1
i−1
Vn
+ Re
n=0
i−1
U n
−
n=0
i−1
Un
i−1
n=0
Vn
.
n=0
Vn i−1 n=0
Vn
i−1 n=0
Vn
−
i−1
Vn
n=0
i−1
Vn
n=0
The successive linearisation method (see [9, 38, 39, 44, 45, 50]) is based on iteratively solving the linear components of equations (3.10) and (3.11) which are given as (1) U i + d1,i−1 U i + d2,i−1 U i + d3,i−1 Vi + d4,i−1 Vi = si−1 ,
(3.13)
(2) Viiv + d1,i−1 Vi + d2,i−1 Vi + e1,i−1 Vi + e2,i−1 Vi = si−1 ,
(3.14)
subject to boundary conditions 3.12. Equations (3.13)–(3.14) can be inverted using any numerical method of solutions. In some cases, it may be possible to solve the equations using analytical means. In the case of the ISHAM, an attempt is made to solve the full set of equations (3.10)–(3.11) without linearizing. Starting from the initial approximations U 0 (y), V0 (y) the subsequent solutions U i , Vi (i ≥ 1) are obtained by recursively solving equations (3.10)–(3.12) using the SHAM [36, 37, 41–43]. To find the SHAM solutions of (3.10)–(3.12) we begin by defining the following linear operators L1 [U˜ i (y; q), V˜ i (y; q)] =
∂ 2 U˜ i ∂ y2
˜
˜
+ d1,i−1 ∂∂Uy i + d2,i−1 U˜ i + d3,i−1 ∂∂Vy i + d4,i−1 V˜ i , (3.15)
L2 [V˜ i (y; q)] =
∂ 4 V˜ i ∂ y4
3 ˜ 2 ˜ ˜ + d1,i−1 ∂∂ yV3 i + d2,i−1 ∂∂ yV2 i + e1,i−1 ∂∂Vy i + e2,i−1 V˜ i ,
(3.16) where q ∈ [0, 1] is an embedding parameter and U˜ i (y; q), V˜ i (y; q) are unknown functions. If we also define auxiliary nonlinear operators N1 [U˜ i (y; q), V˜ i (y; q)] = L1 [U˜ i (y; q), V˜ i (y; q)] + Re(V˜ i U˜ i − U˜ i V˜ i ),
(3.17)
Numer Algor (2012) 60:463–481
469
N2 [V˜ i (y; q)] = L2 [V˜ i (y; q)] + Re(V˜ i V˜ i − V˜ i V˜ i ),
(3.18)
then the zero order deformation equations for this problem are
(1 − q)L1 [U˜ i (y; q), V˜ i (y; q)] − [U i,0 (y), Vi,0 (y)]
(1) = q N1 [U˜ i (y; q), V˜ i (y; q)] − si−1 ,
(3.19)
(2) , (1 − q)L2 [V˜ i (y; q) − Vi,0 (y)] = q N2 [V˜ i (y; q)] − si−1
(3.20)
where U i,0 (y), Vi,0 (y) are initial approximations for the unknown functions U(y), V(y) respectively and = 0 is the converging controlling auxiliary parameter. The corresponding m-th order deformation equations are L1 [U i,m , Vi,m ] − χm [U i,m−1 , Vi,m−1 ] = Sm,1 ,
(3.21)
L2 [Vi,m − χm Vi,m−1 ] = Sm,2 ,
(3.22)
subject to boundary conditions (0) = U i,m (1) = Vi,m (1) = Vi,m (1) = 0. U i (0) = Vi,m (0) = Vi,m
where
χm =
0, 1,
⎛ Sm,1 = L1 [U i,m−1 , Vi,m−1 ] + Re ⎝
m−1
m≤1 m>1
−(1 − and
(3.24)
Vi, jU i,m−1− j−
j=0
m−1
⎞ ⎠ U i, j Vi,m−1− j
j=0
(1) , χm )si−1
⎛
Sm,2 = L2 [Vi,m−1 ] + Re ⎝
(3.23)
(3.25)
m−1 j=0
Vi, j Vi,m−1− j−
m−1
⎞ (2) ⎠ Vi, j Vi,m−1− j − (1 − χm )si−1 ,
j=0
(3.26) The unknown functions U i,0 (y), Vi,0 (y) are obtained as solutions to the boundary value problem (1) + d1,i−1 U i,0 + d2,i−1 U i,0 + d3,i−1 Vi,0 + d4,i−1 Vi,0 = si−1 , U i,0
(3.27)
(2) iv Vi,0 + d1,i−1 Vi,0 + d2,i−1 Vi,0 + e1,i−1 Vi,0 + e2,i−1 Vi,0 = si−1 ,
(3.28)
470
Numer Algor (2012) 60:463–481
subject to boundary conditions U i,0 (0) = Vi,0 (0) = Vi,0 (0) = U i,0 (1) = Vi,0 (1) = Vi,0 (1) = 0.
(3.29)
We use the Chebyshev pseudospectral method to solve equations (3.27)– (3.29). Before using the Chebyshev pseudospectral method to solve this problem it is convenient to first use the mapping η = 2y − 1; 0 ≤ y ≤ 1 to transform the problem domain from [0, 1] into the domain [−1, 1] on which the Chebyshev pseudospectral method is implemented. The unknown functions U i,0 (η) and Vi,0 (η) are approximated as truncated series of Chebyshev polynomials of the form U i,0 (η) =
N
U i,0 (ηk )Tk (η j) and Vi,0 (η) =
k=0
N
Vi,0 (ηk )Tk (η j),
(3.30)
k=0
where j = 0, 1, . . . , N, Tk is the kth Chebyshev polynomial, and η0 , η1 , . . . , η N are Gauss-Lobatto collocation points (see [14]) defined by η j = cos
πj , N
j = 0, 1, . . . , N.
(3.31)
Derivatives of the functions U i,0 (y) and Vi,0 (η) at the collocation points are represented as ds U i,0 s = 2 Dsjk U i,0 (ηk ), dys N
k=0
ds Vi,0 s = 2 Dsjk Vi,0 (ηk ), dys N
(3.32)
k=0
where s is the order of differentiation and D is the Chebyshev spectral differentiation matrix (see for example [14, 51]). Substituting equations (3.30)– (3.32) in (3.27)–(3.29) yields
B11 B12 O B22
(1) Ui,0 S = i−1 (2) Vi,0 Si−1
(3.33)
subject to boundary conditions
U i,0 (η N ) = Vi,0 (η N ) =
N
D Nk Vi,0 (η N ) = U i,0 (η0 )
k=0
= Vi,0 (η0 ) =
N k=0
D0k Vi,0 (η0 ) = 0,
(3.34)
Numer Algor (2012) 60:463–481
471
where B11 = 24 D2 + 2d1,i−1 D + d2,i−1 , B12 = 2d3,i−1 D + d4,i−1 , B22 = 24 D4 + 23 d1,i−1 D3 + 22 d2,i−1 D2 + 2e1,i−1 D + e2,i−1 , Ui,0 = [U i,0 (η0 ), U i,0 (η1 ), . . . , U i,0 (η N )]T , Vi,0 = [Vi,0 (η0 ), Vi,0 (η1 ), . . . , Vi,0 (η N )]T , (1) (1) (1) (1) Si−1 = [Si−1 (η0 ), Si−1 (η1 ), . . . , Si−1 (η N )]T , (2) (2) (2) (2) Si−1 = [Si−1 (η0 ), Si−1 (η1 ), . . . , Si−1 (η N )]T , (2) dr,i−1 = diag(dr,i−1 (η1 ), dr,i−1 (η2 ), . . . , dr,i−1 (η N )); r = 1, 2, 3, 4, (2) er,i−1 = diag(er,i−1 (η1 ), er,i−1 (η2 ), . . . , er,i−1 (η N )); r = 1, 2,
and O is an (N + 1) × (N + 1) zero matrix. Once boundary conditions (3.34) are imposed on (3.33), the linear system is solved for Ui,0 , Vi,0 to reveal U i,0 (y), Vi,0 (y). Applying the Chebyshev spectral method in a similar manner to equations (3.21)–(3.23) gives B11 B12 Ui,m B11 B12 Ui,m−1 = (χm + ) O B22 Vi,m O B22 Vi,m−1 (1) S P (3.35) + i,m−1 − (1 − χm ) i−1 (2) Qi,m−1 Si−1 subject to boundary conditions U i,m (η N ) = Vi,m (η N ) =
N
D Nk Vi,m (η N ) = U i,m (η0 )
k=0
= Vi,m (η0 ) =
N
D0k Vi,m (η0 ) = 0.
(3.36)
k=0
where Ui,m = [U i,m (η0 ), U i,m (η1 ), . . . , U i,m (η N )]T , Vi,m = [Vi,m (η0 ), Vi,m (η1 ), . . . , Vi,m (η N )]T , ⎛ ⎞ m−1 m−1 Pi,m−1 = Re ⎝ Vi, jDUi,m−1− j − Ui, jDVi,m−1− j⎠ , ⎛ Qi,m−1 = Re ⎝
j=0 m−1 j=0
j=0
Vi, jD3 Vi,m−1− j −
m−1 j=0
⎞ DVi, jD2 Vi,m−1− j⎠ .
472
Numer Algor (2012) 60:463–481
Once boundary conditions are imposed, the solution to equations (3.35)–(3.36) is obtained as ˜ i,m−1 + B−1 [Ci,m−1 − (1 − χm )Si,m−1 ]; m = 1, 2, . . . , Zi,m = (χm + )B−1 BZ (3.37) ˜ are coefficient matrices in (3.35) modified to incorporate where B and B boundary conditions, and
Zi,m
(1) Ui,m Pi,m−1 S = , Ci,m−1 = , Si−1 = i−1 (2) Vi,m Qi,m−1 Si−1
(3.38)
Thus, starting from the initial approximation, which is obtained from (3.33), higher order approximations U i,m (y), Vi,m (η), for m ≥ 1, can be obtained through the recursive formula (3.37). The solutions for U i (y), Vi (y) are then generated using the solutions for U i,m , Vi,m as follows U i = U i,0 + U i,1 + U i,2 + U i,3 + · · · + U i,m ,
(3.39)
Vi = Vi,0 + Vi,1 + Vi,2 + Vi,3 + · · · + Vi,m .
(3.40)
The [i, m] approximate solutions for U(y) and V(y) are then obtained by substituting U i and Vi (obtained from (3.39)–(3.40)) in (3.8). We remark that the recursive scheme (3.37) reduces to the original spectral homotopy analysis method [36, 37, 41–43] when i = 1 and the successive linearisation method [9, 38, 39, 44, 45, 50] when m = 0.
4 Results and discussion The governing equations (2.4)–(2.7) were solved using the ISHAM approach described in the previous section. In this section we present the results showing the effect of the governing parameters on the important properties of the flow. To check the accuracy of the proposed ISHAM, comparison was made with numerical solutions obtained using the MATLAB routine bvp4c [27, 49], which is an adaptive Lobatto quadrature scheme. In computing all the results presented in this work the number of collocation points used was N = 60. This value N was found to be appropriate for all the physical parameters used in this section and gave results that matched the results of the MATLAB routine bvp4c when the specified relative tolerance, RelTol, and the specified absolute tolerance, AbsTol were both less than e−9 . We note that the iterative scheme (3.37) that generates the approximate solutions U(y), V(y) depends on the convergence controlling auxiliary parameter . Hence, a careful selection of must be made. In the case of the original homotopy analysis method (HAM), it was pointed out by Liao [28],
Numer Algor (2012) 60:463–481
473 14.041 m=2 m=4 m = 10 m = 16
14.0405 m=2 m=4 m = 10 m = 16
14.04 14.0395 14.039 14.0385 14.038
0
14.0375
0
Fig. 1 −curve for U (0) and V (0) when Re = 10, and Ha = 0
that the accuracy and convergence of the HAM series solution depends on the careful selection of the auxiliary parameter . It was suggested in [28] that the valid values of where the series converges correspond to the horizontal segment of each -curve, which is a plot of some property of the independent variable (or it’s derivatives) at a certain point in it’s governing domain. In this work, the optimal curve was identified as the value of at the critical point of the second order −curves generated using U (0) and V (0). Using numerical experimentation with values taken along the horizontal segment of the curves, it was found that the critical values of the second order curve resulted in the most convergent series. An sample curve is shown in Fig. 1. Tables 1 and 2 give a comparison of the SHAM (order [1, m] of the ISHAM), the homotopy analysis method (HAM) results reported in [56], and numerical results obtained using the Matlab in-built boundary value problem solver bvp4c. In particular, Tables 1 and 2 give a snap shot of the mean velocity values at different locations within the channel. We note that these
Table 1 Comparison of the numerical results against the SHAM (order [1, m] on the ISHAM) approximate solutions for U(y) and V(y) when Ha = 0 and Re = 1 with N = 60 and = −1 y U(y) 0.00 0.25 0.50 0.75 1.00 V(y) 0.00 0.25 0.50 0.75 1.00
1st order
2nd order
3rd order
bvp4c
Ref.[56]
1.000000000 0.661815551 0.380389261 0.164599897 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815280 0.380388861 0.164599661 0.000000000
1.000000000 0.661815298 0.380388853 0.164599643 0.000000000
0.000000000 0.164565728 0.514940157 0.852203849 1.000000000
0.000000000 0.164565723 0.514940146 0.852203844 1.000000000
0.000000000 0.164565723 0.514940146 0.852203844 1.000000000
0.000000000 0.164565723 0.514940146 0.852203845 1.000000000
0.000000000 0.164565724 0.514940140 0.852203842 1.000000000
474
Numer Algor (2012) 60:463–481
Table 2 Comparison of the numerical results, the SHAM (order [1, m] on the ISHAM) solutions for U(y) and V(y) when Re = 1, Ha = 1 and = −1 y U(y) 0.00 0.25 0.50 0.75 1.00 V(y) 0.00 0.25 0.50 0.75 1.00
1st order
2nd order
3rd order
bvp4c
Ref. [56]
1.000000000 0.624376612 0.342902417 0.144030095 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376267 0.342901944 0.144029840 0.000000000
1.000000000 0.624376291 0.342901931 0.144029842 0.000000000
0.000000000 0.165179106 0.514478084 0.851140074 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179100 0.514478076 0.851140071 1.000000000
0.000000000 0.165179101 0.514478095 0.851140064 1.000000000
results rapidly converge to the bvp4c results, with full convergence achieved at order 2. This is a good indication of the improved accuracy of the new solution approach. It is also worth mentioning here that, in their paper, Ziabakhsh and Domairry [56] did not disclose the order of approximation used to generate the HAM results (although we estimate it cannot be anything less than the 10th order) nor do they reveal the value of used. Reproducing their results for purposes of comparing with the present approximations is therefore impractical. Nonetheless, their approximations do not yield the same level of accuracy as the current method. Tables 3 and 4 demonstrate the effect of the variation of Reynolds number on U (0) and V (0), respectively, when the Hartman number is fixed. The tables give a comparison of the SHAM ([1, m]), ISHAM ([i, m] i > 1) and Table 3 ISHAM results for U (0) when Ha = 1 using various values of Re Re
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
10
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.91 −0.91 −0.91 −0.91 −0.91 (SLM →) −0.68 −0.68 −0.68 −0.68 −0.68 (SLM →) −0.55 −0.55 −0.55 −0.55 −0.55
−3.67855093 −3.57357902 −3.57381695 −3.57381721 −3.57381721 −3.57381721 −9.09456068 −7.49394063 −7.48966570 −7.48908565 −7.48908410 −7.48908409 −14.45969153 −10.51500928 −10.46272106 −10.45596459 −10.45586600 −10.45586393
−3.57397714 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.53495863 −7.48908437 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.68732625 −10.45589567 −10.45586397 −10.45586388 −10.45586388 −10.45586388
−3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.48910428 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.45635538 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388
−3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −3.57381721 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −7.48908409 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388 −10.45586388
50
100
Numer Algor (2012) 60:463–481
475
Table 4 ISHAM results for V (0) when Ha = 1 using various values of Re Re
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
10
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.91 −0.91 −0.91 −0.91 −0.91 (SLM →) −0.68 −0.68 −0.68 −0.68 −0.68 (SLM →) −0.55 −0.55 −0.55 −0.55 −0.55
10.01969746 10.06518626 10.06486048 10.06485955 10.06485955 10.06485955 19.10301189 19.20789938 19.18030717 19.17886834 19.17885842 19.17885835 26.33729449 26.34830439 26.24922129 26.23538456 26.23501735 26.23500654
10.06478254 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.16893291 19.17885931 19.17885836 19.17885836 19.17885836 19.17885836 26.20696262 26.23503458 26.23500621 26.23500616 26.23500616 26.23500616
10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.17885554 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 26.23496043 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616
10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 10.06485955 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 19.17885836 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616 26.23500616
50
100
SLM ([i, 0]) generated results. The results presented in the tables were verified against the bvp4c numerical results. It can be seen that for large Reynolds numbers more terms in the SHAM solution series ar e required to achieve a match with the numerical results than in the case of the SLM. This shows that the SLM converges much faster than the SHAM. Essentially, the ISHAM results are a combination between and i steps of the SLM and m steps of the SHAM. The computing times required to generate the results are of the order of fractions of one second for all three methods. Hence, meaningful Table 5 ISHAM results for U (0) when Re = 10 using various values of Ha Ha 0
5
10
Order (m) 0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.90 −0.90 −0.90 −0.90 −0.90 (SLM →) −0.95 −0.95 −0.95 −0.95 −0.95 (SLM →) −0.97 −0.97 −0.97 −0.97 −0.97
[1,m](SHAM) −3.56623375 −3.46024971 −3.46050840 −3.46050933 −3.46050933 −3.46050933 −5.90262607 −5.82294691 −5.82300329 −5.82300332 −5.82300332 −5.82300332 −10.39354751 −10.33060703 −10.33062021 −10.33062025 −10.33062025 −10.33062025
[2,m] −3.46067933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82305801 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33063143 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
[3,m] −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
[4,m] −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −3.46050933 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −5.82300332 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025 −10.33062025
476
Numer Algor (2012) 60:463–481
Table 6 ISHAM results for V (0) when Re = 10 using various values of Ha Ha
Order (m)
[1,m](SHAM)
[2,m]
[3,m]
[4,m]
0
0 2 4 8 12 16 0 2 4 8 12 16 0 2 4 8 12 16
(SLM →) −0.90 −0.90 −0.90 −0.90 −0.90 (SLM →) −0.95 −0.95 −0.95 −0.95 −0.95 (SLM →) −0.97 −0.97 −0.97 −0.97 −0.97
9.98779654 10.03712703 10.03677402 10.03677181 10.03677181 10.03677181 10.94299138 10.94069274 10.94062918 10.94062917 10.94062917 10.94062917 14.04090221 14.03847281 14.03848032 14.03848033 14.03848033 14.03848033
10.03668989 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94059825 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03847155 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.03677182 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 10.94062917 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033 14.03848033
5
10
comparison between the three methods is done using the total number of iteration steps required to achieve full convergence to the numerical solutions. For the SHAM and SLM approaches the total number of iteration steps is equivalent to their respective orders m and i. Since for i iterations there are m SHAM steps of the ISHAM, it follows that the total number of iterations of the ISHAM is i × m. Tables 3 and 4 show that the SLM converges much faster than the SHAM for all values of Re. Comparing the SLM and ISHAM results, it can be noted that the SLM converges faster than the ISHAM. For example, when Re = 100 we note from Table 3 that at the ISHAM order [2, 4] (total of 8 iterations),
1
1
0.9
Ha = 0 Ha = 5 Ha = 20
0.8
0.9
Ha = 0 Ha = 5 Ha = 20
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2 0.2 0.1 0.1
0 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 2 The comparison between the numerical results and [2,2] order ISHAM solution for U(y), V(y) and V (y), for different values of Ha, when Re = 1, N = 60, = −1
Numer Algor (2012) 60:463–481
477 1
1 Re = 1 Re = 10 Re = 20
0.9
Re = 1 Re = 10 Re = 20
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 3 The comparison between the numerical results and [2,2] order ISHAM solution for U(y), V(y) and V (y), for different values of Re, when Ha = 1, N = 60, = −1
U (0) is −10.45586397 compared to the converged solution of −10.45586388 which is obtained after only 4 iterations of the SLM. The same observation can be made by considering Table 4 for V (0). By the same argument, it can be seen that the ISHAM converges more rapidly than then SHAM approach. Tables 5 and 6 demonstrate the effect of the variation of Hartman number on U (0) and V (0), respectively, when the Reynolds number is fixed. Increasing the Hartmann number has no noticeable effect on the convergence rate of the SHAM. Again, it can be seen from Tables 5 and 6 that the SLM converges much faster than ISHAM which in turn converges much faster than the SHAM. Figures 1, 2, 3 and 4 give a comparison between the [2,2] order ISHAM approximations and the numerical results for different Hartmann numbers and increasing Reynolds numbers. The parameter values used are the same as those used by Ziabakhsh and Domairry [56] in order to facilitate a comparison
1.6
1.6 Ha = 0 Ha = 5 Ha = 20
1.4
Re = 1 Re = 10 Re = 20
1.4
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
Fig. 4 Comparison between the numerical results and [2,2] order ISHAM solution for V (y) when Re = 1, Ha = 1 and N = 60, = −1
478
Numer Algor (2012) 60:463–481
of the results. We remark here that reasonable accuracy and matching of the numerical results and the ISHAM approximations is achieved at the lower orders of approximation of the SHAM, SLM and ISHAM. When Ha = 0 and Re = 1 the results presented here match exactly with the HAM results given by Ziabakhsh and Domairry [56]. As expected, increasing the Hartmann number has the effect of suppressing the flow of the fluid due to a resistive Lorentz force now slowing down the motion of the fluid.
5 Conclusion In this study, a new numerical perturbation scheme for solving complex nonlinear boundary value problems is presented. The method is based on a blend between the spectral homotopy analysis method (SHAM) and the successive linearisation method (SLM). The implementation of the improved spectral homotopy analysis method (ISHAM), was demonstrated by solving a system of nonlinear differential equations that govern the problem of laminar viscous flow in a semi-porous channel subject to a transverse magnetic field. Validation of the results obtained in this study was done against numerical results obtained using the bvp4c Matlab routine. The numerical experimentation conducted in this study suggests that the proposed method converges much more rapidly to the numerical when compared to the original SHAM. However, comparison of the convergence rates between the ISHAM and the SLM for this particular set of equations showed that the SLM convergence rate much more rapidly than the ISHAM. This preliminary study shows that the ISHAM is an enhancement and improvement of the accuracy and utility of the SHAM. Consequently, the proposed method would be expected to perform better, with more rapid convergence and better accuracy than the standard HAM and other recent semi-analytical methods when solving nonlinear boundary value problems that arise in science and engineering. Finally, it is worth noting that although in this study we have used the socalled -curves to find the value of the convergence-control parameter , it might be more worthwhile to explore the use of the so-called optimal HAM, Liao [32]. Acknowledgements The authors acknowledge the financial support from the National Research Foundation (NRF). We wish to thank the anonymous reviewers for their helpful comments and suggestions.
References 1. Abbasbandy, S.: Improving Newton-Raphson method for nonlinear equations by modified Adomian decomposition method. Appl. Math. Comput. 145, 887–893 (2003) 2. Abbasbandy, S.: The application of homotopy analysis method to solve a generalized HirotaSatsuma coupled KdV equation. Phys. Lett. A 361, 478–483 (2007)
Numer Algor (2012) 60:463–481
479
3. Abbasbandy, S., Shivanian, E.: Multiple solutions of mixed convection in a porous medium on semi-infinte interval using pseudo-spectral collocation method. Commun. Nonlinear Sci. Numer. Simul. 16, 2745–2752 (2011) 4. Abbaoui, K., Cherruault, Y.: Convergence of Adomian’s method applied to non-linear equations. Math. Comput. Model. 20, 69–73 (1994) 5. Abdulaziz, O., Hashim, I., Saif, A.: Series solutions of time-fractional PDEs by homotopy analysis method, differential equations and nonlinear mechanics, vol. 2008, Article ID 686512. doi:10.1155/2008/686512 (2008) 6. Adomian, G.: Nonlinear stochastic differential equations. J Math. Anal. Appl. 55, 441–452 (1976) 7. Adomian, G.: A review of the decomposition method and some recent results for nonlinear equations. Comput. Math. Appl. 21, 101–127 (1991) 8. Rafiq, A., Javeria, A.: New iterative methods for solving nonlinear equations by using the modified homotopy perturbation method. Acta Univ. Apulensis 18, 129–137 (2009) 9. Awad, F.G., Sibanda, P., Motsa, S.S., Makinde, O.D.: Convection from an inverted cone in a porous medium with cross-diffusion effects. Comput. Math. Appl. 61, 1431–1441 (2011) 10. Babolian, E., Biazar, J.: Solution of non-linear equations by modified Adomian decomposition method. Appl. Math. Comput. 132, 167–172 (2002) 11. Babolian, E., Biazar, J., Vahidi, A.R.: On the decomposition method for system of linear equations and system of linear Volterra integral equations. Appl. Math. Comput. 147, 19–27 (2004) 12. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: Modified homotopy analysis method for solving systems of second-order BVPs. Commun. Nonlinear Sci. Numer. Simul. 14, 430–442 (2009) 13. Bataineh, A.S., Noorani, M.S.M., Hashim, I.: On a new reliable modification of homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 14, 409–423 (2009) 14. Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods in Fluid Dynamics. Springer, Berlin (1988) 15. Chen, X., L, Zheng, Zhang, X.: Convergence of the homotopy decomposition method for solving nonlinear equations. Adv. Dyn. Syst. Appl. 2(1), 59–64 (2007) 16. Cherruault, Y.: Convergence of Adomian’s method. Kybernetes 18, 31–38 (1989) 17. Desseaux, A.: Influence of a magnetic field over a laminar viscous flow in a semi-porous channel. Int. J. Eng. Sci. 37, 1781–1794 (1999) 18. Don, W.S., Solomonoff, A.: Accuracy and speed in computing the Chebyshev Collocation Derivative. SIAM J. Sci. Comput. 16(6), 1253–1268 (1995) 19. Ganji, D.D., Nourollahi, M., Rostamian, M.: A comparison of variational iteration method with Adomian’s decompostion method in some highly nonlinear equations. Int. J. Sci. Technol. 2(2), 179–188 (2007) 20. Golbabaia, A., Javidi, M.: Application of homotopy perturbation method for solving eighthorder boundary value problems. Appl. Math. Comput. 191, 334–346 (2007) 21. He, J.H.: Variational iteration method: a kind of non-linear analytical technique: some examples. Int. J. Non-Linear Mech. 34(4), 699–708 (1999) 22. He, J.H.: Variational iteration method for autonomous ordinary differential systems. Appl. Math. Comput. 114(2–3), 115–123 (2000) 23. He, J.H.: Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 178, 257– 262 (1999) 24. He, J.H.: A coupling method of a homotopy technique and a perturbation technique for nonlinear problems. Int. J. Non-linear Mech. 35, 37–43 (2000) 25. Liao, S.J.: The proposed homotopy analysis technique for the solution of nonlinear problems. PhD thesis, Shanghai Jiao Tong University (1992) 26. Karmishin, A.M., Zhukov, A.I., Kolsov, V.G.: Methods of Dynamics Calculation and Testing for Thin-walled Structures. Mashinostroyenie, Moscow (1990) 27. Kierzenka, J., Shampine, L.: A BVP solver based on residual control and the Matlab PSE. ACM Trans. Math. Softw. 27, 299–316 (2001) 28. Liao, S.J.: Beyond Perturbation: Introduction to Homotopy Analysis Method. Chapman & Hall/CRC Press (2003) 29. Liao, S.J.: Comparison between the homotopy analysis method and the homotopy perturbation method. Appl. Math. Comput. 169, 1186–1194 (2005)
480
Numer Algor (2012) 60:463–481
30. Liao, S.J., Tan, Y.: A general approach to obtain series solutions of nonlinear differential equations. Stud. Appl. Math. 119, 297–354 (2007) 31. Liao, S.J.: Notes on the homotopy analysis method: Some definitions and theores. Commun. Nonlinear Sci. Numer. Simul. 14, 983–997 (2009) 32. Liao, S.J.: An optimal homotopy-analysis approach for strongly nonlinear differential equations. Commun. Nonlinear Sci. Numer. Simul. 15, 2003–2016 (2010) 33. Li, J.L.: Adomian’s decomposition method and homotopy perturbation method in solving nonlinear equations. J. Comput. Appl. Math. 228, 168–173 (2009) 34. Liang, S., Jeffrey, D.J.: Comparison of homotopy analysis method and homotopy perturbation method through an evolution equation. Commun. Nonlinear Sci. Numer. Simul. 14, 4057–4064 (2009) 35. Liang, S., Jeffrey, D.J.: An efficient analytical approach for solving fourth order boundary value problems. Comput. Phys. Commun. 180, 2034–2040 (2009) 36. Makukula, Z., Sibanda, P., Motsa, S.S.: A note on the solution of the Von Kármán equations using series and Chebyshev spectral methods. BVP 2010, Article ID 471793, 1–17 (2010). doi:10.1155/2010/471793 37. Makukula, Z., Sibanda, P., Motsa, S.S.: A novel numerical technique for two-dimensional Laminar flow between two moving porous walls. Math. Probl. Eng. 2010, Article ID 528956, 1–15 (2010). doi:10.1155/2010/528956 38. Makukula, Z.G., Sibanda, P., Motsa, S.S.: On new solutions for heat transfer in a visco-elastic fluid between parallel plates. M3AS 4(4), 221–230 (2010) 39. Makukula, Z., Motsa, S.S., Sibanda, P.: On a new solution for the viscoelastic squeezing flow between two parallel plates. JARAM 2, 31–38 (2010) 40. Chun-Mei, C., Gao, F.: A few numerical methods for solving nonlinear equations. Int. Math. Forum 3(29), 1437–1443 (2008) 41. Motsa, S.S., Sibanda, P., Shateyi, S.: A new spectral-homotopy analysis method for solving a nonlinear second order BVP. Commun. Nonlinear Sci. Numer. Simul. 15, 2293–2302 (2010) 42. Motsa, S.S., Sibanda, P., Awad, F.G., Shateyi, S.: A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem. Comp. Fluid. 39, 1219–1225 (2010) 43. Motsa, S.S., Sibanda, P.: On the solution of MHD flow over a nonlinear stretching sheet by an efficient semi-analytical technique. Int. J. Numer. Methods Fluids doi:10.1002/fld.2541 44. Motsa, S.S., Shateyi, S.: A new approach for the solution of three-dimensional magnetohydrodynamic rotating flow over a shrinking sheet. Math. Probl. Eng. 2010, Article ID 586340, 1–15 (2010). doi:10.1155/2010/586340 45. Motsa, S.S., Shateyi, S.: Successive linearisation solution of free convection non-Darcy flow with heat and mass transfer. In: El-Amin, M. (ed.) Advanced Topics in Mass Transfer, pp. 425– 438. InTech Open Access Publishers (2011) 46. Mohyud-Din, S.T., Noor, M.A.: Homotopy perturbation method and pad approximants for solving Flierl-Petviashivili equation. Appl. Appl. Math. 3(2), 224–234 (2008) 47. Osterle, J.F., J Young, F.: Natural convection between heated vertical plates in a horizontal magnetic field. J. Fluid Mech. 11, 512–518 (1961) 48. Sajida, M., Hayat, T.: Comparison of HAM and HPM methods in nonlinear heat conduction and convection equations. Nonlinear Anal.: Real World Appl. 9, 2296–2301 (2008) 49. Shampine, L.F., Gladwell, I., Thompson, S.: Solving ODEs with MATLAB. Cambridge University Press, Cambridge (2003) 50. Shateyi, S., Motsa, S.S.: Variable viscosity on magnetohydrodynamic fluid flow and heat transfer over an unsteady stretching surface with Hall effect. BVP 2010, Article ID 257568, 1–20 (2010). doi:10.1155/2010/257568 51. Trefethen, L.N.: Spectral Methods in MATLAB. SIAM (2000) 52. Umavathi, J.C.: A note on magnetoconvection in a vertical enclosure. Int. J. Non-Linear Mech. 31(3), 371–376 (1996) 53. Wu, T.M.: A new formula of solving nonlinear equations by Adomian and homotopy methods. Appl. Math. Comput. 172, 903–907 (2006)
Numer Algor (2012) 60:463–481
481
54. Wu, T.M.: A study of convergence on the Newton-homotopy continuation method. Appl. Math. Comput. 168, 1169–1174 (2005) 55. Yildirim, A.: Solution of BVPs for fourth-order integro-differential equations by using homotopy perturbation method. Comput. Math. Appl. 56, 3175–3180 (2008) 56. Ziabakhsh, Z., Domairry, G.: Solution of the laminar viscous flow in a semi-porous channel in the presence of a uniform magnetic field by using the homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 14, 1284–1294 (2009)
