About this Journal Submit a Manuscript Table of Contents
Mathematical Problems in Engineering
Volume 2011 (2011), Article ID 489217, 19 pages
http://dx.doi.org/10.1155/2011/489217
Research Article

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

1School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa
2Department of Mathematics and Applied Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa

Received 29 July 2011; Revised 22 August 2011; Accepted 23 August 2011

Academic Editor: Oleg V. Gendelman

Copyright © 2011 Z. G. Makukula et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 [13]. A shrinking film is useful in the packaging of bulk products since it can be unwrapped easily with adequate heat [47]. 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) [1417], and Homotopy perturbation method (HPM) [1823].

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. [2730]) 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 components 𝑢, 𝑣, and 𝑤 are (see Muhaimin et al. [5])𝜕𝑢+𝜕𝑥𝜕𝑣+𝜕𝑦𝜕𝑤𝑢𝜕𝑧=0,𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦+𝑤𝜕𝑢=1𝜕𝑧𝜌𝜕𝑝𝜕𝜕𝑥+𝜈2𝑢𝜕𝑥2+𝜕2𝑢𝜕𝑦2+𝜕2𝑢𝜕𝑧2𝑢𝜎𝐵20𝜌𝜈𝑢𝐾,𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦+𝑤𝜕𝑣=1𝜕𝑧𝜌𝜕𝑝𝜕𝜕𝑦+𝜈2𝑣𝜕𝑥2+𝜕2𝑣+𝜕𝜕𝑦2𝑣𝜕𝑧2𝑣𝜎𝐵20𝜌𝜈𝑣𝐾,𝑢𝜕𝑤𝜕𝑥+𝑣𝜕𝑤𝜕𝑦+𝑤𝜕𝑤=1𝜕𝑧𝜌𝜕𝑝𝜕𝜕𝑧+𝜈2𝑤𝜕𝑥2+𝜕2𝑤𝜕𝑦2+𝜕2𝑤𝜕𝑧2,𝑢𝜕𝑇𝜕𝑥+𝑣𝜕𝑇𝜕𝑦+𝑤𝜕𝑇𝜕𝜕𝑧=𝛼2𝑇𝜕𝑥2+𝜕2𝑇𝜕𝑦2+𝜕2𝑇𝜕𝑧2,𝑢𝜕𝐶𝜕𝑥+𝑣𝜕𝐶𝜕𝑦+𝑤𝜕𝐶𝜕𝜕𝑧=𝐷2𝐶𝜕𝑥2+𝜕2𝐶𝜕𝑦2+𝜕2𝐶𝜕𝑧2𝑘1𝐶,(2.1)

where 𝛼 is the thermal conductivity of the fluid, B0 is the magnetic field, 𝜅 is the thermal viscosity, 𝐾 is the permeability of the porous medium, 𝑘1 is the rate of chemical reaction, 𝜈=𝜇/𝜌 is the kinetic viscosity, 𝜇 is the dynamic viscosity, and 𝜎 is the electrical conductivity.

The applicable boundary conditions are𝑢=𝑎𝑥,𝑣=𝑎(𝑚1)𝑦,𝑤=𝑊,𝑇=𝑇𝑤,𝐶=𝐶𝑤,at𝑦=0,𝑢0,𝐶𝐶,𝑇𝑇as𝑦,(2.2) where 𝑎>0 is the shrinking constant and 𝑊 is the suction velocity. The cases 𝑚=1 and 𝑚=2 correspond to shrinking sheets in the 𝑥- and 𝑦-directions, respectively.

Using the similarity transformations (see Sajid and Hayat [32]): 𝑢=𝑎𝑥𝐹(𝜂),𝑣=𝑎(𝑚1)𝑦𝐹(𝜂),𝑤=𝑎𝜈𝑚𝐹(𝜂),𝜂=𝑎𝜈𝑧,𝜃(𝜂)=𝑇𝑇𝑇𝑤𝑇,𝜙(𝜂)=𝐶𝐶𝐶𝑤𝐶,(2.3)(2.1) are transformed to the system of nonlinear equations𝑓𝑀2+𝜆𝑃𝑟𝑓𝑓2+𝑚𝑓𝑓𝜃=0,(2.4)𝑃𝑟𝑓𝜃+𝑚𝑃𝑟𝑓𝜃𝜙=0,(2.5)𝑆𝑐𝑓𝜙+𝑚𝑆𝑐𝑓𝜙𝑆𝑐𝛾𝜙=0,(2.6) subject to𝑓(0)=𝑠,𝑓(0)=1,𝑓()=0,𝜃(0)=1,𝜃()=0,𝜙(0)=1,𝜙()=0,(2.7) where 𝑃𝑟=𝜈/𝜅 is the Prandtl number, 𝑆𝑐=𝜈/𝐷 is the Schmidt number, 𝜆=𝜅/𝐾𝑎 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 𝑓, but the solutions for 𝜃 and 𝜙 directly depend on the solution for 𝑓. 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 𝑃𝑟 and 𝑚 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 106.

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,𝐿], where 𝐿 is a fixed length that is taken to be larger than the thickness of the boundary layer. The domain [0,𝐿] is then transformed to [1,1] using the algebraic mapping𝜉=2𝜂𝐿[].1,𝜉1,1(3.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 𝑓(𝜂), 𝜃(𝜂), and 𝜙(𝜂) may be expanded as𝑓(𝜂)=𝑓𝑖(𝜂)+𝑖1𝑚=0𝐹𝑚(𝜂),𝜃(𝜂)=𝜃𝑖(𝜂)+𝑖1𝑚=0Θ𝑚(𝜂),𝜙(𝜂)=𝜙𝑖(𝜂)+𝑖1𝑚=0Φ𝑚(𝜂),𝑖=1,2,3,,(3.2) where 𝑓𝑖, 𝜃𝑖 and 𝜙𝑖, are unknown functions and 𝐹𝑚, Θ𝑚, and Φ𝑚(𝑚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 𝑓𝑖+𝑎1,𝑖1𝑓𝑖+𝑎2,𝑖1𝑓𝑖+𝑎3,𝑖1𝑓𝑖+𝑚𝑓𝑓𝑓2=𝑟1,𝑖1,𝜃𝑖+𝑏1,𝑖1𝜃𝑖+𝑏2,𝑖1𝜃𝑖+𝑏3,𝑖1𝑓𝑖+𝑏4,𝑖1𝑓𝑖𝑃𝑟𝑓𝜃+𝑚𝑃𝑟𝑓𝜃=𝑟2,𝑖1,𝜙𝑖+𝑐1,𝑖1𝜙𝑖+𝑐2,𝑖1𝜙𝑖+𝑐3,𝑖1𝑓𝑖+𝑐4,𝑖1𝑓𝑖𝑆𝑐𝑓𝜙+𝑚𝑆𝑐𝑓𝜙=𝑟3,𝑖1,(3.3) where the coefficient parameters 𝑎𝑘,𝑖1, 𝑏𝑘,𝑖1, 𝑐𝑘,𝑖1 (𝑘=1,,4), 𝑟1,𝑖1, 𝑟2,𝑖1, and 𝑟3,𝑖1 are defined as𝑎1,𝑖1=𝑚𝑖1𝑚=0𝐹𝑚,𝑎2,𝑖12=𝑖1𝑚=0𝐹𝑚+𝑀2+𝜆𝑃𝑟,𝑎3,𝑖1=𝑚𝑖1𝑚=0𝐹𝑚,𝑏1,𝑖1=𝑚𝑃𝑟𝑖1𝑚=0𝐹𝑚,𝑏2,𝑖1=𝑃𝑟𝑖1𝑚=0𝐹𝑚,𝑏3,𝑖1=𝑃𝑟𝑖1𝑚=0Θ𝑚,𝑏4,𝑖1=𝑚𝑃𝑟𝑖1𝑚=0Θ𝑚,𝑐1,𝑖1=𝑚𝑆𝑐𝑖1𝑚=0𝐹𝑚,𝑐2,𝑖1=𝑆𝑐𝑖1𝑚=0𝐹𝑚𝑐𝑆𝑐𝛾,3,𝑖1=𝑆𝑐𝑖1𝑚=0Φ𝑚,𝑐4,𝑖1=𝑚𝑆𝑐𝑖1𝑚=0Φ𝑚,𝑟1,𝑖1=𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0𝐹𝑚+𝑚𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0𝐹𝑚𝑀2+𝜆𝑃𝑟𝑖1𝑚=0𝐹𝑚,𝑟2,𝑖1=𝑖1𝑚=0Θ𝑚𝑃𝑟𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0Θ𝑚+𝑚𝑃𝑟𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0Θ𝑚,𝑟3,𝑖1=𝑖1𝑚=0Φ𝑚𝑆𝑐𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0Φ𝑚+𝑚𝑆𝑐𝑖1𝑚=0𝐹𝑚𝑖1𝑚=0Φ𝑚𝑆𝑐𝛾𝑖1𝑚=0Φ.(3.4)

Starting from the initial approximations𝐹0(𝜂)=𝑠+𝑒2𝜂𝑒𝜂,Θ0(𝜂)=𝑒𝜂,Φ0(𝜂)=𝑒𝜂,(3.5) which are chosen to satisfy the boundary conditions (2.7), the subsequent solutions for 𝐹𝑚, Θ𝑚, and Φ𝑚, 𝑚1, are obtained by successively solving the linearized form of (3.3) which are𝐹𝑖+𝑎1,𝑖1𝐹𝑖+𝑎2,𝑖1𝐹𝑖+𝑎3,𝑖1𝐹𝑖=𝑟1,𝑖1,Θ𝑖+𝑏1,𝑖1Θ𝑖+𝑏2,𝑖1Θ𝑖+𝑏3,𝑖1𝐹𝑖+𝑏4,𝑖1𝐹𝑖=𝑟2,𝑖1,Φ𝑖+𝑐1,𝑖1Φ𝑖+𝑐2,𝑖1Φ𝑖+𝑐3,𝑖1𝐹𝑖+𝑐4,𝑖1𝐹𝑖=𝑟3,𝑖1,(3.6) subject to the boundary conditions𝐹𝑖(0)=𝐹𝑖(0)=𝐹𝑖()=0,Θ𝑖(0)=0,Φ𝑖(0)=0.(3.7) Once each solution for 𝐹𝑖, Θ𝑖, and Φ𝑖 (𝑖1) has been found from iteratively solving (3.6)-(3.7) for each 𝑖, the approximate solutions for 𝑓(𝜂), 𝜃(𝜂), and 𝜙(𝜂) are obtained as𝑓(𝜂)𝑀𝑚=0𝐹𝑚(𝜂),𝜃(𝜂)𝑀𝑚=0Θ𝑚(𝜂),𝜙(𝜂)𝑀𝑚=0Φ𝑚(𝜂).(3.8) In coming up with (3.8), we have assumed that lim𝑖𝑓𝑖=0,lim𝑖𝜃𝑖=0,lim𝑖𝜙𝑖=0.(3.9)

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𝑦𝑗=cos𝜋𝑗𝑁,𝑗=0,1,,𝑁,(3.10) where 𝑁 is the number of collocation points used. The unknown functions 𝐹𝑖, Θ𝑖, and Φ𝑖 are approximated at the collocation points by𝐹𝑖(𝜉)𝑁𝑘=0𝐹𝑖𝜉𝑘𝑇𝑘𝜉𝑗,Θ𝑖(𝜉)𝑁𝑘=0Θ𝑖𝜉𝑘𝑇𝑘𝜉𝑗,Φ𝑖(𝜉)𝑁𝑘=0Φ𝑖𝜉𝑘𝑇𝑘𝜉𝑗,𝑗=0,1,,𝑁,(3.11) where 𝑇𝑘 is the 𝑘th Chebyshev polynomial defined as𝑇𝑘(𝜉)=cos𝑘cos1.(𝜉)(3.12) The derivatives at the collocation points are represented as𝑑𝑎𝐹𝑖𝑑𝜂𝑎=𝑁𝑘=0𝐃𝑎𝑘𝑗𝐹𝑖𝜉𝑘,𝑑𝑎Θ𝑖𝑑𝜂𝑎=𝑁𝑘=0𝐃𝑎𝑘𝑗Θ𝑖𝜉𝑘,𝑑𝑎Φ𝑖𝑑𝜂𝑎=𝑁𝑘=0𝐃𝑎𝑘𝑗Φ𝑖𝜉𝑘,𝑗=0,1,,𝑁,(3.13) where 𝑎 is the order of differentiation and 𝐃=(2/𝐿)𝒟 with 𝒟 being the Chebyshev spectral differentiation matrix. Substituting (3.13) in (3.6)-(3.7) leads to the matrix equation𝐀𝑖1𝐗𝑖=𝐑𝑖1,(3.14) where 𝐁𝑖1 is a 3(𝑁+1)×3(𝑁+1) square matrix and 𝐗𝑖 and 𝐏𝑖1 are 3(𝑁+1)×1 column vectors defined by𝐀𝑖1=𝐃3+𝐚1𝐃2+𝐚2𝐃+𝐚3𝐛𝐈0𝐈0𝐈3𝐃+𝐛4𝐈𝐃2+𝐛1𝐃+𝐛2𝐜𝐈0𝐈3𝐃+𝐜4𝐈0𝐈𝐃2+𝐜1𝐃+𝐜2𝐈,𝐗𝑖=𝐹𝑚𝜉0,𝐹𝑚𝜉1,,𝐹𝑚𝜉𝑁,Θ𝑚𝜉0,Θ𝑚𝜉1,,Θ𝑚𝜉𝑁,Φ𝑚𝜉0,Φ𝑚𝜉1,,Φ𝑚𝜉𝑁𝑇,𝐑𝑖1=𝑟1,𝑖1𝜂0,𝑟1,𝑖1𝜂1,,𝑟1,𝑖1𝜂𝑁,𝑟2,𝑖1𝜂0,𝑟2,𝑖1𝜂1,,𝑟2,𝑖1𝜂𝑁,𝑟3,𝑖1𝜂0,𝑟3,𝑖1𝜂1,,𝑟3,𝑖1𝜂𝑁𝑇.(3.15)

In the above definitions, 𝐚𝑘,𝑖1, and 𝐛𝑘,𝑖1, 𝐜𝑘,𝑖1 (𝑘=1,,4) are diagonal matrices of size (𝑁+1)×(𝑁+1). After modifying the matrix system (3.14) to incorporate the boundary conditions, the solution is obtained as𝐗𝑖=𝐀1𝑖1𝐑𝑖1.(3.16)

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𝑓(𝜂)=𝐹(𝜉)+𝑓0(𝜂),𝜃(𝜂)=Θ(𝜉)+𝜃0(𝜂),𝜙(𝜂)=Φ(𝜉)+𝜙0(𝜂),(3.17) where 𝑓0(𝜂), 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𝑓2(𝜂)=𝐿𝐹(𝜉)+𝑓0(𝜂),𝑓4(𝜂)=𝐿2𝐹(𝜉)+𝑓0𝑓(𝜂),8(𝜂)=𝐿3𝐹(𝜉)+𝑓0𝜃(𝜂),2(𝜂)=𝐿Θ(𝜉)+𝜃0(𝜂),𝜃4(𝜂)=𝐿2Θ(𝜉)+𝜃0𝜙(𝜂),2(𝜂)=𝐿Φ(𝜉)+𝜙0(𝜂),𝜙4(𝜂)=𝐿2Φ(𝜉)+𝜙0(𝜂).(3.18) Substituting (3.1) and (3.17)-(3.18) in the governing equations and boundary conditions gives𝑎0𝐹+𝑎1𝐹+𝑎2𝐹+𝑎34𝐹+𝐿2𝑚𝐹4𝐹𝐿2𝐹𝐹=𝑟1𝑏(𝜂),0Θ+𝑏1Θ+𝑏2Θ+𝑏3𝐹+𝑏42𝐹𝐿𝑃𝑟𝐹2Θ+𝐿𝑚𝑃𝑟𝐹Θ=𝑟2𝑐(𝜂),0Φ+𝑐1Φ+𝑐2Φ+𝑐3𝐹+𝑐42𝐹𝐿𝑆𝑐𝐹2Φ+𝐿𝑚𝑆𝑐𝐹Φ=𝑟3(𝜂),(3.19) where prime now denotes derivative with respect to 𝜉 and 𝑎0=8𝐿3,𝑎1=4𝐿2𝑚𝑓0,𝑎22=𝐿𝑀2+4+𝜆𝑃𝑟𝐿𝑓0,𝑎3=𝑚𝑓0,𝑟1𝑓(𝜂)=0+𝑚𝑓0𝑓0𝑓0𝑓0𝑀2𝑓+𝜆𝑃𝑟0,𝑏0=4𝐿2,𝑏1=2𝐿𝑚𝑃𝑟𝑓0,𝑏2=𝑃𝑟𝑓0,𝑏32=𝐿𝑃𝑟𝜃0,𝑏4=𝑚𝑃𝑟𝜃0,𝑟2𝜃(𝜂)=𝑃𝑟𝑓0𝜃0+𝑚𝑃𝑟𝑓0𝜃0,𝑐0=4𝐿2,𝑐1=2𝐿𝑚𝑆𝑐𝑓0,𝑐2𝑓=𝑆𝑐0+𝛾,𝑐32=𝐿𝑆𝑐𝜙0,𝑐4=𝑚𝑆𝑐𝜙0,𝑟3𝜙(𝜂)=𝑆𝑐𝑓0𝜙0+𝑚𝑆𝑐𝑓0𝜙0.𝑆𝑐𝛾𝜙(3.20)

The initial guesses used are𝑓0(𝜂)=𝑠+𝑒2𝜂𝜂𝑒𝜂,𝜃0(𝜂)=𝑒𝜂,𝜙0(𝜂)=𝑒𝜂.(3.21) Solving the linear part of the equation system (3.19), that is,𝑎0𝐹0+𝑎1𝐹0+𝑎2𝐹0+𝑎3𝐹0=𝑟1𝑏(𝜂),0Θ0+𝑏1Θ0+𝑏2Θ0+𝑏3𝐹0+𝑏4𝐹0=𝑟2𝑐(𝜂),0Φ0+𝑐1Φ0+𝑐2Φ0+𝑐3𝐹0+𝑐4𝐹0=𝑟3(𝜂),(3.22) subject to𝐹02(1)=𝐿𝐹02(1)=𝐿𝐹0(1)=0,Θ0(1)=Θ0(1)=0,Φ0(1)=Φ0(1)=0,(3.23) will yield the initial SHAM approximate solution. Applying the Chebyshev pseudospectral method on equations (3.22)-(3.23) yields the matrix form𝐁𝐘0=𝐑,(3.24) where𝐚𝐁=0𝒟3+𝐚1𝒟2+𝐚2𝒟+𝐚3𝐛𝐈0𝐈0𝐈3𝒟+𝐛4𝐈𝐛0𝒟2+𝐛1𝒟+𝐛2𝐜𝐈0𝐈3𝒟+𝐜4𝐈0𝐈𝐜0𝒟2+𝐜1𝒟+𝐜2𝐈𝑟,(3.25)𝐑=1𝜂0,𝑟1𝜂1,,𝑟1𝜂𝑁,𝑟2𝜂0,𝑟2𝜂1,,𝑟2𝜂𝑁,𝑟3𝜂0,𝑟3𝜂1,,𝑟3𝜂𝑁𝑇,𝐘0=𝐹0𝜉0,𝐹0𝜉1,,𝐹0𝜉𝑁,Θ0𝜉0,Θ0𝜉1,,Θ0𝜉𝑁,Φ0𝜉0,Φ0𝜉1,,Φ0𝜉𝑁𝑇,𝐚𝑖𝑎=diag𝑖𝜂0,,𝑎𝑖𝜂𝑁1,𝑎𝑖𝜂𝑁,𝐛𝑖𝑏=diag𝑖𝜂0,,𝑏𝑖𝜂𝑁1,𝑏𝑖𝜂𝑁,𝐜𝑖𝑐=diag𝑖𝜂0,,𝑐𝑖𝜂𝑁1,𝑐𝑖𝜂𝑁,𝑖=0,1,2,3,4.(3.26)The superscript 𝑇 denotes the transpose, diag is a diagonal matrix, and 𝐈 is an identity matrix of size (𝑁+1)×(𝑁+1). The boundary conditions (3.23) are implemented in matrix 𝐁 and vector 𝐑 of equation (3.24). The values of [𝑌0(𝜉1),𝑌0(𝜉2),,𝑌0(𝜉𝑁1)] are then determined from the following equation:𝐘0=𝐁1𝐑,(3.27) 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:𝐹𝐹(𝜉;𝑞)=𝑎0𝜕3𝐹𝜕𝜉3+𝑎1𝜕2𝐹𝜕𝜉2+𝑎2𝜕𝐹𝜕𝜉+𝑎3𝐹,Θ𝐹Θ(𝜉;𝑞),(𝜉;𝑞)=𝑏0𝜕2Θ𝜕𝜉2+𝑏1𝜕Θ𝜕𝜉+𝑏2Θ+𝑏3𝜕𝐹𝜕𝜉+𝑏4𝐹,Φ𝐹(𝜉;𝑞),Φ(𝜉;𝑞)=𝑐0𝜕2Φ𝜕𝜉2+𝑐1𝜕Φ𝜕𝜉+𝑐2Φ+𝑐3𝜕𝐹𝜕𝜉+𝑐4𝐹,(3.28) where 𝑞[0,1] is the embedding parameter and 𝐹(𝜉;𝑞), Θ(𝜉;𝑞), and Φ(𝜉;𝑞) are unknown functions. The zeroth-order deformation equations are given by(1𝑞)𝐹𝐹(𝜉;𝑞)𝐹0𝒩(𝜉)=𝑞𝐹𝐹(𝜉;𝑞)𝑟1,(1𝑞)ΘΘ(𝜉;𝑞)Θ0𝒩(𝜉)=𝑞Θ𝐹(𝜉;𝑞),Θ(𝜉;𝑞)𝑟2,(1𝑞)ΦΦ(𝜉;𝑞)Φ0𝒩(𝜉)=𝑞Φ𝐹(𝜉;𝑞),Φ(𝜉;𝑞)𝑟3,(3.29) where is the nonzero convergence controlling auxiliary parameter and 𝒩𝐹, 𝒩Θ, and 𝒩Φ are nonlinear operators given by𝒩𝐹𝐹(𝜉;𝑞)=𝑎0𝜕3𝐹𝜕𝜉3+𝑎1𝜕2𝐹𝜕𝜉2+𝑎2𝜕𝐹𝜕𝜉+𝑎34𝐹+𝐿2𝑚𝐹𝜕2𝐹𝜕𝜉24𝐿2𝜕𝐹𝜕𝐹𝜕𝜉,𝒩𝜕𝜉Θ𝐹Θ(𝜉;𝑞),(𝜉;𝑞)=𝑏0𝜕2Θ𝜕𝜉2+𝑏1𝜕Θ𝜕𝜉+𝑏2Θ+𝑏3𝜕𝐹𝜕𝜉+𝑏4𝐹+2𝐿Θ𝜕𝐹𝑃𝑟𝐹𝜕Θ𝜕𝜉+𝑚,𝒩𝜕𝜉Φ𝐹(𝜉;𝑞),Φ(𝜉;𝑞)=𝑐0𝜕2Φ𝜕𝜉2+𝑐1𝜕Φ𝜕𝜉+𝑐2Φ+𝑐3𝜕𝐹𝜕𝜉+𝑐4𝐹+2𝐿Φ𝜕𝐹𝑆𝑐𝐹𝜕Φ𝜕𝜉+𝑚.𝜕𝜉(3.30) The 𝑚-th order deformation equations are given by𝐹𝐹𝑚(𝜉)𝜒𝑚𝐹𝑚1(𝜉)=𝑅𝐹𝑚,ΘΘ𝑚(𝜉)𝜒𝑚Θ𝑚1(𝜉)=𝑅Θ𝑚,ΦΦ𝑚(𝜉)𝜒𝑚Φ𝑚1(𝜉)=𝑅Φ𝑚,(3.31) subject to the boundary conditions𝐹𝑚(1)=𝐹𝑚(1)=𝐹𝑚(1)=0,Θ𝑚(1)=Θ𝑚(1)=0,Φ𝑚(1)=Φ𝑚(1)=0,(3.32) where𝑅𝐹𝑚(𝜉)=𝑎0𝐹𝑚1+𝑎1𝐹𝑚1+𝑎2𝐹𝑚1+𝑎3𝐹𝑚1+4𝐿2𝑚1𝑛=0𝐹𝑛𝐹𝑚1𝑛+𝑚𝐹𝑛𝐹𝑚1𝑛𝑟1(𝜂)1𝜒𝑚,𝑅Θ𝑚(𝜉)=𝑏0Θ𝑚1+𝑏1Θ𝑚1+𝑏2Θ𝑚1+𝑏3𝐹𝑚1+𝑏4𝐹𝑚1+2𝐿𝑃𝑟𝑚1𝑛=0𝐹𝑛Θ𝑚1𝑛+𝑚Θ𝑛𝐹𝑚1𝑛𝑟2(𝜂)1𝜒𝑚,𝑅Φ𝑚(𝜉)=𝑐0Φ𝑚1+𝑐1Φ𝑚1+𝑐2Φ𝑚1+𝑐3𝐹𝑚1+𝑐4𝐹𝑚1+2𝐿𝑆𝑐𝑚1𝑛=0𝐹𝑛Φ𝑚1𝑛+𝑚Φ𝑛𝐹𝑚1𝑛𝑟3(𝜂)1𝜒𝑚,𝜒(3.33)𝑚=0,𝑚11,𝑚>1.(3.34) Applying the Chebyshev pseudospectral transformation to equations (3.31)–(3.33) gives rise to the matrix equation𝐁𝐘𝑚=𝜒𝑚+𝐁𝐘𝑚11𝜒𝑚𝐑+𝐐𝑚1,(3.35) subject to the boundary conditions𝑁𝑘=0𝒟0𝑘𝐹𝑚𝜉𝑘=0,𝑁𝑘=0𝒟𝑁𝑘𝐹𝑚𝜉𝑘=0,𝐹𝑚𝜉𝑁Θ=0,𝑚𝜉0=0,Θ𝑚𝜉𝑁=0,Φ𝑚𝜉0=0,Φ𝑚𝜉𝑁=0,(3.36) where 𝐁 and 𝐑 are as defined in (3.25) and 𝐘𝑚=𝐹𝑚𝜉0,𝐹𝑚𝜉1,,𝐹𝑚𝜉𝑁,Θ𝑚𝜉0,Θ𝑚𝜉1,,Θ𝑚𝜉𝑁,Φ𝑚𝜉0,Φ𝑚𝜉1,,Φ𝑚𝜉𝑁𝑇,𝐐𝑚1=𝑚1𝑛=04𝐿2𝒟𝐹𝑛𝒟𝐹𝑚1𝑛+4𝐿2𝑚𝐹𝑛𝒟2𝐹𝑚1𝑛2𝐿𝑃𝑟𝑚1𝑛=0𝒟𝐹𝑛Θ𝑚1𝑛+𝑚𝒟Θ𝑛𝐹𝑚1𝑛2𝐿𝑆𝑐𝑚1𝑛=0𝒟𝐹𝑛Φ𝑚1𝑛+𝑚𝒟Φ𝑛𝐹𝑚1𝑛.(3.37) Applying the boundary conditions (3.32) on the right-hand side of (3.35) yields the following recursive formula for higher-order approximations 𝑌𝑚(𝜉) for 𝑚1:𝐘𝑚=𝜒𝑚𝐁+1𝐁𝐘𝑚1+𝐁1𝐐𝑚11𝜒𝑚𝐑.(3.38)

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 𝑓(𝜂), 𝜃(𝜂), and 𝜙(𝜂) can be expanded into𝑓(𝜂)=𝐹𝑖(𝜂)+𝑖1𝑚=0𝐹𝑚(𝜂),𝜃(𝜂)=Θ𝑖(𝜂)+𝑖1𝑚=0Θ𝑚(𝜂),𝜙(𝜂)=Φ𝑖(𝜂)+𝑖1𝑚=0Φ𝑚(𝜂),𝑖=1,2,3,,(3.39) where 𝐹𝑖, Θ𝑖, and Φ𝑖 are unknown functions whose solutions are obtained using the SHAM approach at the 𝑖th iteration and 𝐹𝑚, Θ𝑚, and Φ𝑚 (𝑚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𝐹𝑖+𝑎1,𝑖1𝐹𝑖+𝑎2,𝑖1𝐹𝑖+𝑎3,𝑖1𝐹𝑖+𝑚𝐹𝑖𝐹𝑖𝐹𝑖𝐹𝑖=𝑟1,𝑖1Θ(𝜂),𝑖+𝑏1Θ𝑖+𝑏2Θ𝑖+𝑏3𝐹𝑖+𝑏4𝐹𝑖𝑃𝑟𝐹𝑖Θ𝑖+𝑚𝑃𝑟𝐹𝑖Θ𝑖=𝑟2,𝑖1Φ(𝜂),𝑖+𝑐1Φ𝑖+𝑐2Φ𝑖+𝑐3𝐹𝑖+𝑐4𝐹𝑖𝑆𝑐𝐹𝑖Φ𝑖+𝑚𝑆𝑐𝐹𝑖Φ𝑖=𝑟3,𝑖1(𝜂),(3.40) subject to the boundary conditions𝐹𝑖(0)=0,𝐹𝑖(0)=𝐹𝑖()=0,Θ𝑖(0)=0,Θ𝑖()=0,Φ𝑖(0)=0,Φ𝑖()=0.(3.41) The coefficient parameters 𝑎𝑘,𝑖1, 𝑏𝑘,𝑖1, 𝑐𝑘,𝑖1 (𝑘=0,,4), 𝑟1,𝑖1, 𝑟2,𝑖1, and 𝑟3,𝑖1 are as defined in equation (3.4). Starting from the initial guesses (3.5), the subsequent solutions 𝐹𝑖, Θ𝑖, and Φ𝑖 (𝑖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𝐹𝑖𝜕𝜂3+𝑎1,𝑖1𝜕2𝐹𝑖𝜕𝜂2+𝑎2,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑎3,𝑖1𝐹𝑖,Θ𝐹𝑖Θ(𝜂;𝑞),𝑖=𝜕(𝜂;𝑞)2Θ𝑖𝜕𝜂2+𝑏1,𝑖1𝜕Θ𝑖𝜕𝜂+𝑏2,𝑖1Θ𝑖+𝑏3,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑏4,𝑖1𝐹𝑖,Φ𝐹𝑖Φ(𝜂;𝑞),𝑖=𝜕(𝜂;𝑞)2Φ𝑖𝜕𝜂2+𝑐1,𝑖1𝜕Φ𝑖𝜕𝜂+𝑐2,𝑖1Φ𝑖+𝑐3,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑐4,𝑖1𝐹𝑖.(3.42) The zeroth-order deformation equations are given by(1𝑞)𝐹𝐹𝑖(𝜂;𝑞)𝐹𝑖,0𝒩(𝜂)=𝑞𝐹𝐹𝑖(𝜂;𝑞)𝑟1,𝑖1,(1𝑞)ΘΘ𝑖(𝜂;𝑞)Θ𝑖,0𝒩(𝜂)=𝑞Θ𝐹𝑖Θ(𝜂;𝑞),𝑖(𝜂;𝑞)𝑟2,𝑖1,(1𝑞)ΦΦ𝑖(𝜂;𝑞)Φ𝑖,0𝒩(𝜂)=𝑞Φ𝐹𝑖Φ(𝜂;𝑞),𝑖(𝜂;𝑞)𝑟3,𝑖1,(3.43)𝒩𝐹, 𝒩Θ, and 𝒩Φ are nonlinear operators given by𝒩𝐹𝐹𝑖=𝜕(𝜂;𝑞)3𝐹𝑖𝜕𝜂3+𝑎1,𝑖1𝜕2𝐹𝑖𝜕𝜂2+𝑎2,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑎3,𝑖1𝐹𝑖𝐹+𝑚𝑖𝜕2𝐹𝑖𝜕𝜂2𝜕𝐹𝑖𝜕𝐹𝜕𝜂𝑖,𝒩𝜕𝜂Θ𝐹𝑖Θ(𝜂;𝑞),𝑖=𝜕(𝜂;𝑞)2Θ𝑖𝜕𝜂2+𝑏1,𝑖1𝜕Θ𝑖𝜕𝜂+𝑏2,𝑖1Θ𝑖+𝑏3,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑏4,𝑖1𝐹𝑖Θ+𝑃𝑟𝑖𝜕𝐹𝑖𝐹𝜕𝜂+𝑚𝑖𝜕Θ𝑖,𝒩𝜕𝜂Φ𝐹𝑖Φ(𝜂;𝑞),𝑖=𝜕(𝜂;𝑞)2Φ𝑖𝜕𝜂2+𝑐1,𝑖1𝜕Φ𝑖𝜕𝜂+𝑐2,𝑖1Φ𝑖+𝑐3,𝑖1𝜕𝐹𝑖𝜕𝜂+𝑐4,𝑖1𝐹𝑖Φ+𝑆𝑐𝑖𝜕𝐹𝑖𝐹𝜕𝜂+𝑚𝑖𝜕Φ𝑖.𝜕𝜂(3.44) The 𝑚th order deformation equations are𝐹𝐹𝑖,𝑚(𝜂)𝜒𝑚𝐹𝑖,𝑚1(𝜂)=𝑅𝐹𝑖,𝑚,ΘΘ𝑖,𝑚(𝜂)𝜒𝑚Θ𝑖,𝑚1(𝜂)=𝑅Θ𝑖,𝑚,ΦΦ𝑖,𝑚(𝜂)𝜒𝑚Φ𝑖,𝑚1(𝜂)=𝑅Φ𝑖,𝑚,(3.45) subject to the boundary conditions𝐹𝑖,𝑚(0)=𝐹𝑖,𝑚(0)=𝐹𝑖,𝑚()=0,Θ𝑖,𝑚(0)=Θ𝑖,𝑚()=0,Φ𝑖,𝑚(0)=Φ𝑖,𝑚()=0,(3.46) where𝑅𝐹𝑖,𝑚(𝜂)=𝐹𝑖,𝑚1+𝑎1,𝑖1𝐹𝑖,𝑚1+𝑎2,𝑖1𝐹𝑖,𝑚1+𝑎3,𝑖1𝐹𝑖,𝑚1+𝑚1𝑛=0𝐹𝑖,𝑛𝐹𝑖,𝑚1𝑛+𝑚𝐹𝑖,𝑛𝐹𝑖,𝑚1𝑛𝑟1,𝑖1(𝜂)1𝜒𝑚,𝑅Θ𝑖,𝑚(𝜂)=Θ𝑖,𝑚1+𝑏1,𝑖1Θ𝑖,𝑚1+𝑏2,𝑖1Θ𝑖,𝑚1+𝑏3,𝑖1𝐹𝑖,𝑚1+𝑏4,𝑖1𝐹𝑖,𝑚1+𝑃𝑟𝑚1𝑛=0𝐹𝑖,𝑛Θ𝑖,𝑚1𝑛+𝑚Θ𝑖,𝑛𝐹𝑖,𝑚1𝑛𝑟2,𝑖1(𝜂)1𝜒𝑚,𝑅Φ𝑖,𝑚(𝜂)=Φ𝑖,𝑚1+𝑐1,𝑖1Φ𝑖,𝑚1+𝑐2,𝑖1Φ𝑖,𝑚1+𝑐3,𝑖1𝐹𝑖,𝑚1+𝑐4,𝑖1𝐹𝑖,𝑚1+𝑆𝑐𝑚1𝑛=0𝐹𝑖,𝑛Φ𝑖,𝑚1𝑛+𝑚Φ𝑖,𝑛𝐹𝑖,𝑚1𝑛𝑟3,𝑖1(𝜂)1𝜒𝑚.(3.47) The initial approximations 𝐹𝑖,0, Θ𝑖,0, and Φ𝑖,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𝐹𝑖,0+𝑎1,𝑖1𝐹𝑖,0+𝑎2,𝑖1𝐹𝑖,0+𝑎3,𝑖1𝐹𝑖,0=𝑟1,𝑖1,Θ𝑖,0+𝑏1,𝑖1Θ𝑖,0+𝑏2,𝑖1Θ𝑖,0+𝑏3,𝑖1𝐹𝑖,0+𝑏4,𝑖1𝐹𝑖,0=𝑟2,𝑖1,Φ𝑖,0+𝑐1,𝑖1Φ𝑖,0+𝑐2,𝑖1Φ𝑖,0+𝑐3,𝑖1𝐹𝑖,0+𝑐4,𝑖1𝐹𝑖,0=𝑟3,𝑖1,(3.48) with the boundary conditions𝐹𝑖,0(0)=𝐹𝑖,0(0)=𝐹𝑖,0()=0,Θ𝑖,0(0)=0,Θ𝑖,0()=0,Φ𝑖,0Φ(0)=0,𝑖,0()=0.(3.49) In a similar manner, we apply the spectral methods to solve for the initial approximate solutions 𝐹𝑖,0, Θ𝑖,0, and Φ𝑖,0, and the higher-order deformation equations (3.45)–(3.47) for higher order approximate solutions 𝐹𝑖,𝑚, Θ𝑖,𝑚, and Φ𝑖,𝑚 for 𝑚1. The solutions for 𝐹𝑖, Θ𝑖, and Φ𝑖 are then generated using the solutions for 𝐹𝑖,𝑚, Θ𝑖,𝑚, and Φ𝑖,𝑚 as follows:𝐹𝑖=𝐹𝑖,0+𝐹𝑖,1+𝐹𝑖,2+𝐹𝑖,3++𝐹𝑖,𝑚,Θ𝑖=Θ𝑖,0+Θ𝑖,1+Θ𝑖,2+Θ𝑖,3++Θ𝑖,𝑚,Φ𝑖=Φ𝑖,0+Φ𝑖,1+Φ𝑖,2+Φ𝑖,3++Φ𝑖,𝑚.(3.50) The [𝑖,𝑚] approximate solutions for 𝑓(𝜂), 𝜃(𝜂), and 𝜙(𝜂) are then obtained by substituting 𝐹𝑖, Θ𝑖, and Φ𝑖 from (3.50) into (3.39), where 𝑖 is the 𝑖th iteration of the higher-order deformation equation and 𝑚 is the 𝑚th 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 13.

tab1
Table 1: 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.
tab2
Table 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.
tab3
Table 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.

Tables 13 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.

Table 1 shows an increase in the surface shear stress 𝑓(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 (𝑠=3), suggesting an increase in temperature and concentration gradients with increasing porosity.

Figure 1 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 𝑠 and the Hartmann number 𝑀 on the velocity profiles 𝑓(𝜂).

fig1
Figure 1: On the comparison between the 2nd-order ISHAM solution (figures) and the bvp4c numerical solution (solid line) for f(η) and θ(η) at different values of 𝜆 when 𝑀=1, 𝑚=1, 𝑃𝑟=3, 𝜆1=2, 𝑠=1, 𝐿=30, and 𝑁=150.

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 (𝑚=1) and show that increasing the suction parameter 𝑠 and the Hartmann number 𝑀 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.

Acknowledgments

The authors wish to acknowledge financial support from the University of KwaZulu-Natal, University of Venda, and the National Research Foundation (NRF).

References

  1. N. Bachok, A. Ishak, and I. Pop, “Unsteady three-dimensional boundary layer flow due to a permeable shrinking sheet,” Applied Mathematics and Mechanics, vol. 31, no. 11, pp. 1421–1428, 2010. View at Publisher · View at Google Scholar
  2. T.-G. Fang, J. Zhang, and S. S. Yao, “Slip magnetohydrodynamic viscous flow over a permeable shrinking sheet,” Chinese Physics Letters, vol. 27, no. 12, article 124702, 2010. View at Publisher · View at Google Scholar
  3. T. Fang, S. Yao, J. Zhang, and A. Aziz, “Viscous flow over a shrinking sheet with a second order slip flow model,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 7, pp. 1831–1842, 2010. View at Publisher · View at Google Scholar
  4. T. Fang and J. Zhang, “Closed-form exact solutions of MHD viscous flow over a shrinking sheet,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 7, pp. 2853–2857, 2009. View at Publisher · View at Google Scholar · View at Scopus
  5. Muhaimin, R. Kandasamy, and A. B. Khamis, “Effects of heat and mass transfer on nonlinear MHD boundary layer flow over a shrinking sheet in the presence of suction,” Applied Mathematics and Mechanics, vol. 29, no. 10, pp. 1309–1317, 2008. View at Publisher · View at Google Scholar · View at Scopus
  6. N. F. M. Noor, S. Awang Kechil, and I. Hashim, “Simple non-perturbative solution for MHD viscous flow due to a shrinking sheet,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 2, pp. 144–148, 2010. View at Publisher · View at Google Scholar
  7. N. F. Mohd and I. Hashim, “MHD flow and heat transfer adjacent to a permeable shrinking sheet embedded in a porous medium,” Sains Malaysiana, vol. 38, no. 4, pp. 559–565, 2009. View at Google Scholar · View at Scopus
  8. S. Nadeem and A. Hussain, “MHD flow of a viscous fluid on a nonlinear porous shrinking sheet with homotopy analysis method,” Applied Mathematics and Mechanics, vol. 30, no. 12, pp. 1569–1578, 2009. View at Publisher · View at Google Scholar
  9. Muhaimin, R. Kandasamy, I. Hashim, and A. B. Khamis, “On the effect of chemical reaction, heat and mass transfer on nonlinear MHD boundary layer past a porous shrinking sheet with suction,” Theoretical and Applied Mechanics, vol. 36, no. 2, pp. 101–117, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. S. Abbasbandy, “A numerical solution of Blasius equation by Adomian's decomposition method and comparison with homotopy perturbation method,” Chaos, Solitons and Fractals, vol. 31, no. 1, pp. 257–260, 2007. View at Publisher · View at Google Scholar
  11. N. S. Elgazery, “Numerical solution for the Falkner-Skan equation,” Chaos, Solitons and Fractals, vol. 35, no. 4, pp. 738–746, 2008. View at Publisher · View at Google Scholar · View at Scopus
  12. B. L. Kuo, “Heat transfer analysis for the Falkner-Skan wedge flow by the differential transformation method,” International Journal of Heat and Mass Transfer, vol. 48, no. 23-24, pp. 5036–5046, 2005. View at Publisher · View at Google Scholar · View at Scopus
  13. A. Wazwaz, “The variational iteration method for solving two forms of Blasius equation on a half-infinite domain,” Applied Mathematics and Computation, vol. 188, no. 1, pp. 485–491, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  14. S. J. Liao, Beyond Perturbation, Introduction to Homotopy Analysis Method, Chapman & Hall, Boca Raton, Fla, USA, 2003.
  15. B. Yao and J. Chen, “A new analytical solution branch for the Blasius equation with a shrinking sheet,” Applied Mathematics and Computation, vol. 215, no. 3, pp. 1146–1153, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. B. Yao and J. Chen, “Series solution to the Falkner-Skan equation with stretching boundary,” Applied Mathematics and Computation, vol. 208, no. 1, pp. 156–164, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  17. A. A. Joneidi, G. Domairry, M. Babaelahi, and M. Mozaffari, “Analytical treatment on magnetohydrodynamic (MHD) flow and heat transfer due to a stretching hollow cylinder,” International Journal for Numerical Methods in Fluids, vol. 63, no. 5, pp. 548–563, 2010. View at Google Scholar
  18. A. R. Noiey, N. Haghparast, M. Miansari, and D. D. Ganji, “Application of homotopy perturbation method to the MHD pipe flow of a fourth grade fluid,” Journal of Physics, vol. 96, no. 1, 2008. View at Publisher · View at Google Scholar
  19. M. Jalaal and D. D. Ganji, “An analytical study on motion of a sphere rolling down an inclined plane submerged in a Newtonian fluid,” Powder Technology, vol. 198, no. 1, pp. 82–92, 2010. View at Publisher · View at Google Scholar · View at Scopus
  20. M. Jalaal, D. D. Ganji, and G. Ahmadi, “An analytical study on settling of non-spherical particles,” Asia-Pacific Journal of Chemical Engineering, 2010. View at Publisher · View at Google Scholar
  21. M. Jalaal, D. D. Ganji, and G. Ahmadi, “Analytical investigation on acceleration motion of a vertically falling spherical particle in incompressible Newtonian media,” Advanced Powder Technology, vol. 21, no. 3, pp. 298–304, 2010. View at Publisher · View at Google Scholar · View at Scopus
  22. M. Jalaal and D. D. Ganji, “On unsteady rolling motion of spheres in inclined tubes filled with incompressible Newtonian fluids,” Advanced Powder Technology, vol. 22, pp. 58–67, 2011. View at Google Scholar
  23. S. M. Moghimi, D. D. Ganji, H. Bararnia, M. Hosseini, and M. Jalaal, “Homotopy perturbation method for nonlinear MHD Jeffery-Hamel problem,” Computers & Mathematics with Applications, vol. 61, no. 8, pp. 2213–2216, 2011. View at Publisher · View at Google Scholar
  24. S. S. Motsa and S. Shateyi, “A new approach for the solution of three-dimensional magnetohydrodynamic rotating flow over a shrinking sheet,” Mathematical Problems in Engineering, vol. 2010, Article ID 586340, 15 pages, 2010. View at Publisher · View at Google Scholar
  25. S. S. Motsa, P. Sibanda, and S. Shateyi, “A new spectral-homotopy analysis method for solving a nonlinear second order BVP,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 9, pp. 2293–2302, 2010. View at Publisher · View at Google Scholar
  26. S. S. Motsa, P. Sibanda, F. G. Awad, and S. Shateyi, “A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem,” Computers & Fluids, vol. 39, no. 7, pp. 1219–1225, 2010. View at Publisher · View at Google Scholar
  27. Z. G. Makukula, P. Sibanda, and S. Motsa, “On a new solution for the viscoelastic squeezing flow between two parallel plates,” Journal of Advanced Research in Applied Mathematics, vol. 2, no. 4, pp. 31–38, 2010. View at Google Scholar
  28. Z. G. Makukula, P. Sibanda, and S. S. Motsa, “A note on the solution of the von Kármán equations using series and Chebyshev spectral methods,” Boundary Value Problems, vol. 2010, Article ID 471793, 17 pages, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  29. Zodwa G. Makukula, Precious Sibanda, and Sandile S. Motsa, “A novel numerical technique for two-dimensional laminar flow between two moving porous walls,” Mathematical Problems in Engineering, vol. 2010, Article ID 528956, 15 pages, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  30. Z. G. Makukula, P. Sibanda, and S. Motsa, “On a new quasi-linearisation method for heat transfer in a visco-elastic fluid between parallel plates,” in New Aspects of Fluid Mechanics, Heat Transfer and Environment, N. Mastorakis, V. Mladenov, and Z. Bojkovic, Eds., 8th iasme/wseas international conference on fluid mechanics & Aerodynamics (FMA '10), pp. 178–186, Tapei, Taiwan, August 2010.
  31. Z. G. Makukula, P. Sibanda, and S. S. Motsa, “On new solutions for heat transfer in a visco-elastic fluid between parallel plates,” International Journal of Mathematical Models and Methods in Applied Sciences, vol. 4, no. 4, pp. 221–230, 2010. View at Google Scholar · View at Scopus
  32. M. Sajid and T. Hayat, “The application of homotopy analysis method for MHD viscous flow due to a shrinking sheet,” Chaos, Solitons and Fractals, vol. 39, no. 3, pp. 1317–1323, 2009. View at Publisher · View at Google Scholar · View at Scopus
  33. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, Germany, 1988.
  34. L. N. Trefethen, Spectral Methods in MATLAB, vol. 10, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2000.
  35. H. S. Takhar, A. J. Chamkha, and G. Nath, “Natural convection on a thin vertical cylinder moving in a high-porosity ambient medium,” International Journal of Engineering Science, vol. 41, no. 16, pp. 1935–1950, 2003. View at Publisher · View at Google Scholar · View at Scopus