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

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

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

The applicable boundary conditions are𝑢=−𝑎𝑥,𝑣𝑢=⟶−𝑎0(𝑚,−𝐶1⟶)𝑦,𝐶𝑤_∞=−𝑊,𝑇=𝑇_𝑤,𝐶=𝐶_𝑤,at𝑦=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𝑓^−𝜃𝑀²+𝜆𝑃𝑟𝑓′−𝑓′²+𝑚𝑓𝑓^=0,(2.4)_𝜙−𝑃𝑟𝑓^𝜃+𝑚𝑃𝑟𝑓𝜃^=0,(2.5)_−𝑆𝑐𝑓^𝜙+𝑚𝑆𝑐𝑓𝜙^−𝑆𝑐𝛾𝜙=0,(2.6) subject to𝑓(0)=𝑠,𝑓^(0)=−𝜙1(,0)=𝑓1,^(∞)𝜙=(∞0,)=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 10⁻⁶.

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_𝑚₌₀𝐹_𝑚(𝜂),𝜃(𝜂)=𝜃_𝑖(𝜂)+^𝑖−1_𝑚₌₀𝜙(𝜂)=𝜙_𝑖Θ_𝑚(𝜂),(𝜂)+^𝑖−1_𝑚₌₀Φ_𝑚(𝜂),𝑖=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_𝑚₌₀𝐹_𝑚,𝑎_2,𝑖−1=−2^𝑖−1_𝑚₌₀𝐹^𝑚+𝑀²+𝜆𝑃𝑟,𝑎_3,𝑖−1=𝑚^𝑖−1_𝑚₌₀𝑏_1,𝑖−1𝐹^𝑚,=𝑚𝑃𝑟^𝑖−1_𝑚₌₀𝐹_𝑚,𝑏_2,𝑖−1=−𝑃𝑟^𝑖−1_𝑚₌₀𝐹^𝑚,𝑏_3,𝑖−1=−𝑃𝑟^𝑖−1_𝑚₌₀𝑏_4,𝑖−1Θ_𝑚,=𝑚𝑃𝑟^𝑖−1_𝑚₌₀Θ^𝑚,𝑐_1,𝑖−1=𝑚𝑆𝑐^𝑖−1_𝑚₌₀𝐹_𝑚,𝑐_2,𝑖−1=−𝑆𝑐^𝑖−1_𝑚₌₀𝑐𝐹^𝑚−𝑆𝑐𝛾,_3,𝑖−1=−𝑆𝑐^𝑖−1_𝑚₌₀Φ_𝑚,𝑐_4,𝑖−1=𝑚𝑆𝑐^𝑖−1_𝑚₌₀𝑟_1,𝑖−1Φ^𝑚,=−^𝑖−1_𝑚₌₀𝐹^{𝑚}−^𝑖−1_𝑚₌₀𝐹^{𝑚𝑖−1}_𝑚₌₀𝐹^𝑚+𝑚^𝑖−1_𝑚₌₀𝐹^{𝑚𝑖−1}_𝑚₌₀𝐹_𝑚−𝑀²+𝜆𝑃𝑟^𝑖−1_𝑚₌₀𝐹^𝑚𝑟_2,𝑖−1,=−^𝑖−1_𝑚₌₀Θ^𝑚−𝑃𝑟^𝑖−1_𝑚₌₀𝐹^{𝑚𝑖−1}_𝑚₌₀Θ_𝑚+𝑚𝑃𝑟^𝑖−1_𝑚₌₀𝐹_{𝑚𝑖−1𝑚}₌₀Θ^𝑚𝑟_3,𝑖−1,=−^𝑖−1_𝑚₌₀Φ^𝑚−𝑆𝑐^𝑖−1_𝑚₌₀𝐹^{𝑚𝑖−1}_𝑚₌₀Φ_𝑚+𝑚𝑆𝑐^𝑖−1_𝑚₌₀𝐹_{𝑚𝑖−1𝑚}₌₀Φ^𝑚−𝑆𝑐𝛾^𝑖−1_𝑚₌₀Φ.(3.4)

Starting from the initial approximations𝐹₀(𝜂)=𝑠+𝑒^−2𝜂−𝑒^−𝜂,Θ₀(𝜂)=𝑒^−𝜂,Φ₀(𝜂)=𝑒^−𝜂,(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𝑓(𝜂)≈^𝑀_𝑚₌₀𝐹_𝑚(𝜂),𝜃(𝜂)≈^𝑀_𝑚₌₀Θ_𝑚(𝜂),𝜙(𝜂)≈^𝑀_𝑚₌₀Φ_𝑚(𝜂).(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𝑘cos⁻¹(𝜉).(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=⎛⎜⎜⎜⎜⎝𝐃³+𝐚𝐛₁𝐃²+𝐚₂𝐃+𝐚₃𝐈0𝐈0𝐈₃𝐜𝐃+𝐛₄𝐈𝐃²+𝐛₁𝐃+𝐛₂𝐈0𝐈₃𝐗_𝑖𝐃+𝐜₄𝐈0𝐈𝐃²+𝐜₁𝐃+𝐜₂𝐈⎞⎟⎟⎟⎟⎠,=𝐹Φ_𝑚𝑚𝜉₀,𝐹_𝑚𝜉₁,…,𝐹_𝑚𝜉_𝑁,Θ_𝑚𝜉₀,Θ_𝑚𝜉₁,…,Θ_𝑚𝜉_𝑁,𝐑_𝑖−1𝜉₀,Φ_𝑚𝜉₁,…,Φ_𝑚𝜉_𝑁^𝑇,=𝑟𝑟_{1,𝑖−13,𝑖−1}𝜂₀,𝑟_1,𝑖−1𝜂₁,…,𝑟_1,𝑖−1𝜂_𝑁,𝑟_2,𝑖−1𝜂₀,𝑟_2,𝑖−1𝜂₁,…,𝑟_2,𝑖−1𝜂_𝑁,𝜂₀,𝑟_3,𝑖−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𝑓(𝜂)=𝐹(𝜉)+𝑓₀(𝜂),𝜃(𝜂)=Θ(𝜉)+𝜃₀(𝜂),𝜙(𝜂)=Φ(𝜉)+𝜙₀(𝜂),(3.17) where 𝑓₀(𝜂), and 𝜃₀(𝜂), 𝜙₀(𝜂) 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𝐿²𝑓𝐹^(𝜉)+𝑓^0(𝜂),_(𝜂)=8𝐿³𝜃𝐹^(𝜉)+𝑓^0(𝜂),_(𝜂)=2𝐿Θ^(𝜉)+𝜃^0(𝜂),𝜃^(𝜂)=4𝐿²𝜙Θ^(𝜉)+𝜃^0(𝜂),_(𝜂)=2𝐿Φ^(𝜉)+𝜙^0(𝜂),𝜙^(𝜂)=4𝐿²Φ^(𝜉)+𝜙^0(𝜂).(3.18) Substituting (3.1) and (3.17)-(3.18) in the governing equations and boundary conditions gives𝑎₀𝐹^+𝑎₁𝐹^+𝑎₂𝐹^+𝑎₃𝐹+4𝐿²𝑚𝐹^𝐹−4𝐿²𝑏𝐹^𝐹^=𝑟₁(𝜂),₀Θ^+𝑏₁Θ^+𝑏₂Θ+𝑏₃𝐹^+𝑏₄𝐹−2𝐿𝑃𝑟𝐹^Θ+2𝐿𝑚𝑃𝑟𝐹Θ^𝑐=𝑟₂(𝜂),₀Φ^+𝑐₁Φ^+𝑐₂Φ+𝑐₃𝐹^+𝑐₄𝐹−2𝐿𝑆𝑐𝐹^Φ+2𝐿𝑚𝑆𝑐𝐹Φ^=𝑟₃(𝜂),(3.19) where prime now denotes derivative with respect to 𝜉 and 𝑎₀=8𝐿³,𝑎₁=4𝐿²𝑚𝑓₀,𝑎₂=−2𝐿𝑀²+𝜆𝑃𝑟+4𝐿𝑓^0,𝑎³𝑟₁=𝑚𝑓^0,𝑏₀(𝜂)=−𝑓^0+𝑚𝑓₀𝑓^0−𝑓^0𝑓^0−𝑀²+𝜆𝑃𝑟𝑓^0,=4𝐿²,𝑏₁=2𝐿𝑚𝑃𝑟𝑓⁰,𝑏₂=−𝑃𝑟𝑓^0,𝑏₃=−2𝐿𝑃𝑟𝜃⁰𝑟₂,𝑏₄=𝑚𝑃𝑟𝜃^0,𝑐₀(𝜂)=−𝜃^−𝑃𝑟𝑓^0𝜃₀+𝑚𝑃𝑟𝑓₀𝜃₀,=4𝐿²,𝑐₁=2𝐿𝑚𝑆𝑐𝑓⁰,𝑐₂=−𝑆𝑐𝑓^0+𝛾,𝑐₃=−2𝐿𝑆𝑐𝜙⁰𝑟₃,𝑐₄=𝑚𝑆𝑐𝜙^0,(𝜂)=−𝜙^−𝑆𝑐𝑓^0𝜙₀+𝑚𝑆𝑐𝑓₀𝜙^0−𝑆𝑐𝛾𝜙.(3.20)

The initial guesses used are𝑓₀(𝜂)=𝑠+𝑒^−2𝜂−𝜂𝑒^−𝜂,𝜃₀(𝜂)=𝑒^−𝜂,𝜙₀(𝜂)=𝑒^−𝜂.(3.21) Solving the linear part of the equation system (3.19), that is,𝑏𝑎₀𝐹^0+𝑎₁𝐹^0+𝑎₂𝐹^0+𝑎₃𝐹₀=𝑟₁(𝜂),₀𝑐Θ^0+𝑏₁Θ^0+𝑏₂Θ₀+𝑏₃𝐹^0+𝑏₄𝐹₀=𝑟₂(𝜂),₀Φ^0+𝑐₁Φ^0+𝑐₂Φ₀+𝑐₃𝐹^0+𝑐₄𝐹₀=𝑟₃(𝜂),(3.22) subject to𝐹₀(−1)=2𝐿𝐹^0(−1)=2𝐿𝐹^0(1)=0,Θ₀(−1)=Θ₀(1)=0,Φ₀(−1)=Φ₀(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𝐁𝐘₀=𝐑,(3.24) where𝐁=⎛⎜⎜⎜⎜⎝𝐚⁰𝒟³+𝐛𝐚₁𝒟²+𝐚₂𝒟+𝐚₃𝐈0𝐈0𝐈₃𝐜𝒟+𝐛₄𝐈𝐛₀𝒟²+𝐛₁𝒟+𝐛₂𝐈0𝐈₃𝐑=𝑟₁𝒟+𝐜₄𝐈0𝐈𝐜₀𝒟²+𝐜₁𝒟+𝐜₂𝐈⎞⎟⎟⎟⎟⎠,(3.25)𝐘₀𝜂₀,𝑟₁𝜂₁,…,𝑟₁𝜂_𝑁,𝑟₂𝜂₀,𝑟₂𝜂₁,…,𝑟₂𝜂_𝑁,𝑟₃𝜂₀,𝑟₃𝜂₁,…,𝑟₃𝜂_𝑁^𝑇,𝐚_𝑖=𝐹₀𝜉₀,𝐹₀𝜉₁,…,𝐹₀𝜉_𝑁,Θ₀𝜉₀,Θ₀𝜉₁,…,Θ₀𝜉_𝑁,Φ₀𝜉₀,Φ₀𝜉₁,…,Φ₀𝜉_𝑁^𝑇,𝐜=diag𝑎_𝑖𝜂₀,…,𝑎_𝑖𝜂_𝑁−1,𝑎_𝑖𝜂_𝑁,𝐛_𝑖=diag𝑏_𝑖𝜂₀,…,𝑏_𝑖𝜂_𝑁−1,𝑏_𝑖𝜂_𝑁_𝑖,=diag𝑐_𝑖𝜂₀,…,𝑐_𝑖𝜂_𝑁−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 [𝑌₀(𝜉₁),𝑌₀(𝜉₂),…,𝑌₀(𝜉_𝑁−1)] are then determined from the following equation:𝐘₀=𝐁⁻¹𝐑,(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:ℒ_𝐹𝐹(𝜉;𝑞)=𝑎⁰𝜕³𝐹𝜕𝜉³+𝑎₁𝜕²𝐹𝜕𝜉²+𝑎₂𝜕𝐹𝜕𝜉+𝑎³ℒ𝐹,_Θ𝐹(𝜉;𝑞),Θ(𝜉;𝑞)=𝑏⁰𝜕²Θ𝜕𝜉²+𝑏₁𝜕Θ𝜕𝜉+𝑏²Θ+𝑏₃𝜕𝐹𝜕𝜉+𝑏⁴ℒ𝐹,_Φ𝐹(𝜉;𝑞),Φ(𝜉;𝑞)=𝑐⁰𝜕²Φ𝜕𝜉²+𝑐₁𝜕Φ𝜕𝜉+𝑐²Φ+𝑐₃𝜕𝐹𝜕𝜉+𝑐⁴𝐹,(3.28) where 𝑞∈[0,1] is the embedding parameter and 𝐹(𝜉;𝑞), Θ(𝜉;𝑞), and Φ(𝜉;𝑞) are unknown functions. The zeroth-order deformation equations are given by((11−−𝑞𝑞))ℒℒ_𝐹Θ𝐹(𝜉;𝑞)−𝐹⁰(𝜉)=𝑞ℏ𝒩_𝐹𝐹(𝜉;𝑞)−𝑟¹,(1−𝑞)ℒ_ΦΘ(𝜉;𝑞)−Θ⁰(𝜉)=𝑞ℏ𝒩_Θ𝐹(𝜉;𝑞),Θ(𝜉;𝑞)−𝑟²,Φ(𝜉;𝑞)−Φ⁰(𝜉)=𝑞ℏ𝒩_Φ𝐹(𝜉;𝑞),Φ(𝜉;𝑞)−𝑟³,(3.29) where ℏ is the nonzero convergence controlling auxiliary parameter and 𝒩_𝐹, 𝒩_Θ, and 𝒩_Φ are nonlinear operators given by𝒩_𝐹𝐹(𝜉;𝑞)=𝑎⁰𝜕³𝐹𝜕𝜉³+𝑎₁𝜕²𝐹𝜕𝜉²+𝑎₂𝜕𝐹𝜕𝜉+𝑎³𝐹+4𝐿²𝑚𝐹𝜕²𝐹𝜕𝜉²−4𝐿²𝜕𝐹𝒩𝜕𝜉𝜕𝜕𝐹𝜉_Θ,𝐹(𝜉;𝑞),Θ(𝜉;𝑞)=𝑏⁰𝜕²Θ𝜕𝜉²+𝑏₁𝜕Θ𝜕𝜉+𝑏²Θ+𝑏₃𝜕𝐹𝜕𝜉+𝑏⁴+2𝐹𝒩𝐿𝑃𝑟−Θ𝜕𝜕𝐹𝜉+𝑚𝐹𝜕𝜕Θ𝜉_Φ,𝐹(𝜉;𝑞),Φ(𝜉;𝑞)=𝑐⁰𝜕²Φ𝜕𝜉²+𝑐₁𝜕Φ𝜕𝜉+𝑐²Φ+𝑐₃𝜕𝐹𝜕𝜉+𝑐⁴+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𝑅^𝐹𝑚(𝜉)=𝑎₀𝐹^{𝑚−1}+𝑎₁𝐹^{𝑚−1}+𝑎₂𝐹^𝑚−1+𝑎₃𝐹_𝑚−1+4𝐿^2𝑚−1_𝑛=0𝑅_Θ𝑚−𝐹^𝑛𝐹^{𝑚−1−𝑛}+𝑚𝐹_𝑛𝐹^{𝑚−1−𝑛}−𝑟₁(𝜂)1−𝜒_𝑚,(𝜉)=𝑏₀Θ^{𝑚−1}+𝑏₁Θ^𝑚−1+𝑏₂Θ_𝑚−1+𝑏₃𝐹^𝑚−1+𝑏₄𝐹_𝑚−1+2𝐿𝑃𝑟^𝑚−1_𝑛=0−𝐹^𝑛𝑅_Φ𝑚Θ_{𝑚−1−𝑛}+𝑚Θ^𝑛𝐹_{𝑚−1−𝑛}−𝑟₂(𝜂)1−𝜒_𝑚,(𝜉)=𝑐₀Φ^{𝑚−1}+𝑐₁Φ^𝑚−1+𝑐₂Φ_𝑚−1+𝑐₃𝐹^𝑚−1+𝑐₄𝐹_𝑚−1+2𝐿𝑆𝑐^𝑚−1_𝑛=0−𝐹^𝑛𝜒Φ_{𝑚−1−𝑛}+𝑚Φ^𝑛𝐹_{𝑚−1−𝑛}−𝑟₃(𝜂)1−𝜒_𝑚,(3.33)_𝑚=01,,𝑚𝑚≤>11.(3.34) Applying the Chebyshev pseudospectral transformation to equations (3.31)–(3.33) gives rise to the matrix equation𝐁𝐘_𝑚=𝜒_𝑚+ℏ𝐁𝐘_𝑚−1−ℏ1−𝜒_𝑚𝐑+ℏ𝐐_𝑚−1,(3.35) subject to the boundary conditions^𝑁_𝑘=0𝒟_0𝑘𝐹_𝑚𝜉_𝑘=0,^𝑁_𝑘=0Θ𝒟_𝑁𝑘𝐹_𝑚𝜉_𝑘=0,𝐹_𝑚𝜉_𝑁=0,_𝑚𝜉₀=0,Θ_𝑚𝜉_𝑁=0,Φ_𝑚𝜉₀=0,Φ_𝑚𝜉_𝑁=0,(3.36) where 𝐁 and 𝐑 are as defined in (3.25) and 𝐘_𝑚=𝐹Φ_𝑚𝑚𝜉₀,𝐹_𝑚𝜉₁,…,𝐹_𝑚𝜉_𝑁,Θ_𝑚𝜉₀,Θ_𝑚𝜉₁,…,Θ_𝑚𝜉_𝑁,𝐐_𝑚−1𝜉₀,Φ_𝑚𝜉₁,…,Φ_𝑚𝜉_𝑁^𝑇,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝^𝑚−1_𝑛=0−𝐿4²𝒟𝐹_𝑛𝒟𝐹_{𝑚−1−𝑛}+𝐿4₂2𝑚𝐹_𝑛𝒟²𝐹_{𝑚−1−𝑛}𝐿𝑃𝑟^𝑚−1_𝑛=0−𝒟𝐹_𝑛Θ_{𝑚−1−𝑛}2+𝑚𝒟Θ_𝑛𝐹_{𝑚−1−𝑛}𝐿𝑆𝑐^𝑚−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−𝜒_𝑚𝐑.(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_𝑚₌₀𝐹_𝑚(𝜂),𝜃(𝜂)=Θ_𝑖(𝜂)+^𝑖−1_𝑚₌₀𝜙(𝜂)=Φ_𝑖Θ_𝑚(𝜂),(𝜂)+^𝑖−1_𝑚₌₀Φ_𝑚(𝜂),𝑖=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(𝜂),_𝑖Φ+𝑏₁Θ^𝑖+𝑏₂Θ_𝑖+𝑏₃𝐹^𝑖+𝑏₄𝐹_𝑖−𝑃𝑟𝐹^𝑖Θ_𝑖+𝑚𝑃𝑟𝐹_𝑖Θ^𝑖=𝑟_2,𝑖−1(𝜂),_𝑖+𝑐₁Φ^𝑖+𝑐₂Φ_𝑖+𝑐₃𝐹^𝑖+𝑐₄𝐹_𝑖−𝑆𝑐𝐹^𝑖Φ_𝑖+𝑚𝑆𝑐𝐹_𝑖Φ^𝑖=𝑟_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:ℒ_𝐹𝐹^𝑖(𝜂;𝑞)=𝜕³𝐹_𝑖𝜕𝜂³+𝑎_1,𝑖−1𝜕²𝐹_𝑖𝜕𝜂²+𝑎_2,𝑖−1𝜕𝐹_𝑖𝜕𝜂+𝑎^3,𝑖−1ℒ_Θ𝐹_𝑖,𝐹^𝑖(𝜂;𝑞),Θ_𝑖(𝜂;𝑞)=𝜕²Θ_𝑖𝜕𝜂²+𝑏_1,𝑖−1𝜕Θ_𝑖𝜕𝜂+𝑏^2,𝑖−1Θ_𝑖+𝑏_3,𝑖−1𝜕𝐹_𝑖𝜕𝜂+𝑏^4,𝑖−1ℒ_Φ𝐹_𝑖,𝐹^𝑖(𝜂;𝑞),Φ_𝑖(𝜂;𝑞)=𝜕²Φ_𝑖𝜕𝜂²+𝑐_1,𝑖−1𝜕Φ_𝑖𝜕𝜂+𝑐^2,𝑖−1Φ_𝑖+𝑐_3,𝑖−1𝜕𝐹_𝑖𝜕𝜂+𝑐^4,𝑖−1𝐹_𝑖.(3.42) The zeroth-order deformation equations are given by((11−−𝑞𝑞))ℒℒ_𝐹Θ𝐹^𝑖(𝜂;𝑞)−𝐹_𝑖,0(𝜂)=𝑞ℏ𝒩_𝐹𝐹^𝑖(𝜂;𝑞)−𝑟_1,𝑖−1,(1−𝑞)ℒ_ΦΘ^𝑖(𝜂;𝑞)−Θ_𝑖,0(𝜂)=𝑞ℏ𝒩_Θ𝐹^𝑖(𝜂;𝑞),Θ_𝑖(𝜂;𝑞)−𝑟_2,𝑖−1,Φ^𝑖(𝜂;𝑞)−Φ_𝑖,0(𝜂)=𝑞ℏ𝒩_Φ𝐹^𝑖(𝜂;𝑞),Φ_𝑖(𝜂;𝑞)−𝑟_3,𝑖−1,(3.43)𝒩_𝐹, 𝒩_Θ, and 𝒩_Φ are nonlinear operators given by𝒩_𝐹𝐹^𝑖(𝜂;𝑞)=𝜕³𝐹_𝑖𝜕𝜂³+𝑎_1,𝑖−1𝜕²𝐹_𝑖𝜕𝜂²+𝑎_2,𝑖−1𝜕𝐹_𝑖𝜕𝜂+𝑎^3,𝑖−1𝐹_𝑖+𝑚𝐹_𝑖𝜕²𝐹_𝑖𝜕𝜂²−𝜕𝐹^𝑖𝒩𝜕𝜂𝜕𝜕𝐹𝜂^𝑖_Θ,𝐹^𝑖(𝜂;𝑞),Θ_𝑖(𝜂;𝑞)=𝜕²Θ_𝑖𝜕𝜂²+𝑏_1,𝑖−1𝜕Θ_𝑖𝜕𝜂+𝑏^2,𝑖−1Θ_𝑖+𝑏_3,𝑖−1𝜕𝐹_𝑖𝜕𝜂+𝑏^4,𝑖−1𝐹_𝑖+𝑃𝑟−Θ_𝑖𝜕𝐹_𝑖𝜕𝜂+𝑚𝐹^𝑖𝜕Θ_𝑖𝒩𝜕𝜂_Φ,𝐹^𝑖(𝜂;𝑞),Φ_𝑖(𝜂;𝑞)=𝜕²Φ_𝑖𝜕𝜂²+𝑐_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𝐹^{𝑖,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,+𝑐_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 1–3.

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.

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.

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 1–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.

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 𝑓^(𝜂).

(a)

(b)

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).