# On the Multidomain Bivariate Spectral Local Linearisation Method for Solving Systems of Nonsimilar Boundary Layer Partial Differential Equations

Authors: Vusi Mpendulo Magagula
Journal: International Journal of Mathematics and Mathematical Sciences
 (2019)

---

## Abstract

In this work, a novel approach for solving systems of nonsimilar boundary layer equations over a large time domain is presented. The method is a multidomain bivariate spectral local linearisation method (MD-BSLLM), Legendre-Gauss-Lobatto grid points, a local linearisation technique, and the spectral collocation method to approximate functions defined as bivariate Lagrange interpolation. The method is developed for a general system ofn nonlinear partial differential equations. The use of the MD-BSLLM is demonstrated by solving a system of nonlinear partial differential equations that describe a class of nonsimilar boundary layer equations. Numerical experiments are conducted to show applicability and accuracy of the method. Grid independence tests establish the accuracy, convergence, and validity of the method. The solution for the limiting case is used to validate the accuracy of the MD-BSLLM. The proposed numerical method performs better than some existing numerical methods for solving a class of nonsimilar boundary layer equations over large time domains since it converges faster and uses few grid points to achieve accurate results. The proposed method uses minimal computation time and its accuracy does not deteriorate with an increase in time.


---

## Body

## 1. Introduction

Solutions of nonlinear nonsimilar boundary layer partial differential equations have received considerable attention of research scientist in the past few years. Series solution methods have been used by several researchers. These include asymptotic methods [1, 2], extended series solutions [2], and homotopy analysis method (HAM) [3–8]. To achieve accurate results when using series based methods, a large number of terms are required. A complex nonsimilar boundary layer differential equation requires many terms and thus solving the problem becomes tedious. Simplification of a complex model may introduce errors and may lead to inaccurate solutions.Numerical methods have been used to solve nonsimilar boundary layer equations. Numerical methods like implicit finite difference methods have been used to solve nonsimilar boundary layer equations by Hossain [1, 9–11]. Hussain [2] and Yih [12] used the Keller-Box method. Finite difference methods require many grid points to achieve accurate results. The accuracy of the results deteriorates fast as the time domain increases for parabolic nonlinear partial differential equations. Recently, Magagula et al. [13] showed that pseudospectral methods can achieve accurate results with few grid points compared to Keller-Box method. Results were obtained in a fraction of a second.Pseudospectral methods have been used to solve boundary layer ordinary differential equations [14–19], and boundary layer partial differential equations [20–22] have been solved using spectral local linearisation method (SLLM) and the bivariate spectral local linearisation method (BSLLM), respectively. Other pseudospectral methods have been used to solve differential equations [23–25]. The bivariate spectral local linearisation method was developed and generalized by Motsa [16] for solving systems of boundary layer ordinary differential equations. This method uses quasilinearisation to linearise nonlinear coupled system of differential equations. The result is a system of linear decoupled differential equations which can be solved independently. The SLLM was extended by Motsa [21] to a system of nonlinear coupled boundary layer partial differential equations, where spectral collocation was used in both space and time. The BSLLM method is suitable for nonlinear partial differential equations with small time domain. Increasing the time domain reduces the accuracy of the method. In some cases, the computed numerical solutions fail to converge even when the time step is decreased. Thus, accurate solutions of systems of differential equations are usually determined over short time intervals. The realization that accurate solutions can be obtained in small time intervals have led to a number of numerical methods that uses the domain-decomposition technique. The domain-decomposition approach together with the spectral collocation technique has been used to solve equations arising from various scientific fields. Funaro [26] and Harald [27] solved elliptic differential equations, Kopriva [28] solved hyperbolic systems, Peyret [29] solved stiff problems in fluid mechanics, and recently, it has been used to solve chaotic system of ordinary differential equations by [30–33] and Lane-Emden type of equations [34]. These methods split large intervals into smaller subintervals so that the equations are solved over a sequence of nonoverlapping subintervals of the domain. Continuity conditions are used to advance the solution across the nonoverlapping subintervals.In this paper, we present a multidomain numerical approach for finding solutions of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval. The multidomain bivariate spectral local linearisation method (MD-BSLLM) is an extension of the BSLLM. The proposed approach is developed for solving system ofn coupled nonsimilar boundary layer partial differential equations. The system of equations is solved iteratively using a Jacobi like relaxation approach. The domain is divided into smaller nonoverlapping subintervals on which the Chebyshev spectral collocation method is used to solve the equations. A continuity condition is used to advance the solution across the subintervals. The algorithm is easy to develop and it yields very accurate results using only a few discretization nodes. The accuracy of the method is validated against the series solution for limiting cases. The main aim of the study is to explore the applicability of the multidomain bivariate spectral local linearisation collocation method to systems of nonlinear coupled nonsimilar boundary layer partial differential equations over large time intervals. The results confirm that the method is suitable for solving all types of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval.The article is organized as follows: in Section2, the multidomain bivariate spectral local linearisation method is described for a general system of n partial differential equations. In Section 3, numerical experiments and the series solutions are presented in Section 4. In Section 5, the results are discussed and we conclude in Section 6.


## 2. Multidomain Bivariate Spectral Local Linearisation Method

In this section, the multidomain bivariate spectral local linearisation method (MD-BSLLM) for approximating solutions of system of nonlinear partial differential equations that model nonsimilar boundary layer equations is introduced and explained in detail. The multidomain bivariate spectral local linearisation method is developed for a general system ofn nonlinear partial differential equations. Consider a system of n nonlinear partial differential equations of the form,(1)Γ1H1,H2,…,Hn=0,(2)Γ2H1,H2,…,Hn=0,(3)⋮(4)ΓnH1,H2,…,Hn=0,where(5)H1=f1,∂f1∂x,∂2f1∂x2,…,∂pf1∂xp,∂f1∂τ,∂∂τ∂f1∂x,(6)H2=f2,∂f2∂x,∂2f2∂x2,…,∂pf2∂xp,∂f2∂τ,∂∂τ∂f2∂x,(7)⋮(8)Hn=fn,∂fn∂x,∂2fn∂x2,…,∂pfn∂xp,∂fn∂τ,∂∂τ∂fn∂x.The order of differentiation is denoted by p; the required solutions by fk(x,τ) for k=1,2,…,n and Γk for k=1,2,…,n are the nonlinear operators containing all the spatial derivatives and time derivatives of fk(x,τ). The Legendre-Gauss-Lobatto grid points and the corresponding differentiation matrices are defined in the interval [-1,1]. Therefore, the time interval τ∈[0,T] and space region x∈[a,b] are transformed to [-1,1] using linear transformations ζ=T(τ+1)/2 and η=(1/2)(b-a)x+(1/2)(b+a), respectively. Thus fk=fk(η,ζ). In order to apply the multidomain Legendre-Gauss-Lobatto bivariate spectral local linearisation method on the system of equations under consideration, the interval ζ is first divided into z nonoverlapping subintervals Ωl=[ζl-1,ζl] for l=1,2,…,z. Each interval Ωl is further divided into s not necessarily equal divisions. The system of nonlinear partial differential equations under consideration is solved in each subinterval [ζl-1,ζl]. The solutions in each subinterval [ζl-1,ζl] are denoted by fk(l)(η,ζ), for k=1,2,…,n. In the first interval [ζ0,ζ1], the solutions fk(1)(η,ζ), for k=1,2,…,n are obtained subject to the “initial" conditions fk(1)(η,0), for k=1,2,…,n, respectively. For each l>1, at each subinterval [ζl-1,ζl], the continuity conditions(9)fklη,ζl-1=fkl-1η,ζl-1,fork=1,2,…,nandl=1,2,…,z,are used to implement the MD-BSQLM over the interval [ζl-1,ζl]. This procedure is repeated to generate a sequence of solutions fk(l)(η,ζ), for l=1,2,…,z and k=1,2,…,n. It is assumed that, at each subinterval [ζl-1,ζl], the solution can be approximated by a bivariate Lagrange interpolation polynomial of the form(10)fklη,ζ≈∑i=0Nη∑j=0Nζfklηi,ζjLiηLjζ,for k=1,2,…,n and l=1,2,…,z. The bivariate Lagrange interpolation polynomial interpolates fk(l)(η,ζ) at selected points (ηi,ζj) in both the η and ζ directions, for i=0,1,2,…,Nη and j=0,1,2,…,Nζ. The grid points in time are given by the zeros of the polynomial(11)pζ=1-ζ2PNζ′ζ,where PNζ′(ζ) is the first derivative of the Legendre polynomial given by(12)PNζ′ζ=NζNζ+12Nζ+1PNζ-1t-PNζ+1ζ1-ζ2.The zeros of PNζ′(ζ) can be computed using any root finding method. In this work, we use Newton’s method. Hence given ζj0 and for m≥0, (1≤j≤Nζ-1), we have(13)ζjm+1=ζjm-PNζ′ζjmPNζ′′ζjm.We note that(14)1-ζ2PNζ′′ζ-2ζPNζ′ζ+ζζ+1PNζζ=0,and thus(15)PNζ′′ζ=2ζPNζ′ζ-ζζ+1PNζζ1-ζ2.Substituting (15) into (13), we get(16)ζjm+1=ζjm-1-ζjm2PNζ′ζjm2ζjmPNζ′ζjm-ζjmζjm+1PNζζjm.Similarly, the grid points in space are given by the zeros of the polynomial(17)pη=1-η2PNη′η,and hence the zeros of PNη′(η) can be obtained by using the modified Newton method(18)ηim+1=ηim-1-ηim2PNη′ηim2ηimPNη′ηim-ηimηim+1PNηηim.The zeros of (11) and (17) are called Legendre-Gauss-Lobatto points [35, 36]. The function Li(η) is the characteristic Lagrange cardinal polynomial with the Legendre-Gauss-Lobatto points [35, 36](19)Liη=∏i=0i≠kNηη-ηkηi-ηk,where(20)Liηk=δik=0ifi≠k1ifi=k.An alternative convenient expression for (19) is given by(21)Liη=1NηNη+1PNηηiη2-1PNη′ηη-ηi.The nonlinear operators Γk, for k=1,2,3,…,n, are first linearised using the quasilinearisation technique as defined by Bellman and Kalaba [37]. The quasilinearisation method is a technique based on the Taylor series expansion of Γk about some previous iteration. It is assumed that the differences between its previous and current solution and all its derivatives are small. Applying the quasilinearisation method yields the following:(22)ΓkH1,H2,…,Hn≈H1,r+1-H1,r,H2,r+1-H2,r,…,Hn,r+1-Hn,r·∇ΓkH1,r,H2,r,…,Hn,r+ΓkH1,r,H2,r,…,Hn,r, where r and r+1 denote previous and current iterations, respectively, and ∇ is a vector of the partial derivatives which is defined as(23)∇=∇f1,∇f2,…,∇fn. We define(24)∇f1=∂∂f1,∂∂f1′,∂∂f1′′,…,∂∂f1p,∂∂∂f1/∂ζ,∂∂∂f1′/∂ζ,(25)∇f2=∂∂f2,∂∂f2′,∂∂f2′′,…,∂∂f2p,∂∂∂f2/∂ζ,∂∂∂f2′/∂ζ,(26)⋮(27)∇fn=∂∂fn,∂∂fn′,∂∂fn′′,…,∂∂fnp,∂∂∂fn/∂ζ,∂∂∂fn′/∂ζ,where the prime denotes differentiation with respect to η. The linearised equation (22) can be expressed in a compact form as (28)H1,r+1·∇f1Γ1H1,r,H2,r,…,Hn,r=H1,r·∇f1Γ1H1,r,H2,r,…,Hn,r-Γ1H1,r,H2,r,…,Hn,r(29)H2,r+1·∇f2Γ2H1,r+1,H2,r,…,Hn,r=H2,r·∇f2Γ2H1,r+1,H2,r,…,Hn,r-Γ2H1,r+1,H2,r,…,Hn,r(30)⋮(31)Hn,r+1·∇fnΓnH1,r+1,H2,r+1,…,Hn-1,r+1,Hn,r=Hn,r·∇fnΓnH1,r+1,H2,r+1,…,Hn-1,r+1,Hn,r-ΓnH1,r+1,H2,r+1,…,Hn-1,r+1,Hn,rEquations (28)-(31) form a system of n decoupled linear partial differential equations. They are solved iteratively for f1(l)(η,t),f2(l)(η,t),…,fn(l)(η,t). Equations (28)-(31) can further be expressed as follows:(32)∑s=0pαs,rl,1η,ζf1,r+1s,l+βrl,1η,ζ∂f1,r+10,l∂ζ+γrl,1η,ζ∂f1,r+11,l∂ζ=R1lη,ζ,(33)∑s=0pαs,rl,2η,ζf2,r+1s,l+βrl,2η,ζ∂f2,r+10,l∂ζ+γrl,2η,ζ∂f2,r+11,l∂ζ=R2lη,ζ,(34)⋮(35)∑s=0pαs,rl,nη,ζfn,r+1s,l+βrl,nη,ζ∂fn,r+10,l∂ζ+γrl,nη,ζ∂fn,r+11,l∂ζ=Rnlη,ζ,where αs,r(l,k)(η,ζ), βr(l,k)(η,ζ), and γr(l,k)(η,ζ) are variable coefficients of fk,r+1(s,l), ∂fk,r+1(0,l)/∂ζ, and ∂fk,r+1(1,l)/∂ζ, respectively, in the lth subinterval of the multidomain approach, for k=1,2,…,n and s=0,1,2,…,p. These variable coefficients correspond to the kth equation, for k=1,2,…,n. The constant p denotes the order of differentiation. Thus, we have(36)αs,rl,kη,ζ=∂Γk∂fk,rs,l,(37)βrl,kη,ζ=∂Γk∂∂fk,r0,l/∂ζ,(38)γrl,kη,ζ=∂Γk∂∂fk,r1,l/∂ζ.In general, the kth right hand side is given by(39)Rklη,ζ=∑s=0pαs,rl,kη,ζfk,rs,l+βrl,kη,ζ∂fk,r0,l∂ζ+γrl,kη,ζ∂fk,r1,l∂ζ)-Γk.Equations (32), (33), and (35) are evaluated at the Legendre-Gauss-Lobatto grid points ζj(j=0,1,…,Nζ) and ηi(i=0,1,…,Nη). The values of the time derivatives are computed at the Legendre-Gauss-Lobatto points (ηi,ζj), as (for j=0,1,2,…,Nζ)(40)∂fnl∂ζηi,ζj=∑s=0Nη∑k=0Nζfnlηs,ζkLsηidLkζjdζ(41)=∑k=0Nζfnlηi,ζkdjk=∑k=0Nζdjkfnlηi,ζkwhere djk=dLk(ζj)/dζ is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size (Nζ+1)×(Nζ+1), given by [36](42)djk=-NζNζ+14ifk=j=0PNζtjPNζζkζj-ζkifk≠j,0≤j,k≤Nζ0if1≤j=k≤Nζ-1NζNζ+14ifk=j=NζThe values of the space derivatives at the Legendre-Gauss-Lobatto points (ηi,ζj) (for i=0,1,2,…,Nη) are similarly computed as(43)∂fnl∂ηηi,ζj=∑s=0Nη∑k=0Nζfnlηs,ζkdLsηidηLkζj(44)=∑s=0Nηfnlηs,ζjDis=∑s=0NηDisfnlηs,ζj,where Dis=dLs(ηi)/dη is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size (Nη+1)×(Nη+1) as defined in [36]. The entries of the space first derivative matrix with respect to the Legendre-Gauss-Lobatto points are exactly the same as those given by (42). Higher nth order derivatives are defined as(45)∂pfnl∂ηpηi,tj=∑s=0NηDispfnlηs,ζj=DpFn,jl,i=0,1,2,…,Nη,where the vector Fn,j(l) is defined as(46)Fn,jl=fnlη0,ζjl,fnlη1,ζjl,…,fnlηNη,ζjlT.and the superscript T denotes matrix transpose. Substituting (44) and (45) into (32), (33), and (35) yields(47)A1,1lF1,il+βrl,1∑j=0Nζdi,jF1,jl+γrl,1∑j=0Nζdi,jDF1,jl=R1,il,(48)A2,2lF2,il+βrl,2∑j=0Nζdi,jF2,jl+γrl,2∑j=0Nζdi,jDF2,jl=R2,il,(49)⋮(50)An,nlFn,il+βrl,n∑j=0Nζdi,jFn,jl+γrl,n∑j=0Nζdi,jDFn,jl=Rn,il,where(51)A1,1l=∑s=0pαs,rl,1Ds,A2,2l=∑s=0pαs,rl,2Ds,…,An,nl=∑s=0pαs,rl,nDs.The diagonal matrices of the corresponding variable coefficients are given by(52)αs,rl,k=αs,rl,kη0,ζjαs,rl,kη1,ζj⋱αs,rl,kηNη,ζj,(53)βrl,k=βrl,kη0,ζjβrl,kη1,ζj⋱βrl,kηNη,ζj,(54)γrl,k=γrl,kη0,ζjγrl,kη1,ζj⋱γrl,kηNη,ζj.Imposing boundary conditions for i=0,1,…,Nζ-1, (47), (48), and (50) can be expressed as the following Nζ(Nη+1)×Nζ(Nη+1) matrix system(55)B0,0l,kB0,1l,k⋯B0,Nζ-1l,kB1,0l,kB1,1l,k⋯B1,Nζ-1l,k⋮⋮⋱⋮BNζ-1,0l,kBNζ-1,1l,k⋯BNζ-1,Nζ-1l,kFk,0lFk,1l⋮Fk,Nζ-1l=Rk,0lRk,1l⋮Rk,Nζ-1l,where(56)Bi,il,k=∑s=0pαs,rl,kDs+βrl,kdi,iI+γrl,kdi,iD,fork=1,2,…,n,l=1,2,…,zwheni=j,Bi,jl,k=βrl,kdi,jI+γrl,kdi,jD,fork=1,2,…,n,l=1,2,…,zwheni≠j.The vector Rk,i(l) is defined as (57)Rk,il=Rk,il-βrl,kdi,NζI+γrl,kdi,NζDFk,Nζl for i=0,1,…,Nζ-1 and k=1,2,…,n. The vector Fk,Nζ(l) corresponds to the initial boundary condition which is always prescribed.


## 3. Numerical Experiments

In this section, we present systems of equation for various fluid flow models considered in this study. The MD-BSLLM algorithm is applied to each system of equations as a numerical experiment.


### 3.1. Steady Free Convection Flow Past a Nonisothermal Vertical Porous Cone

In this subsection, the problem of steady two-dimensional laminar free convection flow past a nonisothermal vertical porous cone with variable surface temperature is also considered for numerical experimentation. The governing nonsimilarity system of partial differential equations is expressed in dimensionless form as (see Hossain et. al. [1] for derivation)(58)f′′′+n+74ff′′-n+12f′2+θ+ξf′′=1-n4ξf′∂f′∂ξ-f′′∂f∂ξ,(59)1Prθ′′+n+74fθ′-nf′θ+ξθ′=1-n4ξf′∂θ∂ξ-θ′∂f∂ξ,where Pr=ν/α is the Prandtl number, ξ is the dimensionless suction parameter, and θ(η,ξ) and f(η,ξ) are the dimensionless temperature and stream functions, respectively. The appropriate corresponding boundary conditions are(60)fξ,0=0,f′ξ,0=0,θξ,0=1,f′ξ,∞=0,θξ,∞=0.The skin friction coefficient Cfx and the Nusselt number Nux describe the shear-stress and heat flux rate at the surface, respectively, and are defined by [1] as(61)CfxGrx1/4=f′′ξ,0,NuxGrx1/4=-θ′ξ,0.In this case, the order of differentiation p=3 and the number of system of equations n=2. In each subinterval Ωl, we have the equation(62)A1,1lF1,il+βrl,1∑j=0Nζdi,jF1,jl+γrl,1∑j=0Nζdi,jDF1,jl=R1,il,(63)A2,2lF2,il+βrl,2∑j=0Nζdi,jF2,jl+γrl,2∑j=0Nζdi,jDF2,jl=R2,il,where(64)A1,1l=∑s=03αs,rl,1Ds=α3,rl,1D3+α2,rl,1D2+α1,rl,1D1+α0,rl,1,(65)A2,2l=∑s=03αs,rl,2Ds=α3,rl,2D3+α2,rl,2D2+α1,rl,2D1+α0,rl,2.We note that the equations are in terms of the functions f(η,ζ) and θ(η,ζ). In order to apply our method, we let f(η,ζ)=f1(η,ζ) and θ(η,ζ)=f2(η,ζ). Therefore, we have (66)Γ1=f1,r3,l+n+74f1,r0,lf1,r2,l-n+12f1,r1,l2+f2,r0,l+ζf1,r2,l-1-n4ζf1,r1,l∂f1,r1,l∂ζ-f1,r2,l∂f1,r0,l∂ζ,(67)Γ2=1Prf2,r2,l+n+74f1,r0,lf2,r1,l-nf1,r1,lf2,r0,l+ζf2,r1,l-1-n4ζf1,r1,l∂f2,r0,l∂ζ-f2,r1,l∂f1,r0,l∂ζ.The variable coefficients in (64) are given by(68)α3,rl,1η,ζ=∂Γ1∂f1,r3,l=1,α2,rl,1η,ζ=∂Γ1∂f1,r2,l=14n+7f1,r1,l+ζ,α1,rl,1η,ζ=∂Γ1∂f1,r1,l=-n+1f1,r1,l,α0,rl,1η,ζ=∂Γ1∂f1,r0,l=14n+7f1,r2,l,βrl,1η,ζ=∂Γ1∂∂f1,r0,l/∂ζ=141-nζf1,r2,l,γrl,1η,ζ=∂Γ1∂∂f1,r1,l/∂ζ=-141-nζf1,r1,l.and the variable coefficients for (65) are given by(69)γrl,2η,ζ=∂Γ2∂∂f2,r1,l/∂ζ=0,α3,rl,2η,ζ=∂Γ2∂f2,r3,l=0,α2,rl,2η,ζ=∂Γ2∂f2,r2,l=1Pr,α1,rl,2η,ζ=∂Γ2∂f2,r1,l=14n+7f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,2η,ζ=∂Γ2∂f2,r0,l=-nf1,r1,l,βrl,2η,ζ=∂Γ2∂∂f2,r0,l/∂ζ=-141-nζf1,r1,l.Equations (62) and (63) can be expressed respectively as the following Nζ(Nη+1)×Nζ(Nη+1) matrix system:(70)B0,0l,1B0,1l,1⋯B0,Nζ-1l,1B1,0l,1B1,1l,1⋯B1,Nζ-1l,1⋮⋮⋱⋮BNζ-1,0l,1BNζ-1,1l,1⋯BNζ-1,Nζ-1l,1F1,0lF1,1l⋮F1,Nζ-1l=R1,0lR1,1l⋮R1,Nζ-1l,(71)B0,0l,2B0,1l,2⋯B0,Nζ-1l,2B1,0l,2B1,1l,2⋯B1,Nζ-1l,2⋮⋮⋱⋮BNζ-1,0l,2BNζ-1,1l,2⋯BNζ-1,Nζ-1l,2F2,0lF2,1l⋮F2,Nζ-1l=R2,0lR2,1l⋮R2,Nζ-1l,where(72)Bi,il,1=∑s=03αs,rl,1Ds+βrl,1di,iI+γrl,1di,iD,forl=1,2,…,zwheni=j,Bi,jl,1=βrl,1di,jI+γrl,1di,jD,for,l=1,2,…,zwheni≠j,(73)Bi,il,2=∑s=03αs,rl,2Ds+βrl,2di,iI+γrl,2di,iD,forl=1,2,…,zwheni=j,Bi,jl,2=βrl,2di,jI+γrl,2di,jD,for,l=1,2,…,zwheni≠j.The vectors R1,i(l) and R2,i(l) are defined as (74)R1,il=R1,il-βrl,1di,NζI+γrl,1di,NζDF1,Nζl,R2,il=R2,il-βrl,2di,NζI+γrl,2di,NζDF2,Nζl.for i=0,1,…,Nζ-1.


### 3.2. Steady Free Convective Flow of a Viscous Incompressible Fluid along a Permeable Vertical Flat Plate

A two-dimensional steady free convective flow of a viscous incompressible fluid along a permeable vertical flat plate in the presence of soluble species is being considered in this subsection. The problem considered here is that of a natural convection boundary layer flow, influenced by the combined buoyancy forces from mass and thermal diffusion from a permeable vertical flat surface with nonuniform surface temperature and nonuniform surface species concentration, but with a uniform rate of suction of fluid through the permeable surface. The transformed boundary layer equations are solved numerically near to and far from the leading edge, using the multidomain bivariate spectral collocation method. The equations were solved by Hussain [2] and are(75)f′′′+14n+3ff′′-12n+1f′2+ζf′′+1-wg+wh=141-nζf′∂f′∂ζ-f′′∂f∂ζ(76)1Prg′′+14n+3fg′+ζg′=141-nζf′∂g∂ζ-g′∂f∂ζ(77)1Sch′′+14n+3fh′+ζh′=141-nζf′∂h∂ζ-h′∂f∂ζThe appropriate boundary conditions are given by(78)fζ,0=0,f′ζ,0=0,gζ,0=hζ,0=1,f′ζ,∞=gζ,∞=hζ,∞=0.In this case, the order of differentiation p=3 and the number of system of equations n=3. Similarly, in each subinterval Ωl, we have(79)A1,1lF1,il+βrl,1∑j=0Nζdi,jF1,jl+γrl,1∑j=0Nζdi,jDF1,jl=R1,il,(80)A2,2lF2,il+βrl,2∑j=0Nζdi,jF2,jl+γrl,2∑j=0Nζdi,jDF2,jl=R2,il,(81)A3,3lF3,il+βrl,3∑j=0Nζdi,jF3,jl+γrl,3∑j=0Nζdi,jDF3,jl=R3,il,where(82)A1,1l=∑s=03αs,rl,1Ds=α3,rl,1D3+α2,rl,1D2+α1,rl,1D1+α0,rl,1,(83)A2,2l=∑s=03αs,rl,2Ds=α3,rl,2D3+α2,rl,2D2+α1,rl,2D1+α0,rl,2,(84)A3,3l=∑s=03αs,rl,3Ds=α3,rl,3D3+α2,rl,3D2+α1,rl,3D1+α0,rl,3.The system of equations is in terms of the functions, f(η,ζ), g(η,ζ), and h(η,ζ) and in order to apply our method, we set f(η,ζ)=f1(η,ζ), g(η,ζ)=f2(η,ζ), and h(η,ζ)=f3(η,ζ). This leads to (85)Γ1=f1,r3,l+14n+3f1,r0,lf1,r2,l-12n+1f1,r1,l2+ζf1,r2,l+1-wf2,r0,l+wf3,r0,l-141-nζf1,r1,l∂f1,r1,l∂ζ-f1,r2,l∂f1,r0,l∂ζ(86)Γ2=1Prf2,r2,l+14n+3f1,r0,lf2,r1,l+ζf2,r1,l-141-nζf1,r1,l∂f2,r0,l∂ζ-f2,r1,l∂f1,r0,l∂ζ(87)Γ3=1Scf3,r2,l+14n+3f1,r0,lf3,r1,l+ζf3,r1,l-141-nζf1,r1,l∂f3,r0,l∂ζ-f3,r1,l∂f1,r0,l∂ζThe variable coefficients for (82) are given by (88)α3,rl,1η,ζ=∂Γ1∂f1,r3,l=1,α2,rl,1η,ζ=∂Γ1∂f1,r2,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α1,rl,1η,ζ=∂Γ1∂f1,r1,l=-n+1f1,r1,l-141-nζ∂f1,r0,l∂ζ,α0,rl,1η,ζ=∂Γ1∂f1,r0,l=14n+3f1,r2,l,βrl,1η,ζ=∂Γ1∂∂f1,r0,l/∂ζ=141-nζf1,r2,l,γrl,1η,ζ=∂Γ1∂∂f1,r1,l/∂ζ=-141-nζf1,r1,l,for (83) are given by(89)γrl,2η,ζ=∂Γ2∂∂f2,r1,l/∂ζ=0,α3,rl,2η,ζ=∂Γ2∂f2,r3,l=0,α2,rl,2η,ζ=∂Γ2∂f2,r2,l=1Pr,α1,rl,2η,ζ=∂Γ2∂f2,r1,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,2η,ζ=∂Γ2∂f2,r0,l=0,βrl,2η,ζ=∂Γ2∂∂f2,r0,l/∂ζ=-141-nζf1,r1,l,and for (84) are given by (90)γrl,3η,ζ=∂Γ3∂∂f3,r1,l/∂ζ=0,α3,rl,3η,ζ=∂Γ3∂f3,r3,l=0,α2,rl,3η,ζ=∂Γ3∂f3,r2,l=1Sc,α1,rl,3η,ζ=∂Γ3∂f3,r1,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,3η,ζ=∂Γ3∂f3,r0,l=0,βrl,3η,ζ=∂Γ3∂∂f3,r0,l/∂ζ=-141-nζf1,r1,l.Equations (79), (80), and (81) can be expressed, respectively, as the following Nζ(Nη+1)×Nζ(Nη+1) matrix system:(91)B0,0l,1B0,1l,1⋯B0,Nζ-1l,1B1,0l,1B1,1l,1⋯B1,Nζ-1l,1⋮⋮⋱⋮BNζ-1,0l,1BNζ-1,1l,1⋯BNζ-1,Nζ-1l,1F1,0lF1,1l⋮F1,Nζ-1l=R1,0lR1,1l⋮R1,Nζ-1l,(92)B0,0l,2B0,1l,2⋯B0,Nζ-1l,2B1,0l,2B1,1l,2⋯B1,Nζ-1l,2⋮⋮⋱⋮BNζ-1,0l,2BNζ-1,1l,2⋯BNζ-1,Nζ-1l,2F2,0lF2,1l⋮F2,Nζ-1l=R2,0lR2,1l⋮R2,Nζ-1l,(93)B0,0l,3B0,1l,3⋯B0,Nζ-1l,3B1,0l,3B1,1l,3⋯B1,Nζ-1l,3⋮⋮⋱⋮BNζ-1,0l,3BNζ-1,1l,3⋯BNζ-1,Nζ-1l,3F3,0lF3,1l⋮F3,Nζ-1l=R3,0lR3,1l⋮R3,Nζ-1l,where(94)Bi,il,1=∑s=03αs,rl,1Ds+βrl,1di,iI+γrl,1di,iD,forl=1,2,…,zwheni=j,Bi,jl,1=βrl,1di,jI+γrl,1di,jD,for,l=1,2,…,zwheni≠j,(95)Bi,il,2=∑s=03αs,rl,2Ds+βrl,2di,iI+γrl,2di,iD,forl=1,2,…,zwheni=j,Bi,jl,2=βrl,2di,jI+γrl,2di,jD,for,l=1,2,…,zwheni≠j,(96)Bi,il,3=∑s=03αs,rl,3Ds+βrl,3di,iI+γrl,3di,iD,forl=1,2,…,zwheni=j,Bi,jl,3=βrl,3di,jI+γrl,3di,jD,for,l=1,2,…,zwheni≠j.The vectors R1,i(l), R2,i(l) and R3,i(l) are defined as (97)R1,il=R1,il-βrl,1di,NζI+γrl,1di,NζDF1,Nζl,R2,il=R2,il-βrl,2di,NζI+γrl,2di,NζDF2,Nζl,R3,il=R3,il-βrl,3di,NζI+γrl,3di,NζDF3,Nζlfor i=0,1,…,Nζ-1.


## 3.1. Steady Free Convection Flow Past a Nonisothermal Vertical Porous Cone

In this subsection, the problem of steady two-dimensional laminar free convection flow past a nonisothermal vertical porous cone with variable surface temperature is also considered for numerical experimentation. The governing nonsimilarity system of partial differential equations is expressed in dimensionless form as (see Hossain et. al. [1] for derivation)(58)f′′′+n+74ff′′-n+12f′2+θ+ξf′′=1-n4ξf′∂f′∂ξ-f′′∂f∂ξ,(59)1Prθ′′+n+74fθ′-nf′θ+ξθ′=1-n4ξf′∂θ∂ξ-θ′∂f∂ξ,where Pr=ν/α is the Prandtl number, ξ is the dimensionless suction parameter, and θ(η,ξ) and f(η,ξ) are the dimensionless temperature and stream functions, respectively. The appropriate corresponding boundary conditions are(60)fξ,0=0,f′ξ,0=0,θξ,0=1,f′ξ,∞=0,θξ,∞=0.The skin friction coefficient Cfx and the Nusselt number Nux describe the shear-stress and heat flux rate at the surface, respectively, and are defined by [1] as(61)CfxGrx1/4=f′′ξ,0,NuxGrx1/4=-θ′ξ,0.In this case, the order of differentiation p=3 and the number of system of equations n=2. In each subinterval Ωl, we have the equation(62)A1,1lF1,il+βrl,1∑j=0Nζdi,jF1,jl+γrl,1∑j=0Nζdi,jDF1,jl=R1,il,(63)A2,2lF2,il+βrl,2∑j=0Nζdi,jF2,jl+γrl,2∑j=0Nζdi,jDF2,jl=R2,il,where(64)A1,1l=∑s=03αs,rl,1Ds=α3,rl,1D3+α2,rl,1D2+α1,rl,1D1+α0,rl,1,(65)A2,2l=∑s=03αs,rl,2Ds=α3,rl,2D3+α2,rl,2D2+α1,rl,2D1+α0,rl,2.We note that the equations are in terms of the functions f(η,ζ) and θ(η,ζ). In order to apply our method, we let f(η,ζ)=f1(η,ζ) and θ(η,ζ)=f2(η,ζ). Therefore, we have (66)Γ1=f1,r3,l+n+74f1,r0,lf1,r2,l-n+12f1,r1,l2+f2,r0,l+ζf1,r2,l-1-n4ζf1,r1,l∂f1,r1,l∂ζ-f1,r2,l∂f1,r0,l∂ζ,(67)Γ2=1Prf2,r2,l+n+74f1,r0,lf2,r1,l-nf1,r1,lf2,r0,l+ζf2,r1,l-1-n4ζf1,r1,l∂f2,r0,l∂ζ-f2,r1,l∂f1,r0,l∂ζ.The variable coefficients in (64) are given by(68)α3,rl,1η,ζ=∂Γ1∂f1,r3,l=1,α2,rl,1η,ζ=∂Γ1∂f1,r2,l=14n+7f1,r1,l+ζ,α1,rl,1η,ζ=∂Γ1∂f1,r1,l=-n+1f1,r1,l,α0,rl,1η,ζ=∂Γ1∂f1,r0,l=14n+7f1,r2,l,βrl,1η,ζ=∂Γ1∂∂f1,r0,l/∂ζ=141-nζf1,r2,l,γrl,1η,ζ=∂Γ1∂∂f1,r1,l/∂ζ=-141-nζf1,r1,l.and the variable coefficients for (65) are given by(69)γrl,2η,ζ=∂Γ2∂∂f2,r1,l/∂ζ=0,α3,rl,2η,ζ=∂Γ2∂f2,r3,l=0,α2,rl,2η,ζ=∂Γ2∂f2,r2,l=1Pr,α1,rl,2η,ζ=∂Γ2∂f2,r1,l=14n+7f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,2η,ζ=∂Γ2∂f2,r0,l=-nf1,r1,l,βrl,2η,ζ=∂Γ2∂∂f2,r0,l/∂ζ=-141-nζf1,r1,l.Equations (62) and (63) can be expressed respectively as the following Nζ(Nη+1)×Nζ(Nη+1) matrix system:(70)B0,0l,1B0,1l,1⋯B0,Nζ-1l,1B1,0l,1B1,1l,1⋯B1,Nζ-1l,1⋮⋮⋱⋮BNζ-1,0l,1BNζ-1,1l,1⋯BNζ-1,Nζ-1l,1F1,0lF1,1l⋮F1,Nζ-1l=R1,0lR1,1l⋮R1,Nζ-1l,(71)B0,0l,2B0,1l,2⋯B0,Nζ-1l,2B1,0l,2B1,1l,2⋯B1,Nζ-1l,2⋮⋮⋱⋮BNζ-1,0l,2BNζ-1,1l,2⋯BNζ-1,Nζ-1l,2F2,0lF2,1l⋮F2,Nζ-1l=R2,0lR2,1l⋮R2,Nζ-1l,where(72)Bi,il,1=∑s=03αs,rl,1Ds+βrl,1di,iI+γrl,1di,iD,forl=1,2,…,zwheni=j,Bi,jl,1=βrl,1di,jI+γrl,1di,jD,for,l=1,2,…,zwheni≠j,(73)Bi,il,2=∑s=03αs,rl,2Ds+βrl,2di,iI+γrl,2di,iD,forl=1,2,…,zwheni=j,Bi,jl,2=βrl,2di,jI+γrl,2di,jD,for,l=1,2,…,zwheni≠j.The vectors R1,i(l) and R2,i(l) are defined as (74)R1,il=R1,il-βrl,1di,NζI+γrl,1di,NζDF1,Nζl,R2,il=R2,il-βrl,2di,NζI+γrl,2di,NζDF2,Nζl.for i=0,1,…,Nζ-1.


## 3.2. Steady Free Convective Flow of a Viscous Incompressible Fluid along a Permeable Vertical Flat Plate

A two-dimensional steady free convective flow of a viscous incompressible fluid along a permeable vertical flat plate in the presence of soluble species is being considered in this subsection. The problem considered here is that of a natural convection boundary layer flow, influenced by the combined buoyancy forces from mass and thermal diffusion from a permeable vertical flat surface with nonuniform surface temperature and nonuniform surface species concentration, but with a uniform rate of suction of fluid through the permeable surface. The transformed boundary layer equations are solved numerically near to and far from the leading edge, using the multidomain bivariate spectral collocation method. The equations were solved by Hussain [2] and are(75)f′′′+14n+3ff′′-12n+1f′2+ζf′′+1-wg+wh=141-nζf′∂f′∂ζ-f′′∂f∂ζ(76)1Prg′′+14n+3fg′+ζg′=141-nζf′∂g∂ζ-g′∂f∂ζ(77)1Sch′′+14n+3fh′+ζh′=141-nζf′∂h∂ζ-h′∂f∂ζThe appropriate boundary conditions are given by(78)fζ,0=0,f′ζ,0=0,gζ,0=hζ,0=1,f′ζ,∞=gζ,∞=hζ,∞=0.In this case, the order of differentiation p=3 and the number of system of equations n=3. Similarly, in each subinterval Ωl, we have(79)A1,1lF1,il+βrl,1∑j=0Nζdi,jF1,jl+γrl,1∑j=0Nζdi,jDF1,jl=R1,il,(80)A2,2lF2,il+βrl,2∑j=0Nζdi,jF2,jl+γrl,2∑j=0Nζdi,jDF2,jl=R2,il,(81)A3,3lF3,il+βrl,3∑j=0Nζdi,jF3,jl+γrl,3∑j=0Nζdi,jDF3,jl=R3,il,where(82)A1,1l=∑s=03αs,rl,1Ds=α3,rl,1D3+α2,rl,1D2+α1,rl,1D1+α0,rl,1,(83)A2,2l=∑s=03αs,rl,2Ds=α3,rl,2D3+α2,rl,2D2+α1,rl,2D1+α0,rl,2,(84)A3,3l=∑s=03αs,rl,3Ds=α3,rl,3D3+α2,rl,3D2+α1,rl,3D1+α0,rl,3.The system of equations is in terms of the functions, f(η,ζ), g(η,ζ), and h(η,ζ) and in order to apply our method, we set f(η,ζ)=f1(η,ζ), g(η,ζ)=f2(η,ζ), and h(η,ζ)=f3(η,ζ). This leads to (85)Γ1=f1,r3,l+14n+3f1,r0,lf1,r2,l-12n+1f1,r1,l2+ζf1,r2,l+1-wf2,r0,l+wf3,r0,l-141-nζf1,r1,l∂f1,r1,l∂ζ-f1,r2,l∂f1,r0,l∂ζ(86)Γ2=1Prf2,r2,l+14n+3f1,r0,lf2,r1,l+ζf2,r1,l-141-nζf1,r1,l∂f2,r0,l∂ζ-f2,r1,l∂f1,r0,l∂ζ(87)Γ3=1Scf3,r2,l+14n+3f1,r0,lf3,r1,l+ζf3,r1,l-141-nζf1,r1,l∂f3,r0,l∂ζ-f3,r1,l∂f1,r0,l∂ζThe variable coefficients for (82) are given by (88)α3,rl,1η,ζ=∂Γ1∂f1,r3,l=1,α2,rl,1η,ζ=∂Γ1∂f1,r2,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α1,rl,1η,ζ=∂Γ1∂f1,r1,l=-n+1f1,r1,l-141-nζ∂f1,r0,l∂ζ,α0,rl,1η,ζ=∂Γ1∂f1,r0,l=14n+3f1,r2,l,βrl,1η,ζ=∂Γ1∂∂f1,r0,l/∂ζ=141-nζf1,r2,l,γrl,1η,ζ=∂Γ1∂∂f1,r1,l/∂ζ=-141-nζf1,r1,l,for (83) are given by(89)γrl,2η,ζ=∂Γ2∂∂f2,r1,l/∂ζ=0,α3,rl,2η,ζ=∂Γ2∂f2,r3,l=0,α2,rl,2η,ζ=∂Γ2∂f2,r2,l=1Pr,α1,rl,2η,ζ=∂Γ2∂f2,r1,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,2η,ζ=∂Γ2∂f2,r0,l=0,βrl,2η,ζ=∂Γ2∂∂f2,r0,l/∂ζ=-141-nζf1,r1,l,and for (84) are given by (90)γrl,3η,ζ=∂Γ3∂∂f3,r1,l/∂ζ=0,α3,rl,3η,ζ=∂Γ3∂f3,r3,l=0,α2,rl,3η,ζ=∂Γ3∂f3,r2,l=1Sc,α1,rl,3η,ζ=∂Γ3∂f3,r1,l=14n+3f1,r0,l+ζ+141-nζ∂f1,r0,l∂ζ,α0,rl,3η,ζ=∂Γ3∂f3,r0,l=0,βrl,3η,ζ=∂Γ3∂∂f3,r0,l/∂ζ=-141-nζf1,r1,l.Equations (79), (80), and (81) can be expressed, respectively, as the following Nζ(Nη+1)×Nζ(Nη+1) matrix system:(91)B0,0l,1B0,1l,1⋯B0,Nζ-1l,1B1,0l,1B1,1l,1⋯B1,Nζ-1l,1⋮⋮⋱⋮BNζ-1,0l,1BNζ-1,1l,1⋯BNζ-1,Nζ-1l,1F1,0lF1,1l⋮F1,Nζ-1l=R1,0lR1,1l⋮R1,Nζ-1l,(92)B0,0l,2B0,1l,2⋯B0,Nζ-1l,2B1,0l,2B1,1l,2⋯B1,Nζ-1l,2⋮⋮⋱⋮BNζ-1,0l,2BNζ-1,1l,2⋯BNζ-1,Nζ-1l,2F2,0lF2,1l⋮F2,Nζ-1l=R2,0lR2,1l⋮R2,Nζ-1l,(93)B0,0l,3B0,1l,3⋯B0,Nζ-1l,3B1,0l,3B1,1l,3⋯B1,Nζ-1l,3⋮⋮⋱⋮BNζ-1,0l,3BNζ-1,1l,3⋯BNζ-1,Nζ-1l,3F3,0lF3,1l⋮F3,Nζ-1l=R3,0lR3,1l⋮R3,Nζ-1l,where(94)Bi,il,1=∑s=03αs,rl,1Ds+βrl,1di,iI+γrl,1di,iD,forl=1,2,…,zwheni=j,Bi,jl,1=βrl,1di,jI+γrl,1di,jD,for,l=1,2,…,zwheni≠j,(95)Bi,il,2=∑s=03αs,rl,2Ds+βrl,2di,iI+γrl,2di,iD,forl=1,2,…,zwheni=j,Bi,jl,2=βrl,2di,jI+γrl,2di,jD,for,l=1,2,…,zwheni≠j,(96)Bi,il,3=∑s=03αs,rl,3Ds+βrl,3di,iI+γrl,3di,iD,forl=1,2,…,zwheni=j,Bi,jl,3=βrl,3di,jI+γrl,3di,jD,for,l=1,2,…,zwheni≠j.The vectors R1,i(l), R2,i(l) and R3,i(l) are defined as (97)R1,il=R1,il-βrl,1di,NζI+γrl,1di,NζDF1,Nζl,R2,il=R2,il-βrl,2di,NζI+γrl,2di,NζDF2,Nζl,R3,il=R3,il-βrl,3di,NζI+γrl,3di,NζDF3,Nζlfor i=0,1,…,Nζ-1.


## 4. Series Solution for the Limiting Case

In this section, the systems of equations presented above are solved using the series solution approach. These series solutions for the limiting case whenξ is very large are used to validate the results obtained using the multidomain bivariate spectral local linearisation method. Equations of the form (58)-(59) and (75)-(87) can be decomposed to a sequence of coupled linear ordinary differential equations when ξ→∞. The solutions of the resulting systems of equations can be determined using elementary methods for solving differential equations.Consider the series solution of (58)-(59). For large values of ξ, the dominant terms in (58) are f′′′ and ξf′′. Balancing the orders of magnitude of these terms gives η=O(ξ-1). Since θ(0)=1, we may infer that θ(η,ξ) is of O(1) as ξ→∞. Assuming that f′′′ has the same order as θ for large ξ, then f=O(ξ-3). Thus for large values of ξ, we define the following transformations:(98)η¯=ξη,f=ξ-3Fη¯,ξ,θ=Θη¯,ξ. The transformation equation (98) reduces (58)-(59) to (99)F′′′+F′′+Θ+1+nξ-4FF′′-nξ-4F′2=1-n4ξ-3F′∂F′∂ξ-F′′∂F∂ξ,(100)1PrΘ′′+Θ′+1+nξ-4FΘ′-nξ-4F′Θ=1-n4ξ-3F′∂Θ∂ξ-Θ′∂F∂ξ,where the primes now denote derivatives with respect to η¯. The corresponding boundary conditions are given by(101)Fξ,0=F′ξ,0=0,Θξ,0=1,F′ξ,∞=Θξ,∞=0.Using the fact that ξ is large, solutions to (99)-(100) are determined in series form as(102)Fξ,η¯=∑k=0∞ξ-4kFkη¯,Θξ,η¯=∑k=0∞ξ-4kΘkη¯.Substituting (102) into (99)-(101) and then equating the coefficients of like powers of ξ, we obtain the equations for k=0 as(103)F0′′′+F0′′+Θ0=0,(104)1PrΘ0′′+Θ0′=0,subject to the following boundary conditions:(105)F00=F0′0=0,Θ00=1,F0′∞=Θ0∞=0. The analytical solutions of the zeroth order equations (103)-(105) are(106)F0η¯=e-Prη¯Pr-1Pr2-1Pr-1Pr2-e-η¯Pr-1Pr+1Pr-1Pr,(107)Θ0η¯=e-Prη¯.For k≥1, we sequentially solve the following system of ordinary differential equations:(108)Fk′′′+Fk′′+Θk=n∑i=0k-1Fk-1-i′Fi′-1+n∑i=0k-1Fk-1-iFi′′+1-n∑i=0k-1iFk-1-i′′Fi-Fk-1-i′Fi′,(109)1PrΘk′′+Θk′=n∑i=0k-1Fk-1-i′Θi-1+n∑i=0k-1Fk-1-iΘi′+1-n∑i=0k-1iΘk-1-i′Fi-Fk-1-i′Θi,subject to the following boundary conditions(110)Fk0=Fk′0=0,Θk0=0,Fk′∞=Θk∞=0.The series solution of (75)-(87) is now considered. For large values of ζ, the dominant terms in (75) are f′′′ and ζf′′. Balancing terms in (75)-(87), for ζ→∞, we get(111)η=Oζ-1,f=Oζ-3,g=O1,h=O1.This leads to the following transformations:(112)f=ζ-3Fη,ζ,η¯=ζη,g=Gη,ζ,h=Hη,ζ.Using the transformations (111)-(112) in (75)-(87), we obtain(113)F′′′+F′′+1-wG+wH+nζ-4FF′′-F′2=1-n4ζ-3F′∂F′∂ζ-F′′∂F∂ζ,(114)1PrG′′+G′+nζ-4FG′=1-n4ζ-3F′∂G∂ζ-G′∂F∂ζ,(115)1ScH′′+H′+nζ-4FH′=1-n4ζ-3F′∂H∂ζ-H′∂F∂ζ.The corresponding boundary conditions are given by(116)Fζ,0=F′ζ,0=0,Gζ,0=1,Hζ,0=1,F′ζ,∞=Gζ,∞=Hζ,∞=0.For large ζ, we seek solutions to (113)-(116) in series form as(117)Fζ,η¯=∑k=0∞ζ-4kFkη¯,Gζ,η¯=∑k=0∞ζ-4kGkη¯,Hζ,η¯=∑k=0∞ζ-4kHkη¯.Substituting (117) into (113)-(116) and then equating the coefficients of like powers of ζ, we obtain the equations for k=0 as(118)F0′′′+F0′′+1-wG0+wH0=0,(119)1PrG0′′+G0′=0,(120)1ScH0′′+H0′=0,subject to the following boundary conditions(121)F00=F0′0=0,G00=1,H00=1,F0′∞=G0∞=H0∞=0.The sequence of equations obtained when k≥1 is given as (122)Fk′′′+Fk′′+1-wGk+wHk=n∑i=0k-1Fk-1-i′Fi′-Fk-1-iFi′′+1-n∑i=0k-1iFk-1-i′′Fi-Fk-1-i′Fi′,(123)1PrGk′′+Gk′=-n∑i=0k-1Fk-1-iGi′+1-n∑i=0k-1iGk-1-i′Fi-Fk-1-i′Gi,(124)1ScHk′′+Hk′=-n∑i=0k-1Fk-1-iHi′+1-n∑i=0k-1iHk-1-i′Fi-Fk-1-i′Hi,subject to the following boundary conditions:(125)Fk0=Fk′0=0,Gk0=1,Hk0=1,Fk′∞=Gk∞=Hk∞=0.


## 5. Results and Discussion

In this section, we present and analyze the results obtained using both the multidomain bivariate spectral local linearisation method (MD-BSLLM) and the series solution. The MD-BSLLM was implemented using MATLAB 2013b. The series solutions were obtained using Mathematica’sNDSolve. The parameters to generate all the solutions are given in captions for the tables and graphs. In this article, the physical quantities of interest in the field of fluid dynamics, the skin friction, Sherwood number, and Nusselt number are presented. These physical quantities were obtained using both methods. We present the values of these quantities for different large values of ξ and ζ. The values of these quantities are in excellent agreement for both methods. This in turn implies that the new approach (MD-BSLLM) presented here can be used to solve other nonsimilar boundary layer partial differential equations. The effect of number of intervals in the solution is also presented in this section. Graphs showing the velocity profiles, temperature profiles, and concentration profiles are also presented for various large values of ξ and ζ. These profiles are in excellent agreement with the results presented by Hussain [2].We first present the solutions of (58) and (59). The physical quantities, skin friction, and the Nusselt number results obtained by solving the governing equations (58)-(59) are given in Table 1. These results were obtained directly by solving the sequences of ordinary differential equations, for the large ξ limiting case, using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 1. Both solutions are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. We observe that an increase in ξ results in a decrease in both the skin friction and the Nusselt number. It is remarkable that this numerical method was able to produce accurate results for large values of ξ with minimum number of grid points. These results were obtained using Nξ=5, Nη=60, and converged after five iterations.Table 1
Comparison of the MD-BSLLM solution and the series solution forf′′(0,ξ) and θ′(0,ξ), with Pr=0.7, n=1/2, and z=40.
Multi-Domain BSLLM Series Solution ξ f ′ ′ ( 0 , ξ ) θ ′ ( 0 , ξ ) f ′ ′ ( 0 , ξ ) θ′(0,ξ) 5 0.2842719 -3.5066677 0.2842719 -3.5066677 10 0.1428115 -7.0008399 0.1428115 -7.0008399 15 0.0952321 -10.5002490 0.09523208 -10.5002490 20 0.0714271 -14.0001050 0.07142714 -14.0001050 25 0.0571424 -17.5000538 0.05714239 -17.5000538 30 0.0476189 -21.0000311 0.04761886 -21.0000311 35 0.0408162 -24.5000196 0.04081624 -24.5000196 40 0.0357142 -28.0000131 0.03571424 -28.0000131Figure1 shows the velocity profiles for different values of ξ. An increase in ξ leads to a decrease in the velocity. The velocity profile is parabolic and hence cannot be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. This velocity profile is in excellent agreement with that obtained by Hussain [2]; thus the MD-BSLLM method can be used as a numerical method for solving large parameter partial differential equations.Figure 1
Velocity profiles at different values ofξ with Pr=0.7, n=1/2, and p=40.Figure2 shows the temperature profile for various values of ξ. The temperature profile behaves like an exponential decay curve, that is, an increase in the pseudosimilarity variable η results in a decrease in the temperature. Temperature profiles were also obtained for large ξ values. Decreasing ξ leads to an increase in the temperature profile. Figures 1 and 2 were obtained using Nξ=5 and Nη=60.Figure 2
Temperature profiles at different values ofξ with Pr=0.7, n=1/2, and z=40.For the second problem, the skin friction, Sherwood number, and Nusselt number results obtained by solving the governing equations (75)-(87) are given in Table 2. Similarly, these results were obtained directly by solving the sequences of ordinary differential equations, for the large ξ, using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 2. The MD-BSLLM and series solutions of (75)-(87) are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. Increasing the transpiration parameter ζ results in a decrease of the skin friction, Sherwood number, and Nusselt number. These results were obtained using few grid points, Nξ=5 and Nη=60.Table 2
Comparison of the multidomain solution and the series solution forf′′(0,ζ), -g′(0,ζ), and -h′(0,ζ): n=w=1/2, Pr=0.7, Sc=0.6, and z=40.
Multi-Domain BSLLM Series Solution ζ f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) 5 0.3088214 3.5018961 3.0018658 0.3088214 3.5018961 3.0018658 10 0.1547399 7.0002370 6.0002332 0.1547399 7.0002370 6.0002332 15 0.1031717 10.5000702 9.0000691 0.1031717 10.5000702 9.0000691 20 0.0773803 14.0000296 12.0000292 0.0773803 14.0000296 12.0000292 25 0.0619045 17.5000152 15.0000149 0.0619045 17.5000152 15.0000149 30 0.0515872 21.0000088 18.0000086 0.0515872 21.0000088 18.0000086 35 0.0442176 24.5000055 21.0000054 0.0442176 24.5000055 21.0000054 40 0.0386905 28.0000037 24.0000036 0.0386905 28.0000037 24.0000036Table3 gives a comparison of the results obtained using the MD-BSLLM and the bivariate spectral local linearisation method (BSLLM) for large values of ζ. The main purpose of the results is to confirm that the accuracy of the BSLLM method deteriorates as ζ increases. We observe that the computed values agree up to an average of three digits. This implies that decomposing the main interval into subintervals and solving the nonsimilar partial differential equations in each subinterval improves the accuracy of the BSLLM method. Table 3 in addition displays the computational time of both methods. For Nη=60, Nζ=5, and z=40, the BSLLM method takes about a fraction of a second to achieve inaccurate results meanwhile the MD-BSLLM method takes about 18 seconds to achieve accurate results.Table 3
Comparison of the MD-BSLLM and the BSLLM solution forf′′(0,ζ), -g′(0,ζ), and -h′(0,ζ): Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.
BSLLM Multi-Domain BSLLM ζ f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) 5 0.3046603 3.5051984 3.0093278 0.3088214 3.5018961 3.0018658 10 0.1546097 7.0003327 6.0003186 0.1547399 7.0002370 6.0002332 15 0.1030889 10.5001437 9.0001125 0.1031717 10.5000702 9.0000691 20 0.0773803 14.0000294 12.0000289 0.0773803 14.0000296 12.0000292 25 0.0619452 17.4999716 14.9999874 0.0619045 17.5000152 15.0000149 30 0.0516135 20.9999793 17.9999897 0.0515872 21.0000088 18.0000086  CPU Time 0.845 0.845 0.845 17.527 17.527 17.527Table4 shows a comparison of the BSLLM method and the MD-BSLLM method with different grid points. For the MD-BSLLM method, we use Nη=60, Nζ=5, and z=40, while for the BSLLM method, we use Nη=40 and Nζ=80 to achieve comparable accurate results. Since the BSLLM method requires more grid points to achieve accurate results over a large domain, then it takes more computational time compared to the MD-BSLLM method over the same interval. In Table 4, the MD-BSLLM method takes about 18 seconds to achieve accurate results; meanwhile the BSLLM method takes about 78 seconds to achieve accurate results. This implies that the MD-BSLLM method takes less computational time compared to the BSLLM method to achieve accurate results.Table 4
Comparison of the MD-BSLLM and the BSLLM solution forf′′(0,ζ), -g′(0,ζ), and -h′(0,ζ): Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.
BSLLM Multi-Domain BSLLM ζ f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) 5 0.3088066 3.5018993 3.0018843 0.3088214 3.5018961 3.0018658 10 0.1547399 7.0002370 6.0002332 0.1547399 7.0002370 6.0002332 15 0.1031717 10.5000702 9.0000691 0.1031717 10.5000702 9.0000691 20 0.0773801 14.0000296 12.0000292 0.0773803 14.0000296 12.0000292 25 0.0619014 17.5000152 15.0000149 0.0619045 17.5000152 15.0000149 30 0.0515660 21.0000087 18.0000086 0.0515872 21.0000088 18.0000086  CPU Time 77.921 77.921 77.921 17.527 17.527 17.527In Table5, we observe that, for small values of ζ, the results from both methods (BSLLM and MD-BSLLM) are in good agreement. The BSLLM method takes about 3 seconds while the MD-BSLLM method takes about 22 seconds to achieve accurate results. Therefore, we can conclude that it is not necessary to use the MD-BSLLM method to solve the partial differential equations in small subdomains of ζ. The multidomain approach is suitable for solving equations with large values of ζ. The results in Table 5 were obtained with Nη=60 and Nζ=15 for the BSLLM method; meanwhile Nη=60, Nζ=5, and z=40 were used for the MD-BSLLM method for ζ=2.Table 5
Comparison of the MD-BSLLM and the BSLLM solution forf′′(0,ζ), -g′(0,ζ), and -h′(0,ζ):n=w=1/2, Pr=0.7, and Sc=0.6.
BSLLM Multi-Domain BSLLM ζ f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) 0.25 0.6346269 0.6230281 0.6046615 0.6346269 0.6230281 0.6046615 0.50 0.6393410 0.7247205 0.6899358 0.6393410 0.7247205 0.6899358 0.75 0.6379544 0.8353765 0.7818927 0.6379544 0.8353765 0.7818927 1.00 0.6307080 0.9546061 0.8803437 0.6307080 0.9546061 0.8803437 1.25 0.6181647 1.0818907 0.9850250 0.6181647 1.0818907 0.9850250 1.50 0.6011394 1.2166019 1.0955986 0.6011394 1.2166019 1.0955986 1.75 0.5806016 1.3580272 1.2116586 0.5806016 1.3580272 1.2116586 2.00 0.5575697 1.5054016 1.3327435 0.5575697 1.5054016 1.3327435  CPU Time 3.088 3.088 3.088 22.436 22.436 22.436Table6 shows the effect of the number of subintervals z. Table 6 shows that increasing the number of subintervals yields more accurate results. We note that when z=1, that is, we have one domain, then the MD-BSLLM is behaves exactly as the BSLLM method. For z=1, the results do not match because the MD-BSLLM is exactly the BSLLM method. However, for ζ=40, only ten subintervals ensure accuracy of the MD-BSLLM method as depicted in Table 6. These results confirm that decomposing the main domain into subintervals improves the accuracy of the BSLLM method.Table 6
Comparison of the multidomain solution and the series solution forf′′(0,ζ), -g′(0,ζ), and -h′(0,ζ): Nη=60, Nζ=5, η∞, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.
Multi-Domain BSLLM Series Solution z f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) f ′ ′ ( 0 , ζ ) - g ′ ( 0 , ζ ) - h ′ ( 0 , ζ ) 1 0.0967112 11.2002848 9.6003085 0.0967241 11.2000579 9.6000569 5 0.0967241 11.2000577 9.6000568 0.0967241 11.2000579 9.6000569 10 0.0967241 11.2000579 9.6000569 0.0967241 11.2000579 9.6000569 15 0.0967241 11.2000579 9.6000569 0.0967241 11.2000579 9.6000569 20 0.0967241 11.2000579 9.6000569 0.0967241 11.2000579 9.6000569Figures3, 4, and 5 show the velocity, temperature, and concentration profiles for various values of ζ, respectively. It is observed that an increase in ζ leads to a decrease in the velocity, temperature, and concentration. The velocity profile is parabolic and hence difficult to be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. The temperature and concentration profiles decay exponentially. The velocity, temperature, and concentration profiles are in excellent agreement with the results by Hossain [1]. This implies that the MD-BSLLM method can be used as a numerical method for solving large parameter nonsimilar, nonlinear partial differential equations.Figure 3
Velocity profiles at different values ofζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.Figure 4
Temperature profiles at different values ofζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.Figure 5
Concentration profiles at different values ofζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.Figures6, 7, and 8 show the convergence error graphs for the velocity, temperature, and concentration profiles for various values of ζ, respectively. It is observed that the method converges after approximately three iterations to a very small error. This validates the accuracy and convergence of the method.Figure 6
Convergence error graph off(ζ,ξ) at different values of ζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.Figure 7
Convergence error graph ofg(ζ,ξ) at different values of ζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.Figure 8
Convergence error graph ofh(ζ,ξ) at different values of ζ, Nη=60, Nζ=5, n=w=1/2, Pr=0.7, Sc=0.6, and z=40.


## 6. Conclusion

In this article, a new numerical approach for solving systems ofn nonsimilar nonlinear partial differential equations over large domains is presented and implemented. The method is generalized for any pth order system of n partial differential equations. The method was tested on two different numerical experiments. The results were validated by using series solution of the numerical experiments. The results of the numerical method and the series solution were in excellent agreement and hence validating the accuracy of the method over large domains. The idea of decomposing the main domain into subintervals increases the accuracy of the method. The more subintervals we have, the more accurate the method is. However, for small domains, it is recommended to use and implement the BSLLM method other than the MD-BSLLM method since the BSLLM method achieves accurate results within few seconds as opposed to the MD-BSLLM method. This project has added to literature a novel approach for solving nonsimilar nonlinear system of partial differential equations over a large domain. In future, we will develop and solve systems of nonlinear partial differential equations arising from other disciplines.



---

*Source: 6423294-2019-06-10.xml*