# On New Numerical Techniques for the MHD Flow Past a Shrinking Sheet with Heat and Mass Transfer in the Presence of a Chemical Reaction

Authors: Z. G. Makukula; P. Sibanda; S. S. Motsa; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We use recent innovative solution techniques to investigate the problem of MHD viscous flow due to a shrinking sheet with a chemical reaction. A comparison is made of the convergence rates, ease of use, and expensiveness (the number of iterations required to give convergent results) of three seminumerical techniques in solving systems of nonlinear boundary value problems. The results were validated using a multistep, multimethod approach comprising the use of the shooting method, the Matlab bvp4c numerical routine, and with results in the literature.


---

## Body

## 1. Introduction

Boundary layer flow over a stretching surface occurs in several engineering processes such as hot rolling, wire drawing, and glass-fibre production. Materials that are manufactured by extrusion processes and heat-treated substances proceeding between a feed roll and a wind-up roll can be classified as a continuously stretching surface [1–3]. A shrinking film is useful in the packaging of bulk products since it can be unwrapped easily with adequate heat [4–7]. Shrinking problems can also be applied to the study of capillary effects in small pores and the hydraulic properties of agricultural clay soils [8]. Studies of flow due to a shrinking sheet with heat transfer and/or mass transfer have been considered by, among others, [7, 9].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 (ADM) [10, 11], differential transform method (DTM) [12], variational iteration method (VIM) [13], homotopy analysis method (HAM) [14–17], and Homotopy perturbation method (HPM) [18–23].Motsa and Shateyi [24] obtained a numerical solution of magnetohydrodynamic (MHD) and rotating flow over a porous shrinking sheet by the new approach known as spectral homotopy analysis method (SHAM). Muhaimin et al. [5] studied magnetohydrodynamic viscous flow due to a shrinking sheet in the presence of suction. The study found out that the shrinking of the sheet has a substantial effect on the flow field and, thus, on the heat and mass transfer rate from the sheet to the fluid.In this paper we provide a qualitative assessment of key features of three recent seminumerical techniques, namely, the successive linearisation method (SLM), the spectral-homotopy analysis method (SHAM), and the improved spectral-homotopy analysis method (ISHAM). The two methods were introduced and used by Motsa and his coworkers (see Motsa et al. [25, 26] and Makukula et al. [27–30]) to solve nonlinear boundary value problems. In Motsa et al. [25, 26, 29] the SHAM approach was tested on simple one-dimensional nonlinear boundary value problems. Later, Makukula et al. [28, 30, 31] extended the application of the SHAM to a system of two coupled nonlinear equations that model the von Kármán fluid flow problem. The SLM method was applied on one-dimensional nonlinear differential equations in Makukula et al. [27]. In this study we solve the nonlinear equations that govern the shrinking sheet problem for purposes of evaluating the efficiency of each method with regards to speed of convergence, ease of use, and expensiveness (in terms of the number of iterations required to give convergent results). We introduce the ISHAM as a method that is meant to improve the accuracy of the standard SHAM approach. The governing equations for the problem are a rather formidable system of three nonlinear differential equations in three unknowns. Parametric study of the effect of different parameters is made and the results compared with previous findings in the literature (see Noor et al.  [6], Mohd  and Hashim  [7], and Muhaimin et al. [5]). The solutions are further compared with results obtained using the shooting method and the bvp4c solver, which is based on Runge-Kutta fourth-order schemes.


## 2. Mathematical Formulation

We investigate the effect of chemical reaction, heat and mass transfer on nonlinear MHD boundary layer past a porous shrinking sheet with suction.  The governing boundary layer equations of momentum, energy, and mass diffusion in terms of the velocity componentsu, v, and w are (see Muhaimin et al. [5])(2.1)∂u∂x+∂v∂y+∂w∂z=0,u∂u∂x+v∂u∂y+w∂u∂z=1ρ∂p∂x+ν{∂2u∂x2+∂2u∂y2+∂2u∂z2}-uσB02ρ-uνK,u∂v∂x+v∂v∂y+w∂v∂z=1ρ∂p∂y+ν{∂2v∂x2+∂2v∂y+∂2v∂z2}-vσB02ρ-vνK,u∂w∂x+v∂w∂y+w∂w∂z=1ρ∂p∂z+ν{∂2w∂x2+∂2w∂y2+∂2w∂z2},u∂T∂x+v∂T∂y+w∂T∂z=α{∂2T∂x2+∂2T∂y2+∂2T∂z2},u∂C∂x+v∂C∂y+w∂C∂z=D{∂2C∂x2+∂2C∂y2+∂2C∂z2}-k1C,whereα is the thermal conductivity of the fluid, B0 is the magnetic field, κ is the thermal viscosity, K is the permeability of the porous medium, k1 is the rate of chemical reaction, ν=μ/ρ is the kinetic viscosity,  μ is the dynamic viscosity, and σ is the electrical conductivity.The applicable boundary conditions are(2.2)u=-ax,v=-a(m-1)y,w=-W,T=Tw,C=Cw,aty=0,u⟶0,C⟶C∞,T⟶T∞asy⟶∞,
where a>0 is the shrinking constant and W is the suction velocity. The cases m=1 and m=2 correspond to shrinking sheets in the x- and y-directions, respectively.Using the similarity transformations (see Sajid and Hayat [32]):
(2.3)u=axF′(η),v=a(m-1)yF′(η),w=-aνmF(η),η=aνz,θ(η)=(T-T∞)(Tw-T∞),ϕ(η)=(C-C∞)(Cw-C∞),(2.1) are transformed to the system of nonlinear equations(2.4)f′′′-(M2+λPr)f′-f′2+mff′′=0,(2.5)θ′′-Prf′θ+mPrfθ′=0,(2.6)ϕ′′-Scf′ϕ+mScfϕ′-Scγϕ=0,
subject to(2.7)f(0)=s,f′(0)=-1,f′(∞)=0,θ(0)=1,θ(∞)=0,ϕ(0)=1,ϕ(∞)=0,
where Pr=ν/κ is the Prandtl number, Sc=ν/D is the Schmidt number, λ=κ/Ka is the porosity, and γ is the chemical reaction parameter. We remark that (2.4) can be solved independently of equations of (2.5)-(2.6) for f,  but the solutions for θ and ϕ directly depend on the solution for f. To demonstrate how robust the proposed methods of solution are, the system of (2.4)-(2.5) is solved simultaneously in the next section. Solving the equations simultaneously is also important when conducting the parametric study because some of the governing parameters such as Pr and m affect all three unknown variables.


## 3. Solution Methods

We solve (2.4)–(2.6) using three recent innovative semi-numerical methods. Validation of the results is done by further solving the equations numerically using a shooting method and the Matlab bvp4c solver. For the last two methods we used a tolerance of 10-6.We begin by transforming the domain[0,∞) to [-1,1], using the domain truncation method, the domain [0,∞) is first approximated by the computational domain [0,L], where L is a fixed length that is taken to be larger than the thickness of the boundary layer. The domain [0,L] is then transformed to [-1,1] using the algebraic mapping(3.1)ξ=2ηL-1,ξ∈[-1,1].


### 3.1. The Successive Linearisation Method (SLM)

The successive linearisation method (see Makukula et al. [27, 28]) is used to solve (2.4)–(2.7). The starting point is to assume that the independent variables f(η), θ(η), and ϕ(η) may be expanded as(3.2)f(η)=fi(η)+∑m=0i-1Fm(η),θ(η)=θi(η)+∑m=0i-1Θm(η),ϕ(η)=ϕi(η)+∑m=0i-1Φm(η),i=1,2,3,…,
where fi, θi and ϕi, are unknown functions and Fm, Θm,  and Φm(m≥1) are approximations that are obtained by recursively solving the linear part of the equation that results from substituting (3.2) in the governing equations (2.4)–(2.7). Substituting (3.2) in the governing equations (2.4)–(2.7) gives    (3.3)fi′′′+a1,i-1fi′′+a2,i-1fi′+a3,i-1fi+mff′′-f′2=r1,i-1,θi′′+b1,i-1θi′+b2,i-1θi+b3,i-1fi′+b4,i-1fi-Prf′θ+mPrfθ′=r2,i-1,ϕi′′+c1,i-1ϕi′+c2,i-1ϕi+c3,i-1fi′+c4,i-1fi-Scf′ϕ+mScfϕ′=r3,i-1,
where the coefficient parameters ak,i-1,  bk,i-1, ck,i-1 (k=1,…,4), r1,i-1, r2,i-1,  and r3,i-1 are defined as(3.4)a1,i-1=m∑m=0i-1Fm,a2,i-1=-(2∑m=0i-1Fm′+M2+λPr),a3,i-1=m∑m=0i-1Fm′′,b1,i-1=mPr∑m=0i-1Fm,b2,i-1=-Pr∑m=0i-1Fm′,b3,i-1=-Pr∑m=0i-1Θm,b4,i-1=mPr∑m=0i-1Θm′,c1,i-1=mSc∑m=0i-1Fm,c2,i-1=-Sc∑m=0i-1Fm′-Scγ,c3,i-1=-Sc∑m=0i-1Φm,c4,i-1=mSc∑m=0i-1Φm′,r1,i-1=-[∑m=0i-1Fm′′′-∑m=0i-1Fm′∑m=0i-1Fm′+m∑m=0i-1Fm′′∑m=0i-1Fm-(M2+λPr)∑m=0i-1Fm′],r2,i-1=-[∑m=0i-1Θm′′-Pr∑m=0i-1Fm′∑m=0i-1Θm+mPr∑m=0i-1Fm∑m=0i-1Θm′],r3,i-1=-[∑m=0i-1Φm′′-Sc∑m=0i-1Fm′∑m=0i-1Φm+mSc∑m=0i-1Fm∑m=0i-1Φm′-Scγ∑m=0i-1Φ].Starting from the initial approximations(3.5)F0(η)=s+e-2η-e-η,Θ0(η)=e-η,Φ0(η)=e-η,
which are chosen to satisfy the boundary conditions (2.7), the subsequent solutions for Fm, Θm,  and Φm, m≥1,  are obtained by successively solving the linearized form of (3.3) which are(3.6)Fi′′′+a1,i-1Fi′′+a2,i-1Fi′+a3,i-1Fi=r1,i-1,Θi′′+b1,i-1Θi′+b2,i-1Θi+b3,i-1Fi′+b4,i-1Fi=r2,i-1,Φi′′+c1,i-1Φi′+c2,i-1Φi+c3,i-1Fi′+c4,i-1Fi=r3,i-1,
subject to the boundary conditions(3.7)Fi(0)=Fi′(0)=Fi′(∞)=0,Θi(0)=0,Φi(0)=0.
Once each solution for Fi, Θi,  and Φi (i≥1) has been found from iteratively solving (3.6)-(3.7) for each i, the approximate solutions for f(η), θ(η),  and ϕ(η) are obtained as(3.8)f(η)≈∑m=0MFm(η),θ(η)≈∑m=0MΘm(η),ϕ(η)≈∑m=0MΦm(η).
In coming up with (3.8), we have assumed that (3.9)lim⁡i→∞fi=0,lim⁡i→∞θi=0,lim⁡i→∞ϕi=0.Equations (3.6)-(3.7) are integrated using the Chebyshev spectral collocation method (Canuto et al. [33] and Trefethen [34]). The unknown functions are defined by the Chebyshev interpolating polynomials with the Gauss-Lobatto points defined as(3.10)yj=cos⁡πjN,j=0,1,…,N,
where N is the number of collocation points used. The unknown functions Fi, Θi,  and Φi are approximated at the collocation points by(3.11)Fi(ξ)≈∑k=0NFi(ξk)Tk(ξj),Θi(ξ)≈∑k=0NΘi(ξk)Tk(ξj),Φi(ξ)≈∑k=0NΦi(ξk)Tk(ξj),j=0,1,…,N,
where Tk is the kth Chebyshev polynomial defined as(3.12)Tk(ξ)=cos⁡[kcos⁡-1(ξ)].
The derivatives at the collocation points are represented as(3.13)daFidηa=∑k=0NDkjaFi(ξk),daΘidηa=∑k=0NDkjaΘi(ξk),daΦidηa=∑k=0NDkjaΦi(ξk),j=0,1,…,N,
where a is the order of differentiation and D=(2/L)𝒟 with 𝒟 being the Chebyshev spectral differentiation matrix. Substituting (3.13) in (3.6)-(3.7) leads to the matrix equation(3.14)Ai-1Xi=Ri-1,
where Bi-1 is a 3(N+1)×3(N+1) square matrix and Xi and Pi-1 are 3(N+1)×1 column vectors defined by(3.15)Ai-1=(D3+a1D2+a2D+a3I0I0Ib3D+b4ID2+b1D+b2I0Ic3D+c4I0ID2+c1D+c2I),Xi=[Fm(ξ0),Fm(ξ1),…,Fm(ξN),Θm(ξ0),Θm(ξ1),…,Θm(ξN),Φm(ξ0),Φm(ξ1),…,Φm(ξN)]T,Ri-1=[r1,i-1(η0),r1,i-1(η1),…,r1,i-1(ηN),r2,i-1(η0),r2,i-1(η1),…,r2,i-1(ηN),r3,i-1(η0),r3,i-1(η1),…,r3,i-1(ηN)]T.In the above definitions,ak,i-1, and bk,i-1, ck,i-1 (k=1,…,4) are diagonal matrices of size (N+1)×(N+1). After modifying the matrix system (3.14) to incorporate the boundary conditions, the solution is obtained as(3.16)Xi=Ai-1-1Ri-1.


### 3.2. Spectral-Homotopy Analysis Method (SHAM)

The spectral-homotopy analysis method (SHAM) has been used by Motsa et al. [25, 26]. It is also convenient to first ensure that the boundary conditions are made homogeneous by using the transformations(3.17)f(η)=F(ξ)+f0(η),θ(η)=Θ(ξ)+θ0(η),ϕ(η)=Φ(ξ)+ϕ0(η),
where f0(η), and  θ0(η), ϕ0(η) are chosen to satisfy the boundary conditions (2.7) of the governing equations (2.4)–(2.6). From (3.1) and the chain rule, we have that(3.18)f′(η)=2LF′(ξ)+f0′(η),f′′(η)=4L2F′′(ξ)+f0′′(η),f′′′(η)=8L3F′′′(ξ)+f0′′′(η),θ′(η)=2LΘ′(ξ)+θ0′(η),θ′′(η)=4L2Θ′′(ξ)+θ0′′(η),ϕ′(η)=2LΦ′(ξ)+ϕ0′(η),ϕ′′(η)=4L2Φ′′(ξ)+ϕ0′′(η).
Substituting (3.1) and (3.17)-(3.18) in the governing equations and boundary conditions gives(3.19)a0F′′′+a1F′′+a2F′+a3F+4L2mF′′F-4L2F′F′=r1(η),b0Θ′′+b1Θ′+b2Θ+b3F′+b4F-2LPrF′Θ+2LmPrFΘ′=r2(η),c0Φ′′+c1Φ′+c2Φ+c3F′+c4F-2LScF′Φ+2LmScFΦ′=r3(η),
where prime now denotes derivative with respect to ξ and (3.20)a0=8L3,a1=4L2mf0,a2=-(2L(M2+λPr)+4Lf0′),a3=mf0′′,r1(η)=-(f0′′′+mf0f0′′-f0′f0′-(M2+λPr)f0′),b0=4L2,b1=2LmPrf0,b2=-Prf0′,b3=-2LPrθ0,b4=mPrθ0′,r2(η)=-(θ′′-Prf0′θ0+mPrf0θ0),c0=4L2,c1=2LmScf0,c2=-Sc(f0′+γ),c3=-2LScϕ0,c4=mScϕ0′,r3(η)=-(ϕ′′-Scf0′ϕ0+mScf0ϕ0′-Scγϕ).The initial guesses used are(3.21)f0(η)=s+e-2η-ηe-η,θ0(η)=e-η,ϕ0(η)=e-η.
Solving the linear part of the equation system (3.19), that is,(3.22)a0F0′′′+a1F0′′+a2F0′+a3F0=r1(η),b0Θ0′′+b1Θ0′+b2Θ0+b3F0′+b4F0=r2(η),c0Φ0′′+c1Φ0′+c2Φ0+c3F0′+c4F0=r3(η),
subject to(3.23)F0(-1)=2LF0′(-1)=2LF0′(1)=0,Θ0(-1)=Θ0(1)=0,Φ0(-1)=Φ0(1)=0,
will yield the initial SHAM approximate solution. Applying the Chebyshev pseudospectral method on equations (3.22)-(3.23) yields the matrix form(3.24)BY0=R,
where(3.25)B=(a0D3+a1D2+a2D+a3I0I0Ib3D+b4Ib0D2+b1D+b2I0Ic3D+c4I0Ic0D2+c1D+c2I),(3.26)R=[r1(η0),r1(η1),…,r1(ηN),r2(η0),r2(η1),…,r2(ηN),r3(η0),r3(η1),…,r3(ηN)]T,Y0=[F0(ξ0),F0(ξ1),…,F0(ξN),Θ0(ξ0),Θ0(ξ1),…,Θ0(ξN),Φ0(ξ0),Φ0(ξ1),…,Φ0(ξN)]T,ai=diag⁡([ai(η0),…,ai(ηN-1),ai(ηN)]),bi=diag⁡([bi(η0),…,bi(ηN-1),bi(ηN)]),ci=diag⁡([ci(η0),…,ci(ηN-1),ci(ηN)]),i=0,1,2,3,4.The superscript T denotes the transpose, diag is a diagonal matrix, and I is an identity matrix of size (N+1)×(N+1). The boundary conditions (3.23) are implemented in matrix B and vector R of equation (3.24). The values of [Y0(ξ1),Y0(ξ2),…,Y0(ξN-1)] are then determined from the  following equation:(3.27)Y0=B-1R,
which provides us with the initial approximation for the solution of the governing equations (3.19). With the initial approximate solution, we then find approximate solutions for the nonlinear equations (3.19). We start by defining the following linear operators:(3.28)LF[F̃(ξ;q)]=a0∂3F̃∂ξ3+a1∂2F̃∂ξ2+a2∂F̃∂ξ+a3F̃,LΘ[F̃(ξ;q),Θ̃(ξ;q)]=b0∂2Θ̃∂ξ2+b1∂Θ̃∂ξ+b2Θ̃+b3∂F̃∂ξ+b4F̃,LΦ[F̃(ξ;q),Φ̃(ξ;q)]=c0∂2Φ̃∂ξ2+c1∂Φ̃∂ξ+c2Φ̃+c3∂F̃∂ξ+c4F̃,
where q∈[0,1] is the embedding parameter and F̃(ξ;q), Θ̃(ξ;q), and Φ̃(ξ;q) are unknown functions. The zeroth-order deformation equations are given by(3.29)(1-q)LF[F̃(ξ;q)-F0(ξ)]=qℏ{NF[F̃(ξ;q)]-r1},(1-q)LΘ[Θ̃(ξ;q)-Θ0(ξ)]=qℏ{NΘ[F̃(ξ;q),Θ̃(ξ;q)]-r2},(1-q)LΦ[Φ̃(ξ;q)-Φ0(ξ)]=qℏ{NΦ[F̃(ξ;q),Φ̃(ξ;q)]-r3},
where ℏ is the nonzero convergence controlling auxiliary parameter and 𝒩F, 𝒩Θ,  and 𝒩Φ are nonlinear operators given by(3.30)NF[F̃(ξ;q)]=a0∂3F̃∂ξ3+a1∂2F̃∂ξ2+a2∂F̃∂ξ+a3F̃+4L2mF̃∂2F̃∂ξ2-4L2∂F̃∂ξ∂F̃∂ξ,NΘ[F̃(ξ;q),Θ̃(ξ;q)]=b0∂2Θ̃∂ξ2+b1∂Θ̃∂ξ+b2Θ̃+b3∂F̃∂ξ+b4F̃+2LPr(-Θ̃∂F̃∂ξ+mF̃∂Θ̃∂ξ),NΦ[F̃(ξ;q),Φ̃(ξ;q)]=c0∂2Φ̃∂ξ2+c1∂Φ̃∂ξ+c2Φ̃+c3∂F̃∂ξ+c4F̃+2LSc(-Φ̃∂F̃∂ξ+mF̃∂Φ̃∂ξ).
The m-th order deformation equations are given by(3.31)LF[Fm(ξ)-χmFm-1(ξ)]=ℏRmF,LΘ[Θm(ξ)-χmΘm-1(ξ)]=ℏRmΘ,LΦ[Φm(ξ)-χmΦm-1(ξ)]=ℏRmΦ,
subject to the boundary conditions(3.32)Fm(-1)=Fm′(-1)=Fm′(1)=0,Θm(-1)=Θm(1)=0,Φm(-1)=Φm(1)=0,
where(3.33)RmF(ξ)=a0Fm-1′′′+a1Fm-1′′+a2Fm-1′+a3Fm-1+4L2∑n=0m-1(-Fn′Fm-1-n′+mFnFm-1-n′′)-r1(η)(1-χm),RmΘ(ξ)=b0Θm-1′′+b1Θm-1′+b2Θm-1+b3Fm-1′+b4Fm-1+2LPr∑n=0m-1(-Fn′Θm-1-n+mΘn′Fm-1-n)-r2(η)(1-χm),RmΦ(ξ)=c0Φm-1′′+c1Φm-1′+c2Φm-1+c3Fm-1′+c4Fm-1+2LSc∑n=0m-1(-Fn′Φm-1-n+mΦn′Fm-1-n)-r3(η)(1-χm),(3.34)χm={0,m≤11,m>1.
Applying the Chebyshev pseudospectral transformation to equations (3.31)–(3.33) gives rise to the matrix equation(3.35)BYm=(χm+ℏ)BYm-1-ℏ(1-χm)R+ℏQm-1,
subject to the boundary conditions(3.36)∑k=0ND0kFm(ξk)=0,∑k=0NDNkFm(ξk)=0,Fm(ξN)=0,Θm(ξ0)=0,Θm(ξN)=0,Φm(ξ0)=0,Φm(ξN)=0,
where B and R are as defined in (3.25) and (3.37)Ym=[Fm(ξ0),Fm(ξ1),…,Fm(ξN),Θm(ξ0),Θm(ξ1),…,Θm(ξN),Φm(ξ0),Φm(ξ1),…,Φm(ξN)]T,Qm-1=(∑n=0m-1[-4L2(DFn)(DFm-1-n)+4L2mFn(D2Fm-1-n)]2LPr∑n=0m-1[-(DFn)Θm-1-n+m(DΘn)Fm-1-n]2LSc∑n=0m-1[-(DFn)Φm-1-n+m(DΦn)Fm-1-n]).
Applying the boundary conditions (3.32) on the right-hand side of (3.35) yields the following recursive formula for higher-order approximations Ym(ξ) for m≥1:(3.38)Ym=(χm+ℏ)B-1B̃Ym-1+ℏB-1[Qm-1-(1-χm)R].


### 3.3. Improved Spectral-Homotopy Analysis Method (ISHAM)

Details of the improved spectral-homotopy analysis method (ISHAM) can be found in Makukula et al. [30]. The main objective is to improve the convergence rate of the spectral-homotopy analysis method by using an optimal initial approximation. Hence, instead of a random solution choice a systematic approach is employed to find the optimal initial approximation. This is achieved by first assuming that the solutions f(η), θ(η),  and ϕ(η) can be expanded into(3.39)f(η)=Fi(η)+∑m=0i-1Fm(η),θ(η)=Θi(η)+∑m=0i-1Θm(η),ϕ(η)=Φi(η)+∑m=0i-1Φm(η),i=1,2,3,…,
where Fi, Θi,  and Φi are unknown functions whose solutions are obtained using the SHAM approach at the ith iteration and Fm, Θm,  and Φm (m≥1) are known from previous iterations. We use the same initial guesses as with the SHAM solution in Sections 3.1 and 3.2. Substituting (3.39) into the governing equations gives(3.40)Fi′′′+a1,i-1Fi′′+a2,i-1Fi′+a3,i-1Fi+mFi′′Fi-Fi′Fi′=r1,i-1(η),Θi′′+b1Θi′+b2Θi+b3Fi′+b4Fi-PrFi′Θi+mPrFiΘi′=r2,i-1(η),Φi′′+c1Φi′+c2Φi+c3Fi′+c4Fi-ScFi′Φi+mScFiΦi′=r3,i-1(η),
subject to the boundary conditions(3.41)Fi(0)=0,Fi′(0)=Fi′(∞)=0,Θi(0)=0,Θi(∞)=0,Φi(0)=0,Φi(∞)=0.
The coefficient parameters ak,i-1,  bk,i-1,  ck,i-1 (k=0,…,4), r1,i-1, r2,i-1,  and r3,i-1 are as defined in equation (3.4). Starting from the initial guesses (3.5), the subsequent solutions Fi, Θi,  and Φi (i≥1) are obtained by recursively solving (3.40) using the SHAM approach. To find the SHAM solutions of (3.40),  we start by defining the following linear operators:(3.42)LF[F̃i(η;q)]=∂3F̃i∂η3+a1,i-1∂2F̃i∂η2+a2,i-1∂F̃i∂η+a3,i-1F̃i,LΘ[F̃i(η;q),Θ̃i(η;q)]=∂2Θ̃i∂η2+b1,i-1∂Θ̃i∂η+b2,i-1Θ̃i+b3,i-1∂F̃i∂η+b4,i-1F̃i,LΦ[F̃i(η;q),Φ̃i(η;q)]=∂2Φ̃i∂η2+c1,i-1∂Φ̃i∂η+c2,i-1Φ̃i+c3,i-1∂F̃i∂η+c4,i-1F̃i.
The zeroth-order deformation equations are given by(3.43)(1-q)LF[F̃i(η;q)-Fi,0(η)]=qℏ{NF[F̃i(η;q)]-r1,i-1},(1-q)LΘ[Θ̃i(η;q)-Θi,0(η)]=qℏ{NΘ[F̃i(η;q),Θ̃i(η;q)]-r2,i-1},(1-q)LΦ[Φ̃i(η;q)-Φi,0(η)]=qℏ{NΦ[F̃i(η;q),Φ̃i(η;q)]-r3,i-1},𝒩F,  𝒩Θ, and 𝒩Φ are nonlinear operators given by(3.44)NF[F̃i(η;q)]=∂3F̃i∂η3+a1,i-1∂2F̃i∂η2+a2,i-1∂F̃i∂η+a3,i-1F̃i+mF̃i∂2F̃i∂η2-∂F̃i∂η∂F̃i∂η,NΘ[F̃i(η;q),Θ̃i(η;q)]=∂2Θ̃i∂η2+b1,i-1∂Θ̃i∂η+b2,i-1Θ̃i+b3,i-1∂F̃i∂η+b4,i-1F̃i+Pr(-Θ̃i∂F̃i∂η+mF̃i∂Θ̃i∂η),NΦ[F̃i(η;q),Φ̃i(η;q)]=∂2Φ̃i∂η2+c1,i-1∂Φ̃i∂η+c2,i-1Φ̃i+c3,i-1∂F̃i∂η+c4,i-1F̃i+Sc(-Φ̃i∂F̃i∂η+mF̃i∂Φ̃i∂η).
The mth order deformation equations are(3.45)LF[Fi,m(η)-χmFi,m-1(η)]=ℏRi,mF,LΘ[Θi,m(η)-χmΘi,m-1(η)]=ℏRi,mΘ,LΦ[Φi,m(η)-χmΦi,m-1(η)]=ℏRi,mΦ,
subject to the boundary conditions(3.46)Fi,m(0)=Fi,m′(0)=Fi,m′(∞)=0,Θi,m(0)=Θi,m(∞)=0,Φi,m(0)=Φi,m(∞)=0,
where(3.47)Ri,mF(η)=Fi,m-1′′′+a1,i-1Fi,m-1′′+a2,i-1Fi,m-1′+a3,i-1Fi,m-1+∑n=0m-1(-Fi,n′Fi,m-1-n′+mFi,nFi,m-1-n′′)-r1,i-1(η)(1-χm),Ri,mΘ(η)=Θi,m-1′′+b1,i-1Θi,m-1′+b2,i-1Θi,m-1+b3,i-1Fi,m-1′+b4,i-1Fi,m-1+Pr∑n=0m-1(-Fi,n′Θi,m-1-n+mΘi,n′Fi,m-1-n)-r2,i-1(η)(1-χm),Ri,mΦ(η)=Φi,m-1′′+c1,i-1Φi,m-1′+c2,i-1Φi,m-1+c3,i-1Fi,m-1′+c4,i-1Fi,m-1+Sc∑n=0m-1(-Fi,n′Φi,m-1-n+mΦi,n′Fi,m-1-n)-r3,i-1(η)(1-χm).
The initial approximations Fi,0, Θi,0, and Φi,0 that are used in the higher-order equations (3.45)–(3.47) are obtained by solving the linear part of (3.40) given by(3.48)Fi,0′′′+a1,i-1Fi,0′′+a2,i-1Fi,0′+a3,i-1Fi,0=r1,i-1,Θi,0′′+b1,i-1Θi,0′+b2,i-1Θi,0+b3,i-1Fi,0′+b4,i-1Fi,0=r2,i-1,Φi,0′′+c1,i-1Φi,0′+c2,i-1Φi,0+c3,i-1Fi,0′+c4,i-1Fi,0=r3,i-1,
with the boundary conditions(3.49)Fi,0(0)=Fi,0′(0)=Fi,0′(∞)=0,Θi,0(0)=0,Θi,0(∞)=0,Φi,0(0)=0,Φi,0(∞)=0.
In a similar manner, we apply the spectral methods to solve for the initial approximate solutions Fi,0, Θi,0,  and Φi,0,  and the higher-order deformation equations (3.45)–(3.47) for higher order approximate solutions Fi,m, Θi,m,  and Φi,m for m≥1. The solutions for Fi, Θi, and Φi are then generated using the solutions for Fi,m, Θi,m,  and Φi,m as follows:(3.50)Fi=Fi,0+Fi,1+Fi,2+Fi,3+⋯+Fi,m,Θi=Θi,0+Θi,1+Θi,2+Θi,3+⋯+Θi,m,Φi=Φi,0+Φi,1+Φi,2+Φi,3+⋯+Φi,m.
The [i,m] approximate solutions for f(η), θ(η), and ϕ(η) are then obtained by substituting Fi, Θi, and Φi from (3.50) into (3.39), where i is the ith iteration of the higher-order deformation equation and m is the mth iteration of the initial approximation.


## 3.1. The Successive Linearisation Method (SLM)

The successive linearisation method (see Makukula et al. [27, 28]) is used to solve (2.4)–(2.7). The starting point is to assume that the independent variables f(η), θ(η), and ϕ(η) may be expanded as(3.2)f(η)=fi(η)+∑m=0i-1Fm(η),θ(η)=θi(η)+∑m=0i-1Θm(η),ϕ(η)=ϕi(η)+∑m=0i-1Φm(η),i=1,2,3,…,
where fi, θi and ϕi, are unknown functions and Fm, Θm,  and Φm(m≥1) are approximations that are obtained by recursively solving the linear part of the equation that results from substituting (3.2) in the governing equations (2.4)–(2.7). Substituting (3.2) in the governing equations (2.4)–(2.7) gives    (3.3)fi′′′+a1,i-1fi′′+a2,i-1fi′+a3,i-1fi+mff′′-f′2=r1,i-1,θi′′+b1,i-1θi′+b2,i-1θi+b3,i-1fi′+b4,i-1fi-Prf′θ+mPrfθ′=r2,i-1,ϕi′′+c1,i-1ϕi′+c2,i-1ϕi+c3,i-1fi′+c4,i-1fi-Scf′ϕ+mScfϕ′=r3,i-1,
where the coefficient parameters ak,i-1,  bk,i-1, ck,i-1 (k=1,…,4), r1,i-1, r2,i-1,  and r3,i-1 are defined as(3.4)a1,i-1=m∑m=0i-1Fm,a2,i-1=-(2∑m=0i-1Fm′+M2+λPr),a3,i-1=m∑m=0i-1Fm′′,b1,i-1=mPr∑m=0i-1Fm,b2,i-1=-Pr∑m=0i-1Fm′,b3,i-1=-Pr∑m=0i-1Θm,b4,i-1=mPr∑m=0i-1Θm′,c1,i-1=mSc∑m=0i-1Fm,c2,i-1=-Sc∑m=0i-1Fm′-Scγ,c3,i-1=-Sc∑m=0i-1Φm,c4,i-1=mSc∑m=0i-1Φm′,r1,i-1=-[∑m=0i-1Fm′′′-∑m=0i-1Fm′∑m=0i-1Fm′+m∑m=0i-1Fm′′∑m=0i-1Fm-(M2+λPr)∑m=0i-1Fm′],r2,i-1=-[∑m=0i-1Θm′′-Pr∑m=0i-1Fm′∑m=0i-1Θm+mPr∑m=0i-1Fm∑m=0i-1Θm′],r3,i-1=-[∑m=0i-1Φm′′-Sc∑m=0i-1Fm′∑m=0i-1Φm+mSc∑m=0i-1Fm∑m=0i-1Φm′-Scγ∑m=0i-1Φ].Starting from the initial approximations(3.5)F0(η)=s+e-2η-e-η,Θ0(η)=e-η,Φ0(η)=e-η,
which are chosen to satisfy the boundary conditions (2.7), the subsequent solutions for Fm, Θm,  and Φm, m≥1,  are obtained by successively solving the linearized form of (3.3) which are(3.6)Fi′′′+a1,i-1Fi′′+a2,i-1Fi′+a3,i-1Fi=r1,i-1,Θi′′+b1,i-1Θi′+b2,i-1Θi+b3,i-1Fi′+b4,i-1Fi=r2,i-1,Φi′′+c1,i-1Φi′+c2,i-1Φi+c3,i-1Fi′+c4,i-1Fi=r3,i-1,
subject to the boundary conditions(3.7)Fi(0)=Fi′(0)=Fi′(∞)=0,Θi(0)=0,Φi(0)=0.
Once each solution for Fi, Θi,  and Φi (i≥1) has been found from iteratively solving (3.6)-(3.7) for each i, the approximate solutions for f(η), θ(η),  and ϕ(η) are obtained as(3.8)f(η)≈∑m=0MFm(η),θ(η)≈∑m=0MΘm(η),ϕ(η)≈∑m=0MΦm(η).
In coming up with (3.8), we have assumed that (3.9)lim⁡i→∞fi=0,lim⁡i→∞θi=0,lim⁡i→∞ϕi=0.Equations (3.6)-(3.7) are integrated using the Chebyshev spectral collocation method (Canuto et al. [33] and Trefethen [34]). The unknown functions are defined by the Chebyshev interpolating polynomials with the Gauss-Lobatto points defined as(3.10)yj=cos⁡πjN,j=0,1,…,N,
where N is the number of collocation points used. The unknown functions Fi, Θi,  and Φi are approximated at the collocation points by(3.11)Fi(ξ)≈∑k=0NFi(ξk)Tk(ξj),Θi(ξ)≈∑k=0NΘi(ξk)Tk(ξj),Φi(ξ)≈∑k=0NΦi(ξk)Tk(ξj),j=0,1,…,N,
where Tk is the kth Chebyshev polynomial defined as(3.12)Tk(ξ)=cos⁡[kcos⁡-1(ξ)].
The derivatives at the collocation points are represented as(3.13)daFidηa=∑k=0NDkjaFi(ξk),daΘidηa=∑k=0NDkjaΘi(ξk),daΦidηa=∑k=0NDkjaΦi(ξk),j=0,1,…,N,
where a is the order of differentiation and D=(2/L)𝒟 with 𝒟 being the Chebyshev spectral differentiation matrix. Substituting (3.13) in (3.6)-(3.7) leads to the matrix equation(3.14)Ai-1Xi=Ri-1,
where Bi-1 is a 3(N+1)×3(N+1) square matrix and Xi and Pi-1 are 3(N+1)×1 column vectors defined by(3.15)Ai-1=(D3+a1D2+a2D+a3I0I0Ib3D+b4ID2+b1D+b2I0Ic3D+c4I0ID2+c1D+c2I),Xi=[Fm(ξ0),Fm(ξ1),…,Fm(ξN),Θm(ξ0),Θm(ξ1),…,Θm(ξN),Φm(ξ0),Φm(ξ1),…,Φm(ξN)]T,Ri-1=[r1,i-1(η0),r1,i-1(η1),…,r1,i-1(ηN),r2,i-1(η0),r2,i-1(η1),…,r2,i-1(ηN),r3,i-1(η0),r3,i-1(η1),…,r3,i-1(ηN)]T.In the above definitions,ak,i-1, and bk,i-1, ck,i-1 (k=1,…,4) are diagonal matrices of size (N+1)×(N+1). After modifying the matrix system (3.14) to incorporate the boundary conditions, the solution is obtained as(3.16)Xi=Ai-1-1Ri-1.


## 3.2. Spectral-Homotopy Analysis Method (SHAM)

The spectral-homotopy analysis method (SHAM) has been used by Motsa et al. [25, 26]. It is also convenient to first ensure that the boundary conditions are made homogeneous by using the transformations(3.17)f(η)=F(ξ)+f0(η),θ(η)=Θ(ξ)+θ0(η),ϕ(η)=Φ(ξ)+ϕ0(η),
where f0(η), and  θ0(η), ϕ0(η) are chosen to satisfy the boundary conditions (2.7) of the governing equations (2.4)–(2.6). From (3.1) and the chain rule, we have that(3.18)f′(η)=2LF′(ξ)+f0′(η),f′′(η)=4L2F′′(ξ)+f0′′(η),f′′′(η)=8L3F′′′(ξ)+f0′′′(η),θ′(η)=2LΘ′(ξ)+θ0′(η),θ′′(η)=4L2Θ′′(ξ)+θ0′′(η),ϕ′(η)=2LΦ′(ξ)+ϕ0′(η),ϕ′′(η)=4L2Φ′′(ξ)+ϕ0′′(η).
Substituting (3.1) and (3.17)-(3.18) in the governing equations and boundary conditions gives(3.19)a0F′′′+a1F′′+a2F′+a3F+4L2mF′′F-4L2F′F′=r1(η),b0Θ′′+b1Θ′+b2Θ+b3F′+b4F-2LPrF′Θ+2LmPrFΘ′=r2(η),c0Φ′′+c1Φ′+c2Φ+c3F′+c4F-2LScF′Φ+2LmScFΦ′=r3(η),
where prime now denotes derivative with respect to ξ and (3.20)a0=8L3,a1=4L2mf0,a2=-(2L(M2+λPr)+4Lf0′),a3=mf0′′,r1(η)=-(f0′′′+mf0f0′′-f0′f0′-(M2+λPr)f0′),b0=4L2,b1=2LmPrf0,b2=-Prf0′,b3=-2LPrθ0,b4=mPrθ0′,r2(η)=-(θ′′-Prf0′θ0+mPrf0θ0),c0=4L2,c1=2LmScf0,c2=-Sc(f0′+γ),c3=-2LScϕ0,c4=mScϕ0′,r3(η)=-(ϕ′′-Scf0′ϕ0+mScf0ϕ0′-Scγϕ).The initial guesses used are(3.21)f0(η)=s+e-2η-ηe-η,θ0(η)=e-η,ϕ0(η)=e-η.
Solving the linear part of the equation system (3.19), that is,(3.22)a0F0′′′+a1F0′′+a2F0′+a3F0=r1(η),b0Θ0′′+b1Θ0′+b2Θ0+b3F0′+b4F0=r2(η),c0Φ0′′+c1Φ0′+c2Φ0+c3F0′+c4F0=r3(η),
subject to(3.23)F0(-1)=2LF0′(-1)=2LF0′(1)=0,Θ0(-1)=Θ0(1)=0,Φ0(-1)=Φ0(1)=0,
will yield the initial SHAM approximate solution. Applying the Chebyshev pseudospectral method on equations (3.22)-(3.23) yields the matrix form(3.24)BY0=R,
where(3.25)B=(a0D3+a1D2+a2D+a3I0I0Ib3D+b4Ib0D2+b1D+b2I0Ic3D+c4I0Ic0D2+c1D+c2I),(3.26)R=[r1(η0),r1(η1),…,r1(ηN),r2(η0),r2(η1),…,r2(ηN),r3(η0),r3(η1),…,r3(ηN)]T,Y0=[F0(ξ0),F0(ξ1),…,F0(ξN),Θ0(ξ0),Θ0(ξ1),…,Θ0(ξN),Φ0(ξ0),Φ0(ξ1),…,Φ0(ξN)]T,ai=diag⁡([ai(η0),…,ai(ηN-1),ai(ηN)]),bi=diag⁡([bi(η0),…,bi(ηN-1),bi(ηN)]),ci=diag⁡([ci(η0),…,ci(ηN-1),ci(ηN)]),i=0,1,2,3,4.The superscript T denotes the transpose, diag is a diagonal matrix, and I is an identity matrix of size (N+1)×(N+1). The boundary conditions (3.23) are implemented in matrix B and vector R of equation (3.24). The values of [Y0(ξ1),Y0(ξ2),…,Y0(ξN-1)] are then determined from the  following equation:(3.27)Y0=B-1R,
which provides us with the initial approximation for the solution of the governing equations (3.19). With the initial approximate solution, we then find approximate solutions for the nonlinear equations (3.19). We start by defining the following linear operators:(3.28)LF[F̃(ξ;q)]=a0∂3F̃∂ξ3+a1∂2F̃∂ξ2+a2∂F̃∂ξ+a3F̃,LΘ[F̃(ξ;q),Θ̃(ξ;q)]=b0∂2Θ̃∂ξ2+b1∂Θ̃∂ξ+b2Θ̃+b3∂F̃∂ξ+b4F̃,LΦ[F̃(ξ;q),Φ̃(ξ;q)]=c0∂2Φ̃∂ξ2+c1∂Φ̃∂ξ+c2Φ̃+c3∂F̃∂ξ+c4F̃,
where q∈[0,1] is the embedding parameter and F̃(ξ;q), Θ̃(ξ;q), and Φ̃(ξ;q) are unknown functions. The zeroth-order deformation equations are given by(3.29)(1-q)LF[F̃(ξ;q)-F0(ξ)]=qℏ{NF[F̃(ξ;q)]-r1},(1-q)LΘ[Θ̃(ξ;q)-Θ0(ξ)]=qℏ{NΘ[F̃(ξ;q),Θ̃(ξ;q)]-r2},(1-q)LΦ[Φ̃(ξ;q)-Φ0(ξ)]=qℏ{NΦ[F̃(ξ;q),Φ̃(ξ;q)]-r3},
where ℏ is the nonzero convergence controlling auxiliary parameter and 𝒩F, 𝒩Θ,  and 𝒩Φ are nonlinear operators given by(3.30)NF[F̃(ξ;q)]=a0∂3F̃∂ξ3+a1∂2F̃∂ξ2+a2∂F̃∂ξ+a3F̃+4L2mF̃∂2F̃∂ξ2-4L2∂F̃∂ξ∂F̃∂ξ,NΘ[F̃(ξ;q),Θ̃(ξ;q)]=b0∂2Θ̃∂ξ2+b1∂Θ̃∂ξ+b2Θ̃+b3∂F̃∂ξ+b4F̃+2LPr(-Θ̃∂F̃∂ξ+mF̃∂Θ̃∂ξ),NΦ[F̃(ξ;q),Φ̃(ξ;q)]=c0∂2Φ̃∂ξ2+c1∂Φ̃∂ξ+c2Φ̃+c3∂F̃∂ξ+c4F̃+2LSc(-Φ̃∂F̃∂ξ+mF̃∂Φ̃∂ξ).
The m-th order deformation equations are given by(3.31)LF[Fm(ξ)-χmFm-1(ξ)]=ℏRmF,LΘ[Θm(ξ)-χmΘm-1(ξ)]=ℏRmΘ,LΦ[Φm(ξ)-χmΦm-1(ξ)]=ℏRmΦ,
subject to the boundary conditions(3.32)Fm(-1)=Fm′(-1)=Fm′(1)=0,Θm(-1)=Θm(1)=0,Φm(-1)=Φm(1)=0,
where(3.33)RmF(ξ)=a0Fm-1′′′+a1Fm-1′′+a2Fm-1′+a3Fm-1+4L2∑n=0m-1(-Fn′Fm-1-n′+mFnFm-1-n′′)-r1(η)(1-χm),RmΘ(ξ)=b0Θm-1′′+b1Θm-1′+b2Θm-1+b3Fm-1′+b4Fm-1+2LPr∑n=0m-1(-Fn′Θm-1-n+mΘn′Fm-1-n)-r2(η)(1-χm),RmΦ(ξ)=c0Φm-1′′+c1Φm-1′+c2Φm-1+c3Fm-1′+c4Fm-1+2LSc∑n=0m-1(-Fn′Φm-1-n+mΦn′Fm-1-n)-r3(η)(1-χm),(3.34)χm={0,m≤11,m>1.
Applying the Chebyshev pseudospectral transformation to equations (3.31)–(3.33) gives rise to the matrix equation(3.35)BYm=(χm+ℏ)BYm-1-ℏ(1-χm)R+ℏQm-1,
subject to the boundary conditions(3.36)∑k=0ND0kFm(ξk)=0,∑k=0NDNkFm(ξk)=0,Fm(ξN)=0,Θm(ξ0)=0,Θm(ξN)=0,Φm(ξ0)=0,Φm(ξN)=0,
where B and R are as defined in (3.25) and (3.37)Ym=[Fm(ξ0),Fm(ξ1),…,Fm(ξN),Θm(ξ0),Θm(ξ1),…,Θm(ξN),Φm(ξ0),Φm(ξ1),…,Φm(ξN)]T,Qm-1=(∑n=0m-1[-4L2(DFn)(DFm-1-n)+4L2mFn(D2Fm-1-n)]2LPr∑n=0m-1[-(DFn)Θm-1-n+m(DΘn)Fm-1-n]2LSc∑n=0m-1[-(DFn)Φm-1-n+m(DΦn)Fm-1-n]).
Applying the boundary conditions (3.32) on the right-hand side of (3.35) yields the following recursive formula for higher-order approximations Ym(ξ) for m≥1:(3.38)Ym=(χm+ℏ)B-1B̃Ym-1+ℏB-1[Qm-1-(1-χm)R].


## 3.3. Improved Spectral-Homotopy Analysis Method (ISHAM)

Details of the improved spectral-homotopy analysis method (ISHAM) can be found in Makukula et al. [30]. The main objective is to improve the convergence rate of the spectral-homotopy analysis method by using an optimal initial approximation. Hence, instead of a random solution choice a systematic approach is employed to find the optimal initial approximation. This is achieved by first assuming that the solutions f(η), θ(η),  and ϕ(η) can be expanded into(3.39)f(η)=Fi(η)+∑m=0i-1Fm(η),θ(η)=Θi(η)+∑m=0i-1Θm(η),ϕ(η)=Φi(η)+∑m=0i-1Φm(η),i=1,2,3,…,
where Fi, Θi,  and Φi are unknown functions whose solutions are obtained using the SHAM approach at the ith iteration and Fm, Θm,  and Φm (m≥1) are known from previous iterations. We use the same initial guesses as with the SHAM solution in Sections 3.1 and 3.2. Substituting (3.39) into the governing equations gives(3.40)Fi′′′+a1,i-1Fi′′+a2,i-1Fi′+a3,i-1Fi+mFi′′Fi-Fi′Fi′=r1,i-1(η),Θi′′+b1Θi′+b2Θi+b3Fi′+b4Fi-PrFi′Θi+mPrFiΘi′=r2,i-1(η),Φi′′+c1Φi′+c2Φi+c3Fi′+c4Fi-ScFi′Φi+mScFiΦi′=r3,i-1(η),
subject to the boundary conditions(3.41)Fi(0)=0,Fi′(0)=Fi′(∞)=0,Θi(0)=0,Θi(∞)=0,Φi(0)=0,Φi(∞)=0.
The coefficient parameters ak,i-1,  bk,i-1,  ck,i-1 (k=0,…,4), r1,i-1, r2,i-1,  and r3,i-1 are as defined in equation (3.4). Starting from the initial guesses (3.5), the subsequent solutions Fi, Θi,  and Φi (i≥1) are obtained by recursively solving (3.40) using the SHAM approach. To find the SHAM solutions of (3.40),  we start by defining the following linear operators:(3.42)LF[F̃i(η;q)]=∂3F̃i∂η3+a1,i-1∂2F̃i∂η2+a2,i-1∂F̃i∂η+a3,i-1F̃i,LΘ[F̃i(η;q),Θ̃i(η;q)]=∂2Θ̃i∂η2+b1,i-1∂Θ̃i∂η+b2,i-1Θ̃i+b3,i-1∂F̃i∂η+b4,i-1F̃i,LΦ[F̃i(η;q),Φ̃i(η;q)]=∂2Φ̃i∂η2+c1,i-1∂Φ̃i∂η+c2,i-1Φ̃i+c3,i-1∂F̃i∂η+c4,i-1F̃i.
The zeroth-order deformation equations are given by(3.43)(1-q)LF[F̃i(η;q)-Fi,0(η)]=qℏ{NF[F̃i(η;q)]-r1,i-1},(1-q)LΘ[Θ̃i(η;q)-Θi,0(η)]=qℏ{NΘ[F̃i(η;q),Θ̃i(η;q)]-r2,i-1},(1-q)LΦ[Φ̃i(η;q)-Φi,0(η)]=qℏ{NΦ[F̃i(η;q),Φ̃i(η;q)]-r3,i-1},𝒩F,  𝒩Θ, and 𝒩Φ are nonlinear operators given by(3.44)NF[F̃i(η;q)]=∂3F̃i∂η3+a1,i-1∂2F̃i∂η2+a2,i-1∂F̃i∂η+a3,i-1F̃i+mF̃i∂2F̃i∂η2-∂F̃i∂η∂F̃i∂η,NΘ[F̃i(η;q),Θ̃i(η;q)]=∂2Θ̃i∂η2+b1,i-1∂Θ̃i∂η+b2,i-1Θ̃i+b3,i-1∂F̃i∂η+b4,i-1F̃i+Pr(-Θ̃i∂F̃i∂η+mF̃i∂Θ̃i∂η),NΦ[F̃i(η;q),Φ̃i(η;q)]=∂2Φ̃i∂η2+c1,i-1∂Φ̃i∂η+c2,i-1Φ̃i+c3,i-1∂F̃i∂η+c4,i-1F̃i+Sc(-Φ̃i∂F̃i∂η+mF̃i∂Φ̃i∂η).
The mth order deformation equations are(3.45)LF[Fi,m(η)-χmFi,m-1(η)]=ℏRi,mF,LΘ[Θi,m(η)-χmΘi,m-1(η)]=ℏRi,mΘ,LΦ[Φi,m(η)-χmΦi,m-1(η)]=ℏRi,mΦ,
subject to the boundary conditions(3.46)Fi,m(0)=Fi,m′(0)=Fi,m′(∞)=0,Θi,m(0)=Θi,m(∞)=0,Φi,m(0)=Φi,m(∞)=0,
where(3.47)Ri,mF(η)=Fi,m-1′′′+a1,i-1Fi,m-1′′+a2,i-1Fi,m-1′+a3,i-1Fi,m-1+∑n=0m-1(-Fi,n′Fi,m-1-n′+mFi,nFi,m-1-n′′)-r1,i-1(η)(1-χm),Ri,mΘ(η)=Θi,m-1′′+b1,i-1Θi,m-1′+b2,i-1Θi,m-1+b3,i-1Fi,m-1′+b4,i-1Fi,m-1+Pr∑n=0m-1(-Fi,n′Θi,m-1-n+mΘi,n′Fi,m-1-n)-r2,i-1(η)(1-χm),Ri,mΦ(η)=Φi,m-1′′+c1,i-1Φi,m-1′+c2,i-1Φi,m-1+c3,i-1Fi,m-1′+c4,i-1Fi,m-1+Sc∑n=0m-1(-Fi,n′Φi,m-1-n+mΦi,n′Fi,m-1-n)-r3,i-1(η)(1-χm).
The initial approximations Fi,0, Θi,0, and Φi,0 that are used in the higher-order equations (3.45)–(3.47) are obtained by solving the linear part of (3.40) given by(3.48)Fi,0′′′+a1,i-1Fi,0′′+a2,i-1Fi,0′+a3,i-1Fi,0=r1,i-1,Θi,0′′+b1,i-1Θi,0′+b2,i-1Θi,0+b3,i-1Fi,0′+b4,i-1Fi,0=r2,i-1,Φi,0′′+c1,i-1Φi,0′+c2,i-1Φi,0+c3,i-1Fi,0′+c4,i-1Fi,0=r3,i-1,
with the boundary conditions(3.49)Fi,0(0)=Fi,0′(0)=Fi,0′(∞)=0,Θi,0(0)=0,Θi,0(∞)=0,Φi,0(0)=0,Φi,0(∞)=0.
In a similar manner, we apply the spectral methods to solve for the initial approximate solutions Fi,0, Θi,0,  and Φi,0,  and the higher-order deformation equations (3.45)–(3.47) for higher order approximate solutions Fi,m, Θi,m,  and Φi,m for m≥1. The solutions for Fi, Θi, and Φi are then generated using the solutions for Fi,m, Θi,m,  and Φi,m as follows:(3.50)Fi=Fi,0+Fi,1+Fi,2+Fi,3+⋯+Fi,m,Θi=Θi,0+Θi,1+Θi,2+Θi,3+⋯+Θi,m,Φi=Φi,0+Φi,1+Φi,2+Φi,3+⋯+Φi,m.
The [i,m] approximate solutions for f(η), θ(η), and ϕ(η) are then obtained by substituting Fi, Θi, and Φi from (3.50) into (3.39), where i is the ith iteration of the higher-order deformation equation and m is the mth iteration of the initial approximation.


## 4. Results and Discussion

Equations (2.4)-(2.6) subject to boundary conditions (2.7) have been solved using three recent semi-numerical techniques as described above. To validate our results, we have compared the skin friction coefficient, the Nusselt number, and the Sherwood number with the theoretical results of Muhaimin et al. [9]. We have further compared our results with the full numerical solutions obtained using the shooting method and the Matlab bvp4c routine. The comparison is given in Tables 1–3.Table 1
Comparison of the approximate solutions off′′(0) at different orders of the SLM, SHAM,  and ISHAM against the numerical solutions at different values of λ when s = 3,  M = 1,    m=1,     Sc = 0.62,   γ =3,    Pr = 1,   λ = 0,    ℏ=−1,    L  = 30,  and  N = 150.
SLMSHAMISHAMShootingbvp4cRef. [9]λorderf′′(0)orderf′′(0)orderf′′(0)f′′(0)f′′(0)f′′(0)113.3206813.30709[1,1]3.338583.302783.302783.30277623.3028323.30338[2,2]3.3027833.3027843.30279[3,3]3.3027843.3027863.30278[4,4]3.30278213.5729213.56643[1,1]3.598832.302782.302783.56155323.5615723.56225[2,2]2.3027833.5615543.56157[3,3]2.3027843.5615563.56155[4,4]2.30278414.0077014.00549[1,1]4.037894.000004.000004.00000024.0000024.00084[2,2]4.0000034.0000044.00002[3,3]4.0000044.0000064.00000[4,4]4.00000Table 2
Comparison of the approximate solutions of-θ′(0) at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of λ when s = 3,  M = 1,  m = 1,  Sc = 0.62,  γ = 3,  Pr = 1,   λ = 0,  ℏ = −1,  L = 30,  and   N = 150.
SLMSHAMISHAMShootingbvp4cRef. [9]λorder-θ′(0)order-θ′(0)order-θ′(0)-θ′(0)-θ′(0)-θ′(0)112.5678312.67511[1,1]2.809122.665542.665542.66553722.6648522.66656[2,2]2.6655432.6655442.66556[3,3]2.6655442.6655462.66554[4,4]2.66554212.5919812.69041[1,1]2.840092.680322.680322.68031522.6798722.68132[2,2]2.6803232.6803242.68034[3,3]2.6803242.6803262.68032[4,4]2.68032412.6291412.71316[1,1]2.884522.702402.702402.70245522.7021522.70336[2,2]2.7024032.7024042.70241[3,3]2.7024042.7024062.70240[4,4]2.70240Table 3
Comparison of the approximate solutions of-ϕ′(0) at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of λ when s = 3,  M = 1,    m = 1,   Sc = 0.62,  γ = 3,    Pr = 1,   λ = 0,   ℏ = −1,    L = 30, and    N = 150.
SLMSHAMISHAMShootingbvp4cRef. [9]λorder-ϕ′(0)order-ϕ′(0)order-ϕ′(0)-ϕ′(0)-ϕ′(0)-ϕ′(0)112.3929412.41413[1,1]2.439762.410292.410292.41028322.4102622.41085[2,2]2.4102932.4102942.41030[3,3]2.4102942.4102962.41029[4,4]2.41029212.4018112.42137[1,1]2.449812.417002.417002.41700022.4169822.41764[2,2]2.4170032.4170042.41702[3,3]2.4170042.4170062.41700[4,4]2.41700412.4155912.43236[1,1]2.464602.427222.427222.42722522.4272122.42797[2,2]2.4272232.4272242.42724[3,3]2.4272242.4272262.42722[4,4]2.42722Tables1–3 give values of the skin friction, heat transfer rate, and the mass transfer rate,  respectively, for different porosity values. The convergence to the two numerical results of the SLM is achieved at the third order of approximation, at the sixth order for the SHAM, and at second order for the ISHAM. Comparison with results reported in Muhaimin et al. [9] shows an excellent agreement.Table1 shows an increase in the surface shear stress f′′(0) with an increase in the porosity parameter λ. The increase in the skin friction with the porosity may be accounted for by the fact that the velocity gradient increases with porosity  (Takhar et al. [35]). Tables 2 and 3 show an increase in the surface heat transfer rate -θ′(0) and the mass transfer rate -θ′(0) with the porosity parameter for large suction values (s=3), suggesting an increase in temperature and concentration gradients with increasing porosity.Figure1 serves two purposes: (a) to give sense of the accuracy of the improved spectral homotopy analysis (ISHAM) by means of a comparison between the numerical results and the second-order improved spectral-homotopy analysis results and (b) to demonstrate the effects of the suction parameter s and the Hartmann number M on the velocity profiles f′(η).On the comparison between the 2nd-order ISHAM solution (figures) and the bvp4c numerical solution (solid line) forf(η) and θ(η) at different values of λ when M=1, m=1,Pr=3,λ1=2,s=1,L=30, and   N=150.
(a)(b)Firstly we observe an excellent agreement between the second-order ISHAM and the numerical bvp4c results for all parameter values used. Secondly we note that these results are qualitatively similar to those reported in Noor et al. [6] for the case of one-direction shrinking (m=1) and show that increasing the suction parameter s and the Hartmann number M leads to an increase in the velocity. This in turn leads to a decrease in the boundary layer thickness as fluid is sucked out of the flow region.


## 5. Conclusions

We have successfully solved the nonlinear system of equations governing MHD boundary layer past a porous shrinking sheet with a chemical reaction and suction.  We demonstrated three recent innovative methods, namely, the successive linearisation method (SLM), the spectral-homotopy analysis method (SHAM), and the improved spectral-homotopy analysis method (ISHAM), and compared the performance of the three methods with regard to the speed of convergence of the solution (the number of iterations required), computational efficiency, and the ease of application of the method. The results were compared with those obtained using the well-known shooting method and the Matlab bvp4c solver. We found that the ISHAM converged at second order. The magnitude of the parameter values used did not affect its performance under the same conditions with the SLM and SHAM. Nevertheless, the ISHAM does not come cheap in terms of the size of the code and computer time, taking about three times as long as the SLM to compute the same result and about double the time taken with the SHAM. The SLM converged at third order, is easy to implement, and has shown a good level of stability when solving highly nonlinear problems. The SHAM gives good convergence under the same conditions but poor convergence with highly nonlinear problems. It is easy to implement but not as easy as with the SLM.Results from simulations revealed an excellent agreement between results from the shooting method and the bvp4c. Our findings indicate that the ISHAM is the best approach of the three methods in terms of the accuracy of the results and speed of convergence. Parametric studies for effects of different parameter values in the problems agreed with results present in the literature.



---

*Source: 489217-2011-11-01.xml*