# A New Spectral Local Linearization Method for Nonlinear Boundary Layer Flow Problems

Authors: S. S. Motsa
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We propose a simple and efficient method for solving highly nonlinear systems of boundary layer flow problems with exponentially decaying profiles. The algorithm of the proposed method is based on an innovative idea of linearizing and decoupling the governing systems of equations and reducing them into a sequence of subsystems of differential equations which are solved using spectral collocation methods. The applicability of the proposed method, hereinafter referred to as the spectral local linearization method (SLLM), is tested on some well-known boundary layer flow equations. The numerical results presented in this investigation indicate that the proposed method, despite being easy to develop and numerically implement, is very robust in that it converges rapidly to yield accurate results and is more efficient in solving very large systems of nonlinear boundary value problems of the similarity variable boundary layer type. The accuracy and numerical stability of the SLLM can further be improved by using successive overrelaxation techniques.


---

## Body

## 1. Introduction

The important properties and applications of the boundary layer flows in many engineering areas are well documented and have been widely researched in the last few decades. The mathematical model of the boundary layer flow system is described, in general, by a set of nonlinear partial differential equations which cannot be solved exactly. A common approach that is employed in the solution of boundary layer flow problems is the reduction of the governing partial differential equations into a set of ordinary differential equations (ODEs) after using appropriate assumptions and a set of transformations based on the so-called similarity variables and using the so-called stream function formulations. The transformed ODEs with the corresponding boundary conditions becomes a system of highly nonlinear boundary value problems (BVP) which, in most cases, does not have closed form solutions. The solution of these nonlinear BVPs continues to fascinate and inspire researchers to develop methods of obtaining solutions to these equations which elucidate the intricate properties of the underlying boundary layer flow problem under different conditions.The standard way of solving the transformed similarity variable boundary layer equations is the numerical approach based on the shooting algorithm with the Runge-Kutta scheme [1, 2]. Other numerical approaches that have been preferred by some researchers include the finite difference method [3], Keller-Box method [4, 5], spectral homotopy analysis method [6, 7], element free Galerkin method [8], and Newton-Raphson based methods such as the quasilinearization method of [9, 10] and the successive linearization method [11–13]. Some analytical approaches have also been found to be very useful in solving boundary layer equations. Examples include the homotopy analysis method [14–16], the homotopy perturbation method [17–19], the differential transformation method [20–22], the Lie-group shooting method [23], the parameter iteration method [24], and the variational iteration method [25], among many other approaches. However, analytical methods are limited in their applications and can only be used in simple systems with few equations. Most real life applications in fluid mechanics involve multiple interacting physical processes and are modelled using complex multiple equation systems. For such problems, it is impractical to use analytical methods because their solution process becomes too cumbersome, and convergence to the true solution can be very slow or not possible at all. For, this reason numerical methods are by far the most practical way of seeking solutions to the boundary layer flow type of highly nonlinear systems. Numerical solution methods also have their own disadvantages such as the stability and convergence issues, difficulty in dealing with singularities, and limiting cases. In addition, some numerical solutions are difficult to interpret such as the occurrence of multiple solutions. The quest for the most optimal method of solving nonlinear problems in fluid mechanics is what drives, ever growing interest in the development of new methods and the modification and improvement of existing analytical and numerical methods.The prime objective of this paper is to present a new numerical method of solving boundary layer equations that seeks to address some of the aforementioned numerical difficulties. We propose a very simple, yet very accurate and convergent iterative algorithm for solving nonlinear systems of equations that model boundary layer flow problems. The proposed method, hereinafter referred to as spectral local linearization method (SLLM), is based on, decoupling and linearizing systems of equations using a combination of a univariate linearization technique and a spectral collocation discretization. The key feature of the SLLM algorithm is that it breaks down a large coupled system of equations into a sequence of smaller subsystems which can be solved sequentially in a very computationally efficient manner. The applicability of the proposed method is tested on the well-known Blasius boundary layer problem and a three-equation coupled system that models the problem of unsteady free convective heat and mass transfer which are used as test cases to compare the performance of the proposed SLLM when compared to other existing methods. The computed SLLM results demonstrate that the method is easy to develop, accurate, convergent, stable, and very efficient when compared with other existing methods of solving some large systems of boundary value problems.


## 2. Description of the Methods of Solution

This section presents a brief description of how the proposed iterative methods of solution are developed for a general system ofm nonlinear ordinary differential equations in m unknown functions.


### 2.1. Spectral Local Linearization Method (SLLM)

Here, we describe the development of the spectral local linearization method (SLLM). Consider a system ofm nonlinear ordinary differential equations in m unknown functions Zi(η)i=1,2,…,m, where η∈[a,b] is the dependent variable. The system can be written in terms of Zi as a sum of its linear (Li) and nonlinear components (Ni) as
(1)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=0,i=1,…,m.To develop the iteration scheme, we apply local linearization ofNi about Zi,r (the previous iteration) to the ith nonlinear equation assuming that all other Zk,r(k≠i) are known. Thus, at the ith equation, Ni is linearized as follows:
(2)Ni[Z1,Z2,…,Zm]=Ni[Z1,r,Z2,r,…,Zm,r]+∂Ni∂Zi[Z1,r,Z2,r,…,Zm,r](Zi-Zi,r).
Thus, at the current iteration with Zi=Zi,r+1, (1) becomes
(3)Li[Z1,r+1,…,Zm,r+1]+∂Ni∂Zi[⋯]Zi,r+1=∂Ni∂Zi[⋯]Zi,r-Ni[Z1,r,…,Zm,r],
where [⋯] denotes [Z1,r,Z2,r,…,Zm,r] and Zi,r+1 and Zi,r are the approximations of Zi at the current and the previous iteration, respectively. To obtain a decoupled iteration scheme, we appeal to the Gauss-Seidel approach of decoupling linear algebraic systems in linear algebra applications. We therefore arrange the equations in a particular order and solve them in a chronological order. In seeking the solution of Zi in the current iteration level, Zi,r+1, we use updated solutions of Zs (s<i) obtained as solutions of the previous i=1,2,…,s equations. Thus, for a system of m equations, the local linearization iteration scheme becomes
(4)L1[Z1,r+1,Z2,r,…,Zm,r]+∂N1∂Z1[⋯]Z1,r+1=∂N1∂Z1[⋯]Z1,r-N1[Z1,r,…,Zm,r],L2[Z1,r+1,Z2,r+1,Z3,r,…,Zm,r]+∂N2∂Z2[⋯]Z2,r+1=∂N2∂Z2[⋯]Z2,r-N2[Z1,r+1,Z2,r,…,Zm,r],⋮Lm[Z1,r+1,Z2,r+1,…,Zm,r+1]+∂Nm∂Zm[⋯]Zm,r+1=∂Nm∂Zm[⋯]Zm,r-Nm[Z1,r+1,…,Zm-1,r+1,Zm,r],
where, at the ith equation, [⋯]≡[Z1,r+1,Z2,r+1,…,Zi-1,r+1,Zi,r,…,Zm,r].Thus, starting from an initial approximationZ1,0,Z2,0,…,Zm,0, the proposed iterative scheme (4) is then solved as a loop until the system converges at a consistent solution for all the variables. To solve the iteration scheme (4), it is convenient to use the Chebyshev pseudospectral method. For this reason the proposed method is referred to as the spectral local linearization iteration method (SLLM) in this work. Spectral methods are now becoming the preferred tools for solving ordinary and partial differential equations because of their elegance and high accuracy in resolving problems with smooth functions.For brevity, we omit the details of the spectral methods and refer interested readers to [26, 27]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [-1,1] on which the spectral method can be implemented. We use the transformation η=(b-a)(τ+1)/2 to map the interval [a,b] on [-1,1]. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix D which is used to approximate the derivatives of the unknown variables Zi(η) at the collocation points as the matrix vector product
(5)dZidη=∑k=0N-‍DlkZi(τk)=DZi,l=0,1,…,N-,
where N-+1 is the number of the collocation points (grid points), D=2D/(b-a), and Z=[z(τ0),z(τ1),…,z(τN)]T is the vector function at the collocation points. Higher order derivatives are obtained as powers of D, that is
(6)Zj(p)=DpZj,
where p is the order of the derivative.


### 2.2. Spectral Quasilinearization Method (SQLM)

In this section, we describe a quasilinearization method that does not use the decoupling approach described in the previous section. This quasilinearization method (QLM) is a generalisation of the Newton-Raphson method and was first proposed by Bellman and Kalaba [28] for solving nonlinear boundary value problems. To obtain the QLM iteration scheme the nonlinear component of a differential equation is linearized using the multivariable Taylor series expansion as opposed to the local linearization approach used in the previous section. The QLM scheme is solved using the Chebyshev pseudospectral method as described in the previous section. Consequently, the method is referred to as the spectral quasilinearization method (SQLM) in this work.To develop the SQLM scheme we consider a system ofm nonlinear ordinary differential equations in m unknowns functions zi(η)i=1,2,…,m, where η is the dependent variable. The system can be written as a sum of its linear ℒ and nonlinear components N as
(7)ℒ[z1(η),z2(η),…,zm(η)]+𝒩[z1(η),z2(η),…,zm(η)]=0,(η)∈(a,b),
subject to the separated boundary conditions
(8)Ai[z1(a),z2(a),…,zm(a)]=Ka,i,Bb[z1(b),z2(b),…,zm(b)]=Kb,i,
where Ai and Bi are linear operators and Ka,i and Kb,i are constants for i=1,2,…,m. Define the vector Zi to be the vector of the derivatives of the variable zi with respect to the dependent variable η, that is,
(9)Zi=[zi(0),zi(1),…,zi(ni)],
where zi(0)=zi,  zi(p) is the pth derivative of zi with respect to η and ni(i=1,2,…,m) is the highest derivative order of the variable zi appearing in the system of equations. In addition, we define Li and Ni to be the linear and nonlinear operators, respectively, that operate on the Zi for i=1,2,…,m. With these definitions, (7) and (8) can be written as
(10)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=∑j=1m∑p=0nj‍αi,j[p]zj(p)+Ni[Z1,Z2,…,Zm]=0,
where αi,j[p] are the constant coefficients of zj(p), and the derivative of zj (j=1,2,…,m) that appears in the ith equation for i=1,2,…,m.Noting that, for each variablezi, the derivatives in the boundary conditions can at most be one less than the highest derivative of zi in the governing system (7), we define the vector Z~i to be the vector of the derivatives of the variable zi with respect to the dependent variable η from 0 up to ni-1, that is,
(11)Z~i=[zi(0),zi(1),…,zi(ni-1)].
The boundary conditions (8) can be written as
(12)Aν[Z~1(a),Z~2(a),…,Z~m(a)]=∑j=1m‍∑p=0nj-1‍βν,j[p]zj(p)(a)=Ka,ν,ν=1,2,…,ma,Bσ[Z~1(b),Z~2(b),…,Z~m(b)]=∑j=1m‍∑p=0nj-1‍γσ,j[p]zj(p)(b)=Kb,σ,σ=1,2,…,mb,
where βν,j[p] (γσ,j[p]) are the constant coefficients of zj(p) in the boundary conditions and ma, mb are the total number of prescribed boundary conditions at x=a and x=b, respectively. We remark that the sum ma+mb is equal to the sum of the highest orders of the derivatives corresponding to the dependent variables zi, that is
(13)ma+mb=∑i=1m‍ni.Assume that the solutionzi(η) of (10) at the (r+1)th iteration is zi,r+1. If the solution at the previous iteration zi,r(η) is sufficiently close to zi,r+1, the nonlinear component Ni of (10) can be linearized using the one-term Taylor series for multiple variables so that (10) can be approximated as
(14)Li[Z1,r+1,…,Zm,r+1]+Ni[⋯]+∑j=1m‍∑p=0nj‍(zj,r+1(p)-zj,r(p))∂Ni∂zj(p)[⋯]=0,
subject to
(15)∑j=1m‍∑p=0nj-1‍βν,j[p]zj,r+1(p)(a)=0,ν=1,2,…,ma,∑j=1m‍∑p=0nj-1‍γσ,j[p]zj,r+1(p)(b)=0,σ=1,2,…,mb,
where
(16)[⋯]=[Z1,r,Z2,r,…,Zm,r].
Equation (14) can be rewritten as
(17)Li[Z1,r+1,…,Zm,r+1]+∑j=1m‍∑p=0nj‍zj,r+1(p)∂Ni∂zj(p)[⋯]=∑j=1m‍∑p=0nj‍zj,r(p)∂Ni∂zj(p)[⋯]-Ni[⋯].The initial approximation,zj,0(η), required to start the iteration scheme (17) is chosen to be a function that satisfies the boundary conditions (8). Applying the Chebyshev spectral collocation on the recursive iteration scheme (17) gives
(18)∑j=1m‍[Λi,j+Πi,j]Zj,r+1=Φi,r,i,j=1,2,…,m,
where Zi,r+1=[zi,r+1(τ0),zi,r+1(τ1),…,zi,r+1(τN)]T, Λi,j, Πi,j and Φi are given by
(19)Λi,j=∑p=0nj‍αi,jpDp,Πi,j=∑p=0nj‍[∂Ni∂zj(p)]dDp,i,j=1,2,…,m,Φi,r=∑j=1m‍∑p=0nj‍zj,r(p)∂Ni∂zj(p)[⋯]-Ni[Z1,r,Z2,r,…,Zm,r],
respectively, where [⋯]d means denote diagonal matrix. Defining Δ=Λ+Π, we can write (18) in matrix form as
(20)[Δ1,1Δ1,2⋯Δ1,mΔ2,1Δ2,2⋯Δ2,m⋮⋮⋮Δm,1Δm,2⋯Δm,m][Z1,r+1Z2,r+1⋮Zm,r+1]=[Φ1,rΦ2,r⋮Φm,r],
where Zi,r, Φi,r are vectors of size (N-+1)×1 and Δi,j are (N-+1)×(N-+1) matrices. Starting from Zi,0, the recursive sequence (20) is solved iteratively for r=0,1,2,3,….Thus, the size of the coefficient matrix in (20) is m(N-+1)×m(N-+1) and the column vector on the right-hand side has the dimension m(N-+1)×1.


## 2.1. Spectral Local Linearization Method (SLLM)

Here, we describe the development of the spectral local linearization method (SLLM). Consider a system ofm nonlinear ordinary differential equations in m unknown functions Zi(η)i=1,2,…,m, where η∈[a,b] is the dependent variable. The system can be written in terms of Zi as a sum of its linear (Li) and nonlinear components (Ni) as
(1)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=0,i=1,…,m.To develop the iteration scheme, we apply local linearization ofNi about Zi,r (the previous iteration) to the ith nonlinear equation assuming that all other Zk,r(k≠i) are known. Thus, at the ith equation, Ni is linearized as follows:
(2)Ni[Z1,Z2,…,Zm]=Ni[Z1,r,Z2,r,…,Zm,r]+∂Ni∂Zi[Z1,r,Z2,r,…,Zm,r](Zi-Zi,r).
Thus, at the current iteration with Zi=Zi,r+1, (1) becomes
(3)Li[Z1,r+1,…,Zm,r+1]+∂Ni∂Zi[⋯]Zi,r+1=∂Ni∂Zi[⋯]Zi,r-Ni[Z1,r,…,Zm,r],
where [⋯] denotes [Z1,r,Z2,r,…,Zm,r] and Zi,r+1 and Zi,r are the approximations of Zi at the current and the previous iteration, respectively. To obtain a decoupled iteration scheme, we appeal to the Gauss-Seidel approach of decoupling linear algebraic systems in linear algebra applications. We therefore arrange the equations in a particular order and solve them in a chronological order. In seeking the solution of Zi in the current iteration level, Zi,r+1, we use updated solutions of Zs (s<i) obtained as solutions of the previous i=1,2,…,s equations. Thus, for a system of m equations, the local linearization iteration scheme becomes
(4)L1[Z1,r+1,Z2,r,…,Zm,r]+∂N1∂Z1[⋯]Z1,r+1=∂N1∂Z1[⋯]Z1,r-N1[Z1,r,…,Zm,r],L2[Z1,r+1,Z2,r+1,Z3,r,…,Zm,r]+∂N2∂Z2[⋯]Z2,r+1=∂N2∂Z2[⋯]Z2,r-N2[Z1,r+1,Z2,r,…,Zm,r],⋮Lm[Z1,r+1,Z2,r+1,…,Zm,r+1]+∂Nm∂Zm[⋯]Zm,r+1=∂Nm∂Zm[⋯]Zm,r-Nm[Z1,r+1,…,Zm-1,r+1,Zm,r],
where, at the ith equation, [⋯]≡[Z1,r+1,Z2,r+1,…,Zi-1,r+1,Zi,r,…,Zm,r].Thus, starting from an initial approximationZ1,0,Z2,0,…,Zm,0, the proposed iterative scheme (4) is then solved as a loop until the system converges at a consistent solution for all the variables. To solve the iteration scheme (4), it is convenient to use the Chebyshev pseudospectral method. For this reason the proposed method is referred to as the spectral local linearization iteration method (SLLM) in this work. Spectral methods are now becoming the preferred tools for solving ordinary and partial differential equations because of their elegance and high accuracy in resolving problems with smooth functions.For brevity, we omit the details of the spectral methods and refer interested readers to [26, 27]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [-1,1] on which the spectral method can be implemented. We use the transformation η=(b-a)(τ+1)/2 to map the interval [a,b] on [-1,1]. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix D which is used to approximate the derivatives of the unknown variables Zi(η) at the collocation points as the matrix vector product
(5)dZidη=∑k=0N-‍DlkZi(τk)=DZi,l=0,1,…,N-,
where N-+1 is the number of the collocation points (grid points), D=2D/(b-a), and Z=[z(τ0),z(τ1),…,z(τN)]T is the vector function at the collocation points. Higher order derivatives are obtained as powers of D, that is
(6)Zj(p)=DpZj,
where p is the order of the derivative.


## 2.2. Spectral Quasilinearization Method (SQLM)

In this section, we describe a quasilinearization method that does not use the decoupling approach described in the previous section. This quasilinearization method (QLM) is a generalisation of the Newton-Raphson method and was first proposed by Bellman and Kalaba [28] for solving nonlinear boundary value problems. To obtain the QLM iteration scheme the nonlinear component of a differential equation is linearized using the multivariable Taylor series expansion as opposed to the local linearization approach used in the previous section. The QLM scheme is solved using the Chebyshev pseudospectral method as described in the previous section. Consequently, the method is referred to as the spectral quasilinearization method (SQLM) in this work.To develop the SQLM scheme we consider a system ofm nonlinear ordinary differential equations in m unknowns functions zi(η)i=1,2,…,m, where η is the dependent variable. The system can be written as a sum of its linear ℒ and nonlinear components N as
(7)ℒ[z1(η),z2(η),…,zm(η)]+𝒩[z1(η),z2(η),…,zm(η)]=0,(η)∈(a,b),
subject to the separated boundary conditions
(8)Ai[z1(a),z2(a),…,zm(a)]=Ka,i,Bb[z1(b),z2(b),…,zm(b)]=Kb,i,
where Ai and Bi are linear operators and Ka,i and Kb,i are constants for i=1,2,…,m. Define the vector Zi to be the vector of the derivatives of the variable zi with respect to the dependent variable η, that is,
(9)Zi=[zi(0),zi(1),…,zi(ni)],
where zi(0)=zi,  zi(p) is the pth derivative of zi with respect to η and ni(i=1,2,…,m) is the highest derivative order of the variable zi appearing in the system of equations. In addition, we define Li and Ni to be the linear and nonlinear operators, respectively, that operate on the Zi for i=1,2,…,m. With these definitions, (7) and (8) can be written as
(10)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=∑j=1m∑p=0nj‍αi,j[p]zj(p)+Ni[Z1,Z2,…,Zm]=0,
where αi,j[p] are the constant coefficients of zj(p), and the derivative of zj (j=1,2,…,m) that appears in the ith equation for i=1,2,…,m.Noting that, for each variablezi, the derivatives in the boundary conditions can at most be one less than the highest derivative of zi in the governing system (7), we define the vector Z~i to be the vector of the derivatives of the variable zi with respect to the dependent variable η from 0 up to ni-1, that is,
(11)Z~i=[zi(0),zi(1),…,zi(ni-1)].
The boundary conditions (8) can be written as
(12)Aν[Z~1(a),Z~2(a),…,Z~m(a)]=∑j=1m‍∑p=0nj-1‍βν,j[p]zj(p)(a)=Ka,ν,ν=1,2,…,ma,Bσ[Z~1(b),Z~2(b),…,Z~m(b)]=∑j=1m‍∑p=0nj-1‍γσ,j[p]zj(p)(b)=Kb,σ,σ=1,2,…,mb,
where βν,j[p] (γσ,j[p]) are the constant coefficients of zj(p) in the boundary conditions and ma, mb are the total number of prescribed boundary conditions at x=a and x=b, respectively. We remark that the sum ma+mb is equal to the sum of the highest orders of the derivatives corresponding to the dependent variables zi, that is
(13)ma+mb=∑i=1m‍ni.Assume that the solutionzi(η) of (10) at the (r+1)th iteration is zi,r+1. If the solution at the previous iteration zi,r(η) is sufficiently close to zi,r+1, the nonlinear component Ni of (10) can be linearized using the one-term Taylor series for multiple variables so that (10) can be approximated as
(14)Li[Z1,r+1,…,Zm,r+1]+Ni[⋯]+∑j=1m‍∑p=0nj‍(zj,r+1(p)-zj,r(p))∂Ni∂zj(p)[⋯]=0,
subject to
(15)∑j=1m‍∑p=0nj-1‍βν,j[p]zj,r+1(p)(a)=0,ν=1,2,…,ma,∑j=1m‍∑p=0nj-1‍γσ,j[p]zj,r+1(p)(b)=0,σ=1,2,…,mb,
where
(16)[⋯]=[Z1,r,Z2,r,…,Zm,r].
Equation (14) can be rewritten as
(17)Li[Z1,r+1,…,Zm,r+1]+∑j=1m‍∑p=0nj‍zj,r+1(p)∂Ni∂zj(p)[⋯]=∑j=1m‍∑p=0nj‍zj,r(p)∂Ni∂zj(p)[⋯]-Ni[⋯].The initial approximation,zj,0(η), required to start the iteration scheme (17) is chosen to be a function that satisfies the boundary conditions (8). Applying the Chebyshev spectral collocation on the recursive iteration scheme (17) gives
(18)∑j=1m‍[Λi,j+Πi,j]Zj,r+1=Φi,r,i,j=1,2,…,m,
where Zi,r+1=[zi,r+1(τ0),zi,r+1(τ1),…,zi,r+1(τN)]T, Λi,j, Πi,j and Φi are given by
(19)Λi,j=∑p=0nj‍αi,jpDp,Πi,j=∑p=0nj‍[∂Ni∂zj(p)]dDp,i,j=1,2,…,m,Φi,r=∑j=1m‍∑p=0nj‍zj,r(p)∂Ni∂zj(p)[⋯]-Ni[Z1,r,Z2,r,…,Zm,r],
respectively, where [⋯]d means denote diagonal matrix. Defining Δ=Λ+Π, we can write (18) in matrix form as
(20)[Δ1,1Δ1,2⋯Δ1,mΔ2,1Δ2,2⋯Δ2,m⋮⋮⋮Δm,1Δm,2⋯Δm,m][Z1,r+1Z2,r+1⋮Zm,r+1]=[Φ1,rΦ2,r⋮Φm,r],
where Zi,r, Φi,r are vectors of size (N-+1)×1 and Δi,j are (N-+1)×(N-+1) matrices. Starting from Zi,0, the recursive sequence (20) is solved iteratively for r=0,1,2,3,….Thus, the size of the coefficient matrix in (20) is m(N-+1)×m(N-+1) and the column vector on the right-hand side has the dimension m(N-+1)×1.


## 3. Numerical Experiments

In this section, we discuss the implementation of the SLLM and SQLM approaches on two examples of boundary layer flow systems. We consider the Blasius boundary layer problem and a three-equation system that models the problem of unsteady free convective heat and mass transfer on a stretching surface in a porous medium in the presence of a chemical reaction.


### 3.1. Blasius Boundary Layer Equation

The Blasius boundary layer flow equation, in a dimensionless form, is expressed in terms of the similarity variableη∈[0,∞) and function f(η) as
(21)f′′′(η)+12f(η)f′′(η)=0,f(0)=0,f′(0)=0,f′(∞)=1,
where the primes denote differentiation with respect to η. To apply the proposed SLLM iterative method on (21), it is convenient to reduce the order of the governing equations by one through setting f′(η)=g(η). This results in the following system:
(22)f′=g,f(0)=0,(23)g′′+12f(η)g′(η)=0,g(0)=0,g(∞)=1.We observe that (23) is essentially a differential equation whose unknown is g(η), but the equation requires f(η) which can be considered to be a known input function. Thus, assuming that f(η) is known at a particular iteration level, the local linearization iteration scheme can be expressed by the following iteration formula;
(24)gr+1′′+12frgr+1′=0,gr+1(0)=0,gr+1(∞)=1,fr+1′=gr+1,fr+1(0)=0.Applying the Chebyshev pseudospectral method on (24), we obtain the following decoupled system of equations:
(25)A1gr+1=B1,gr+1(τN-)=0,gr+1(τ0)=1,Dfr+1=gr+1,fr+1(τN-)=0,
where
(26)A1=D2+12diag[fr]D,B1=O1,
where O1 is an (N-+1)×1 vector on which the boundary conditions are imposed and f and g are the values of functions f and g, respectively, when evaluated at the collocation points. Equation (25) constitutes the SLLM iteration scheme.The initial approximation required to start the iteration process can be chosen as a function that satisfies the governing boundary conditions and from known physical considerations of the flow properties. For most boundary layer flow problems it is well known that the flow velocity tends to the mainstream flow exponentially. Thus, a suitable initial approximation for the boundary layer problem (21) is
(27)f0(η)=η-1+e-η,g0(η)=1-e-η.To obtain the SQLM iterative scheme for the Blasius problem (21), we set
(28)L=f′′′,N=12ff′′,[f,f′,f′′,f′′′]=[z,z(1),z(2),z(3)].
Thus, using (17), we obtain
(29)fr+1′′′+12frfr+1′′+12fr′′fr+1=12fr′′fr,
subject to
(30)fr+1(0)=0,fr+1′(0)=0,fr+1′(∞)=1.Applying the spectral method, with derivative matrices on the iteration scheme (29) and the corresponding boundary conditions, gives the following matrix system:
(31)ArFr+1=Φr,
with the boundary conditions
(32)fr+1(τN)=0,∑k=0N-‍DN-kfr+1(τk)=0,∑k=0N-‍D0kfr+1(τk)=1,
where
(33)Ar=D3+12diag(fr)D2+12diag(fr′′),Φr=12frfr′′.


### 3.2. Unsteady Free Convective Heat and Mass Transfer on a Stretching Surface in a Porous Medium with Suction/Injection

In this section we consider a three-equation system that models the problem of unsteady free convective heat and mass transfer on a stretching surface in a porous medium in the presence of a chemical reaction. The governing equations [29, 30] for this problem are given as the following dimensionless system of equations and boundary conditions:
(34)f′′′+ff′′-(f′)2-Kf′-A(f′+η2f′′)+Grθ+Gcϕ=0,1Prθ′′-f′θ+fθ′-A(θ+12ηθ′)=0,1Scϕ′′-f′ϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0,(35)f(0)=fw,f′(0)=1,θ(0)=1,ϕ(0)=1,f′(∞)=0,θ(∞)=0,ϕ(∞)=0,
where f(η), θ(η), and ϕ(η) are, respectively, the dimensionless velocity, temperature, and concentration, fw is the suction/injection parameter, γ is the chemical reaction constant, Pr is the Prandtl number, Sc is the Schmidt number, K is the permeability parameter, and Gr and Gc are the temperature and concentration dependent Grashof numbers, respectively.To apply the proposed SLLM iterative method on (34), we set f′(η)=g(η) to obtain
(36)f′=g,(37)g′′+fg′′-g2-Kg-A(g+η2g′)+Grθ+Gcϕ=0,(38)1Prθ′′-gθ+fθ′-A(θ+12ηθ′)=0,(39)1Scϕ′′-gϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0.We observe that the energy equation (38) requires only f and g to resolve the solution of θ. The same applies to the mass transfer equation (39) for ϕ. Equation (37), on the other hand, requires the solutions of f,  θ and ϕ. In developing the SLLM iteration scheme we solve (38) and (39) for θ, and ϕ first, then use their updated solutions in the (37) for g. Thus, assuming that f(η) and g(η) are known at a particular iteration level, the local linearization iteration algorithm can be expressed as
(40)1Prθr+1′′-grθr+1+frθr+1′-A(θr+1+η2θr+1′)=0,1Scϕr+1′′-grϕr+1+frϕr+1′-(A+γ)ϕr+1-Aη2ϕr+1′=0,gr+1′′+frgr+1′-(A+K)gr+1+ηA2gr+1′-2grgr+1+Grθr+1+Gcϕr+1=-gr2,fr+1′=gr+1,fr+1(0)=fw,
subject to
(41)gr+1(0)=1,gr+1(∞)=0,θr+1(0)=1,θr+1(∞)=0,ϕr+1(0)=1,ϕr+1(∞)=0.Applying the Chebyshev pseudospectral method on (40)-(41), we obtain the following decoupled matrix system of equations:
(42)A1Θr+1=B1,θr+1(τN-)=1,θr+1(τ0)=0,A2Φr+1=B2,ϕr+1(τN-)=1,ϕr+1(τ0)=0,A3gr+1=B3-GrΘr+1-GcΦr+1,gr+1(τN-)=1,gr+1(τ0)=0,Dfr+1=gr+1,fr+1(τN-)=fw,
where
(43)A1=1PrD2+diag[fr-Aη2]D-diag[gr]-AI,B1=O1,A2=1ScD2+diag[fr-Aη2]D-diag[gr]-(A+γ)I,B2=O1,A3=D2+diag[fr-Aη2]D-diag[2gr]-(K+A)I,B3=-gr2,
where I is an (N-+1)×(N-+1) identity matrix, O1 is an (N-+1)×1 vector, and Θ, Φ,  g,  and  f are the approximate values of θ, ϕ, g,  and  f evaluated at the collocation points. After incorporating the boundary conditions in the matrices in (42), the approximate solutions of θ,  ϕ,  g, and f at each iteration level can be obtained by solving (42) in the order in which they are listed. A suitable initial approximation for the SLLM scheme is
(44)f0(η)=fw+1-e-η,θ0(η)=e-η,ϕ0(η)=e-η.To obtain the SQLM iterative scheme for the system (34), we set
(45)L1=f′′′-A(f′+η2f′′)+Grθ+Gcϕ,N1=ff′′-(f′)2,L2=1Prθ′′-A(θ+12ηθ′),N2=-f′θ+fθ′,L3=1Scϕ′′-A(ϕ+12ηϕ′)-γϕ,N3=-f′ϕ+fϕ′.To use formula (17) we define the following variables:
(46)[f,f′,f′′,f′′′]=[z1,z1(1),z1(2),z1(3)],[θ,θ′,θ′′]=[z2,z2(1),z2(2)],[ϕ,ϕ′,ϕ′′]=[z3,z3(1),z3(2)].Thus, using (45)-(46) in (17), we obtain
(47)fr+1′′′+(fr-Aη2)fr+1′′-(2fr′+A)fr+1′+fr′′fr+1+Grθr+1+Gcϕr+1=frfr′′-fr′2,1Prθr+1′′+(fr-Aη2)θr+1′-(fr′+A)θr+1-θrfr+1′+θr′fr+1=frθr′-θrfr′,1Scϕr+1′′+(fr-Aη2)ϕr+1′-(fr′+A+γ)ϕr+1-ϕrfr+1′+ϕr′fr+1=frϕr′-ϕrfr′,(48)fr+1(0)=fw,fr+1′(0)=1,θr+1(0)=1,ϕr+1(0)=1,fr+1′(∞)=0,θr+1(∞)=0,ϕr+1(∞)=0.Applying the Chebyshev pseudospectral method on (47) we obtain the following SQLM scheme in a matrix form:
(49)[Δ1,1Δ1,2Δ1,3Δ2,1Δ2,2Δ2,3Δ3,1Δ3,2Δ3,3][Z1,r+1Z2,r+1Z3,r+1]=[Φ1,rΦ2,rΦ3,r]
subject to the boundary conditions
(50)fr+1(τN)=0,∑k=0N-‍DN-kfr+1(τk)=0,∑k=0N-‍D0kfr+1(τk)=1,θr+1(τN-)=1,θr+1(τ0)=0,ϕr+1(τN-)=1,ϕr+1(τ0)=0,
where
(51)Δ1,1=Λ1,1+Π1,1=D3-Aη2D2-(K+A)D+[fr′′]d-[2fr′]dD+[fr]dD2,Δ1,2=Λ1,2+Π1,2=GrI,Δ1,3=Λ1,3+Π1,3=GcI,Δ2,1=Λ2,1+Π2,1=[Θr′]d-[Θr]dD,Δ2,2=Λ2,2+Π2,2=1PrD2-Aη2D-AI+[fr]dD-[fr′]d,Δ2,3=Λ2,3+Π2,3=O1,Δ3,1=Λ3,1+Π3,1=[Φr′]d-[Φr]dD,Δ3,2=Λ3,2+Π3,2=O1,Δ3,3=Λ3,3+Π3,3=1ScD2-Aη2D-(A+γ)I+[fr]dD-[fr′]d,Φ1,r=frfr′′-fr′2,Φ2,r=frΘr′-fr′Θ,Φ3,r=frΦr′-fr′Φ,Z1,r=fr,Z2,r=Θr,Z3,r=Φr,
and []d denotes diag(); that is, the vector elements are placed on the main diagonal of a matrix whose entries everywhere else are zero. Starting from the initial approximations (44), the approximate SQLM solutions for f(η),θ(η), and ϕ(η) are obtained by solving (49).


## 3.1. Blasius Boundary Layer Equation

The Blasius boundary layer flow equation, in a dimensionless form, is expressed in terms of the similarity variableη∈[0,∞) and function f(η) as
(21)f′′′(η)+12f(η)f′′(η)=0,f(0)=0,f′(0)=0,f′(∞)=1,
where the primes denote differentiation with respect to η. To apply the proposed SLLM iterative method on (21), it is convenient to reduce the order of the governing equations by one through setting f′(η)=g(η). This results in the following system:
(22)f′=g,f(0)=0,(23)g′′+12f(η)g′(η)=0,g(0)=0,g(∞)=1.We observe that (23) is essentially a differential equation whose unknown is g(η), but the equation requires f(η) which can be considered to be a known input function. Thus, assuming that f(η) is known at a particular iteration level, the local linearization iteration scheme can be expressed by the following iteration formula;
(24)gr+1′′+12frgr+1′=0,gr+1(0)=0,gr+1(∞)=1,fr+1′=gr+1,fr+1(0)=0.Applying the Chebyshev pseudospectral method on (24), we obtain the following decoupled system of equations:
(25)A1gr+1=B1,gr+1(τN-)=0,gr+1(τ0)=1,Dfr+1=gr+1,fr+1(τN-)=0,
where
(26)A1=D2+12diag[fr]D,B1=O1,
where O1 is an (N-+1)×1 vector on which the boundary conditions are imposed and f and g are the values of functions f and g, respectively, when evaluated at the collocation points. Equation (25) constitutes the SLLM iteration scheme.The initial approximation required to start the iteration process can be chosen as a function that satisfies the governing boundary conditions and from known physical considerations of the flow properties. For most boundary layer flow problems it is well known that the flow velocity tends to the mainstream flow exponentially. Thus, a suitable initial approximation for the boundary layer problem (21) is
(27)f0(η)=η-1+e-η,g0(η)=1-e-η.To obtain the SQLM iterative scheme for the Blasius problem (21), we set
(28)L=f′′′,N=12ff′′,[f,f′,f′′,f′′′]=[z,z(1),z(2),z(3)].
Thus, using (17), we obtain
(29)fr+1′′′+12frfr+1′′+12fr′′fr+1=12fr′′fr,
subject to
(30)fr+1(0)=0,fr+1′(0)=0,fr+1′(∞)=1.Applying the spectral method, with derivative matrices on the iteration scheme (29) and the corresponding boundary conditions, gives the following matrix system:
(31)ArFr+1=Φr,
with the boundary conditions
(32)fr+1(τN)=0,∑k=0N-‍DN-kfr+1(τk)=0,∑k=0N-‍D0kfr+1(τk)=1,
where
(33)Ar=D3+12diag(fr)D2+12diag(fr′′),Φr=12frfr′′.


## 3.2. Unsteady Free Convective Heat and Mass Transfer on a Stretching Surface in a Porous Medium with Suction/Injection

In this section we consider a three-equation system that models the problem of unsteady free convective heat and mass transfer on a stretching surface in a porous medium in the presence of a chemical reaction. The governing equations [29, 30] for this problem are given as the following dimensionless system of equations and boundary conditions:
(34)f′′′+ff′′-(f′)2-Kf′-A(f′+η2f′′)+Grθ+Gcϕ=0,1Prθ′′-f′θ+fθ′-A(θ+12ηθ′)=0,1Scϕ′′-f′ϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0,(35)f(0)=fw,f′(0)=1,θ(0)=1,ϕ(0)=1,f′(∞)=0,θ(∞)=0,ϕ(∞)=0,
where f(η), θ(η), and ϕ(η) are, respectively, the dimensionless velocity, temperature, and concentration, fw is the suction/injection parameter, γ is the chemical reaction constant, Pr is the Prandtl number, Sc is the Schmidt number, K is the permeability parameter, and Gr and Gc are the temperature and concentration dependent Grashof numbers, respectively.To apply the proposed SLLM iterative method on (34), we set f′(η)=g(η) to obtain
(36)f′=g,(37)g′′+fg′′-g2-Kg-A(g+η2g′)+Grθ+Gcϕ=0,(38)1Prθ′′-gθ+fθ′-A(θ+12ηθ′)=0,(39)1Scϕ′′-gϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0.We observe that the energy equation (38) requires only f and g to resolve the solution of θ. The same applies to the mass transfer equation (39) for ϕ. Equation (37), on the other hand, requires the solutions of f,  θ and ϕ. In developing the SLLM iteration scheme we solve (38) and (39) for θ, and ϕ first, then use their updated solutions in the (37) for g. Thus, assuming that f(η) and g(η) are known at a particular iteration level, the local linearization iteration algorithm can be expressed as
(40)1Prθr+1′′-grθr+1+frθr+1′-A(θr+1+η2θr+1′)=0,1Scϕr+1′′-grϕr+1+frϕr+1′-(A+γ)ϕr+1-Aη2ϕr+1′=0,gr+1′′+frgr+1′-(A+K)gr+1+ηA2gr+1′-2grgr+1+Grθr+1+Gcϕr+1=-gr2,fr+1′=gr+1,fr+1(0)=fw,
subject to
(41)gr+1(0)=1,gr+1(∞)=0,θr+1(0)=1,θr+1(∞)=0,ϕr+1(0)=1,ϕr+1(∞)=0.Applying the Chebyshev pseudospectral method on (40)-(41), we obtain the following decoupled matrix system of equations:
(42)A1Θr+1=B1,θr+1(τN-)=1,θr+1(τ0)=0,A2Φr+1=B2,ϕr+1(τN-)=1,ϕr+1(τ0)=0,A3gr+1=B3-GrΘr+1-GcΦr+1,gr+1(τN-)=1,gr+1(τ0)=0,Dfr+1=gr+1,fr+1(τN-)=fw,
where
(43)A1=1PrD2+diag[fr-Aη2]D-diag[gr]-AI,B1=O1,A2=1ScD2+diag[fr-Aη2]D-diag[gr]-(A+γ)I,B2=O1,A3=D2+diag[fr-Aη2]D-diag[2gr]-(K+A)I,B3=-gr2,
where I is an (N-+1)×(N-+1) identity matrix, O1 is an (N-+1)×1 vector, and Θ, Φ,  g,  and  f are the approximate values of θ, ϕ, g,  and  f evaluated at the collocation points. After incorporating the boundary conditions in the matrices in (42), the approximate solutions of θ,  ϕ,  g, and f at each iteration level can be obtained by solving (42) in the order in which they are listed. A suitable initial approximation for the SLLM scheme is
(44)f0(η)=fw+1-e-η,θ0(η)=e-η,ϕ0(η)=e-η.To obtain the SQLM iterative scheme for the system (34), we set
(45)L1=f′′′-A(f′+η2f′′)+Grθ+Gcϕ,N1=ff′′-(f′)2,L2=1Prθ′′-A(θ+12ηθ′),N2=-f′θ+fθ′,L3=1Scϕ′′-A(ϕ+12ηϕ′)-γϕ,N3=-f′ϕ+fϕ′.To use formula (17) we define the following variables:
(46)[f,f′,f′′,f′′′]=[z1,z1(1),z1(2),z1(3)],[θ,θ′,θ′′]=[z2,z2(1),z2(2)],[ϕ,ϕ′,ϕ′′]=[z3,z3(1),z3(2)].Thus, using (45)-(46) in (17), we obtain
(47)fr+1′′′+(fr-Aη2)fr+1′′-(2fr′+A)fr+1′+fr′′fr+1+Grθr+1+Gcϕr+1=frfr′′-fr′2,1Prθr+1′′+(fr-Aη2)θr+1′-(fr′+A)θr+1-θrfr+1′+θr′fr+1=frθr′-θrfr′,1Scϕr+1′′+(fr-Aη2)ϕr+1′-(fr′+A+γ)ϕr+1-ϕrfr+1′+ϕr′fr+1=frϕr′-ϕrfr′,(48)fr+1(0)=fw,fr+1′(0)=1,θr+1(0)=1,ϕr+1(0)=1,fr+1′(∞)=0,θr+1(∞)=0,ϕr+1(∞)=0.Applying the Chebyshev pseudospectral method on (47) we obtain the following SQLM scheme in a matrix form:
(49)[Δ1,1Δ1,2Δ1,3Δ2,1Δ2,2Δ2,3Δ3,1Δ3,2Δ3,3][Z1,r+1Z2,r+1Z3,r+1]=[Φ1,rΦ2,rΦ3,r]
subject to the boundary conditions
(50)fr+1(τN)=0,∑k=0N-‍DN-kfr+1(τk)=0,∑k=0N-‍D0kfr+1(τk)=1,θr+1(τN-)=1,θr+1(τ0)=0,ϕr+1(τN-)=1,ϕr+1(τ0)=0,
where
(51)Δ1,1=Λ1,1+Π1,1=D3-Aη2D2-(K+A)D+[fr′′]d-[2fr′]dD+[fr]dD2,Δ1,2=Λ1,2+Π1,2=GrI,Δ1,3=Λ1,3+Π1,3=GcI,Δ2,1=Λ2,1+Π2,1=[Θr′]d-[Θr]dD,Δ2,2=Λ2,2+Π2,2=1PrD2-Aη2D-AI+[fr]dD-[fr′]d,Δ2,3=Λ2,3+Π2,3=O1,Δ3,1=Λ3,1+Π3,1=[Φr′]d-[Φr]dD,Δ3,2=Λ3,2+Π3,2=O1,Δ3,3=Λ3,3+Π3,3=1ScD2-Aη2D-(A+γ)I+[fr]dD-[fr′]d,Φ1,r=frfr′′-fr′2,Φ2,r=frΘr′-fr′Θ,Φ3,r=frΦr′-fr′Φ,Z1,r=fr,Z2,r=Θr,Z3,r=Φr,
and []d denotes diag(); that is, the vector elements are placed on the main diagonal of a matrix whose entries everywhere else are zero. Starting from the initial approximations (44), the approximate SQLM solutions for f(η),θ(η), and ϕ(η) are obtained by solving (49).


## 4. Numerical Convergence, Error, and Stability Analysis of the Iteration Schemes

The convergence and stability of the iteration schemes can be assessed by considering the norm of the difference in the values of functions between two successive iterations. Thus, for each iteration scheme, we define the following maximum error (Ed) at the (r+1)th iteration:
(52)Ed=Max(∥z1,r+1-z1,r∥∞,∥z2,r+1-z2,r∥∞,…,∥zm,r+1-zm,r∥∞),
where zi;  i=1,…,m are the governing unknown functions in the nonlinear system. If the iteration scheme converges, the error Ed is expected to decrease with an increase in the number of iterations. In this paper, the unknowns were calculated, for a given number of collocation points N-, until the following criteria for convergence was fulfilled at iteration r:
(53)Ed≤ε,
where ε is the convergence tolerance level. In this study, the convergence tolerance is set to be ε=10-9. The effect of the number of collocation points N- was examined in order to select the smallest value of N- which gives a consistent solution to the ε error level. This is achieved by repeatedly solving the governing equations using the proposed iteration schemes with different values of N- until the consistent solution is reached.


## 5. Improving the Convergence of the Spectral Local Linearization Method (SLLM)

In numerical linear algebra, successive overrelaxation (SOR) methods are used to accelerate the convergence rates of the Gauss-Seidel method in the solution of the linear systems of equation. In this section we propose a similar method to improve the convergence of the spectral relaxation method (SLLM) described in the previous section. If the SLLM scheme for solving the functionX at the (r+1)th iteration is
(54)AXr+1=B,
then the modified version of the SLLM is defined as
(55)AXr+1=(1-ω)AXr+ωB,
where A, B are matrices and ω is the convergence controlling relaxation parameter. We observe that when ω=1 (55) reduces to the original SLLM method. The results in the next section will show that, for certain values of ω, solving the modified SLLM results in improved efficiency and accuracy.


## 6. Results and Discussion

We solved the governing systems of (21) and (34) we solved using the spectral local linearization method (SLLM) and spectral quasilinearization method (SQLM) as described in the previous section. In this section, we present the results of the numerical computations for the velocity, temperature, and concentration profiles for various input parameters. Values of other flow properties such as skin friction f′′(0), surface heat transfer rate at the θ′(0), and mass transfer rate at the wall ϕ′(0) are also discussed in order to illustrate some special features of the solution and to compare the accuracy, convergence, and efficiency of the SQLM and SLLM algorithms proposed in this study. The accuracy of the present results was verified by comparing them with other results from literature which have been reported to be accurate to within a certain number of decimal digits. Further validation of the solution was established by comparing the present results with the numerical solutions obtained using the MATLAB inbuilt routine bvp4c which is a finite difference code that implements the three-stage Lobatto III formula [33]. In the calculations presented here, the values of the governing physical parameters were chosen deliberately to match some of the results from published literature in order to enable effective comparison.In Figure1 we give a comparison of the SLLM and SQLM for the plot of the Logarithm of the error Ed against the number of iterations in the Blasius boundary layer problem. It can be seen from the figure that the spectral quasi-linearization method converges much faster than the spectral local linearization method. The fast convergence rate of the quasilinearization method was established by many researchers including [9, 10, 34–36] who extended the application of the quasilinearization method to a wide variety of nonlinear boundary value problems and established that the method converges quadratically. The convergence rate of the proposed SLLM is seen to be linear in this case. The effect of introducing successive overrelaxation with ω=1.2 has the beneficial effect of increasing the convergence speed which, in this case, can be characterized by the steeper line. To obtain the optimal value of ω that yields the best accuracy, the results obtained using several values of ω>1 and ω<1 near ω=1 were compared. It was found that underrelaxation (ω<1) slows down the convergence, whereas overrelaxation with 1<ω≤1.2 improves the convergence of the SLLM iteration scheme. The best results for the Blasius boundary layer problem were obtained with ω≈1.2. We remark that the optimal ω was obtained through numerical experimentation.Figure 1
Comparison of Logarithm of errorEd for the Blasius boundary layer computed using the SLLM without SOR, SLLM with SOR (ω=1.2), and SQLM results.To further illustrate how the solution converges with an increase in iterations, we give the SLLM computed results for the skin frictionf′′(0) in Tables 1 and 2 (with SOR). The results are compared with the recently reported sixteen decimal digit accurate results of [31, 32]. Ganapol [32] used an algorithm based on a Maclaurin series with a Wynn-epsilon convergence acceleration and analytical continuation to obtain highly accurate skin friction coefficients for the Blasius layer flows. It can be seen from Table 1 that full convergence to 9-digit accurate results is achieved after 17 iterations. The 17-digit accurate results of [31] are matched exactly after 31 iterations. We remark that this value was reached after 52 iterations as reported in the work of  [32]. This shows that the basic SLLM is more efficient than the method used in [32].Table 1
Comparison between the basic SLLM computed values off′′(0) for the Blasius flow against the results of [31, 32].
Iter.
1
0.403512726
2
0.350208935
3
0.337194863
4
0.333547347
5
0.332492283
6
0.332184532
7
0.332094553
8
0.332068227
9
0.332060524
10
0.332058269
11
0.332057609
12
0.332057416
13
0.332057360
14
0.332057343
15
0.332057338
16
0.332057337
17
0.332057336
18
0.332057336
20
0.33205733621949529
25
0.33205733621520553
26
0.33205733621519900
28
0.33205733621519653
30
0.33205733621519632
31
0.33205733621519630
32
0.33205733621519630
⋮
0.33205733621519630  (see [31])
(52)
0.3320573362151965 (see [32])Table 2
Comparison between the SOR improved SLLM (ω=1.2) values of f′′(0) for the Blasius flow against the results of [31, 32].
1
0.284215262
2
0.350088883
3
0.329960091
4
0.332654418
5
0.331973067
6
0.332078099
7
0.332054014
8
0.332058085
9
0.332057206
10
0.332057364
11
0.332057331
12
0.332057337
13
0.332057336
14
0.332057336
15
0.332057336
16
0.332057336
17
0.332057336
18
0.332057336
19
0.33205733621518464
20
0.33205733621519862
21
0.33205733621519584
22
0.33205733621519639
23
0.33205733621519628
24
0.33205733621519630
25
0.33205733621519630
⋮
0.33205733621519630 (see [31])
(52)
0.3320573362151965 (see [32])Table2 gives the SOR-improved SLLM results (with ω=1.2) for the skin friction at different iterations. The table indicates that 9-digit accurate results are obtained after only 13 iterations and 17-digit accurate results are reached after 24 iterations. This clearly shows that using the SOR improves the convergence and efficiency of the basic SLLM.Table3 gives the SQLM results for the skin friction of the Blasius flow. It can be seen that the SQLM converges very rapidly compared to the SLLM in this example. Full convergence to the 17-digit accurate results reported in [31] is achieved after only 5 iterations. We remark that the results presented in Tables 1–3 were generated using N-=120 and η∞=16. Thus, under the same conditions, the SQLM is more efficient than the SLLM in this case.Table 3
SQLM computed results for the skin frictionf′′(0) of the Blasius flow.
Iter.
f
′
′
(
0
)
1
0.36124527510805664
2
0.33293906079206191
3
0.33205878995514977
4
0.33205733621994973
5
0.33205733621519630
6
0.33205733621519630
⋮
0.33205733621519630 (see [31])
(52)
0.3320573362151965 (see [32])The SLLM Log-error plots for the problem of unsteady free convective heat and mass transfer on a stretching surface are given in Figure2 for varying input parameters. It can be seen that the Logarithm of the error strictly decreases with an increase in the number of iterations when all the parameters are varied. This demonstrates the convergence and stability of the SLLM in solving the three-equation coupled system (34) for varying physical constants. It can also be seen that the convergence speed increases with an increase in A,  fw, and K but decreases with an increase in Gr.Effect of varying parameters on the Logarithm of the SLLM error.
(a)
(b)
(c)
(d)In Table4 we give a comparison of the values of f′′(0) and θ′(0) computed using the SLLM and SQLM approaches at fixed values of the governing input parameters. The table shows the effects of the number of the collocation points N- on the convergence with respect to the benchmark bvp4c numerical results. We remark that the stopping criteria used in both methods are given by (52) where the running of the computations is terminated when the maximum norm of the difference between the values of the governing functions at two successive iterations is less than the prescribed tolerance level (set to be 10-9 in this work). For the selected values of N-, the SLLM converged after 8 iterations. Increasing N- improved the accuracy, and the bvp4c results were matched exactly when N-≥50 for both f′′(0) and θ′(0). The SQLM algorithm was found converge after only 5 iterations when N≤100, but the results did not convergent to the bvp4c results. Increasing the number of collocation points slightly improved the accuracy but full convergence to eight decimal digits was not achieved. Further increase in N- led to the poor accuracy and convergence. This can be seen in the results corresponding to N-≥150 where convergence was not reached even after 100 iterations. These results indicate that the SQLM is not ideal for solving very large systems of equations. The poor accuracy of the SQLM is caused by the numerical difficulties associated with solving large matrix equations. The advantage of using the SLLM in solving large coupled systems of differential equations is that the SLLM algorithm leads to decoupled systems of sub-matrix equations which can easily be solved using standard direct methods in a sequential manner. With the SQLM the coupled system leads to a coupled iteration scheme which after discretization, using the spectral method, leads to large matrices. The dimension of the coefficient matrices and vectors to be solved in the SQLM algorithm is m(N-×1)×m(N-×1), where m is the number of unknowns in the governing systems. For example, in the case when N-=200, the dimension of the coefficient matrix used in the solution for f,θ, and ϕ is 603×603 in the SQLM, compared to 201×201 in the case of the SLLM. Furthermore, the displayed accuracy and stability of the SLLM results for both small and large values of N- suggest that the efficiency of SLLM is due more to the algorithm behind the iteration scheme than to the specific numerical method used to solve the underlying differential equations of the scheme. The SQLM accuracy and stability are strongly linked to the numerical scheme used to solve the differential equations of the scheme. The Chebyshev differentiation matrix is known to be very susceptible to very large round off errors and ill-conditioning for large numbers of collocation points N- (see e.g., [27, 37]), that is why the SQLM results in Table 4 seem to be less accurate when N- becomes very large. Table 4 clearly indicates that the SLLM is more practical to use for large systems of equations than SQLM as it yields more accurate results using only few collocation points, and there is no loss of accuracy when the side of the mesh is refined (when N- becomes large).Table 4
Effect of varyingN- on f′′(0) and θ′(0) using the SLLM, when K=1, Gr=Gc=A=1, and fw=2, γ=0.
Iter.
N
-
f
′
′
(
0
)
θ
′
(
0
)
(SLLM)
8
30
−2.31812851
−2.62227419
8
40
−2.31812936
−2.62228212
8
50
−2.31812942
−2.62228134
8
100
−2.31812942
−2.62228134
8
150
−2.31812942
−2.62228134
8
200
−2.31812942
−2.62228134
(SQLM)
5
30
−2.31691806
−2.62229791
5
40
−2.31816177
−2.62228177
5
50
−2.31813338
−2.62228131
5
100
−2.31812953
−2.62228134
100
150
−2.31811874
−2.62228134
100
200
−2.31821552
−2.62228134
bvp4c
−2.31812942
−2.62228134
[30]
−2.3181294
−2.6222813In Figure3 we give a comparison of the Log-error against iterations plot (Ed) between the SLLM and SQLM results for a fixed number of collocation points. From Figure 3 it can be seen that the SQLM rapidly converges but the accuracy does not improve after 5 iterations. After 5 iterations the error plot seems to plateau at a certain level and remains at that level even when considering more iterations. The SLLM results indicate strict convergence with the error progressively reduced with an increase in the number of iterations. This explains why the SLLM gives more accurate results than the SQLM for larger systems as was observed in Table 4.Figure 3
Comparison of the Logarithm of the errorEd for the unsteady boundary layer problem between the SLLM and SQLM results when N-=150.In Table5 we give a comparison between the basic SLLM results and the results obtained using the SOR accelerated version of the SLLM for the skin friction f′′(0) and wall heat transfer rate θ′(0). It can be noted from Table 3 that introducing the convergence accelerating relaxation parameter ω leads to the improvement of the basic SLLM for some cases. It can be seen that convergence improves with an increase in fw,  A, and K and slows down with an increase in Gr. This result is consistent with the observation made in Figure 2. In cases where convergence is slow, underrelaxation with (0.9≤ω<1) accelerates the convergence. We also observe that the effect of increasing the values of Gr results in an increase in the skin friction and in the absolute values of the wall heat transfer rate. The wall heat transfer rate decreases in absolute value as K increases. An increase in the unsteadiness parameter A and the suction parameter fw results in an increase in the skin friction and the wall heat transfer rate. We remark that the trends in the effect of the mass transfer rate at the wall were found to be similar to the effects observed in the surface heat transfer rate for all varied input parameters. The trends observed in Table 5 are consistent with the observations made in related studies [29, 30].Table 5
Effect of varyingA, fw, Gr, and K on f′′(0) and θ′(0), when Pr=Sc=γ=Gc=1.
A
f
w
Gr
K
Basic SLLM
SLLM with SOR (ω)
Iter.
f
′
′
(
0
)
θ
′
(
0
)
ω
Iter.
f
′
′
(
0
)
θ
′
(
0
)
0
1
1
1
11
−1.31052731
−1.65084042
0.95
7
−1.31052731
−1.65084042
1
1
1
1
10
−1.59363827
−1.89854697
0.95
7
−1.59363827
−1.89854697
2
1
1
1
8
−1.85550919
−2.12936888
0.95
7
−1.85550919
−2.12936888
5
1
1
1
6
−2.49344366
−2.70569628
1.00
6
−2.49344366
−2.70569628
10
1
1
1
6
−3.27938346
−3.44160771
1.00
6
−3.27938346
−3.44160771
1
0
1
1
10
−1.03558899
−1.34595467
0.95
7
−1.03558899
−1.34595467
1
2
1
1
8
−2.34393108
−2.61971104
0.95
6
−2.34393108
−2.61971104
1
4
1
1
6
−4.15359239
−4.35571036
1.00
6
−4.15359239
−4.35571036
1
6
1
1
6
−6.09385439
−6.24419885
1.00
6
−6.09385439
−6.24419885
1
10
1
1
5
−10.05256264
−10.14874633
1.00
5
−10.05256264
−10.14874633
1
1
0
1
8
−1.91528421
−1.87004601
1.00
8
−1.91528421
−1.87004601
1
1
1
1
10
−1.59363827
−1.89854697
0.95
7
−1.59363827
−1.89854697
1
1
5
1
13
−0.41099838
−1.98803325
0.90
9
−0.41099838
−1.98803325
1
1
6
1
15
−0.13374591
−2.00648856
0.90
9
−0.13374591
−2.00648856
1
1
10
1
19
0.92495170
−2.07084699
0.90
10
0.92495170
−2.07084699
1
1
1
0
12
−1.21175687
−1.93465038
0.95
8
−1.21175687
−1.93465038
1
1
1
2
9
−1.90764681
−1.87149682
0.95
7
−1.90764681
−1.87149682
1
1
1
4
8
−2.42217698
−1.83197529
0.95
7
−2.42217698
−1.83197529
1
1
1
6
7
−2.84650750
−1.80348579
0.95
7
−2.84650750
−1.80348579
1
1
1
10
7
−3.54376838
−1.76361252
0.95
7
−3.54376838
−1.76361252Figures4, 5, 6, and 7 illustrate the effects of the input parameters K,  A,  fw and Gr, respectively, on the temperature and concentration profiles. It is noted that the velocity decreases with K, A, and fw but increases with the buoyancy parameter Gr. The temperature decreases with A,  fw, and Gr but increases with an increase in K. Again, the observed trends are consistent with the observations reported in [29, 30]. The effect of the varied parameters on the concentration was found to be similar to the trends observed in the temperature profiles.Effect ofK on the velocity and temperature profiles.
(a)
(b)Effect ofA on the velocity and temperature profiles.
(a)
(b)Effect offw on the velocity and temperature profiles.
(a)
(b)Effect ofGr on the velocity and temperature profiles.
(a)
(b)


## 7. Conclusion

In this work we have introduced a new method for solving systems of nonlinear boundary value equations. The proposed method, called spectral local linearization method (SLLM), is based on a simple idea of the decoupling systems of equations using the linearization of the unknown functions in a sequential manner according to the order of the listing of the governing equations of the system. The linearized equations, which now form a sequence of linear differential equations with variable coefficients, were solved using the Chebyshev spectral collocation method. The applicability of the method was examined on similarity boundary layer problems with exponentially decaying profiles. The Blasius boundary layer problem and a three-equation coupled system that models the problem of unsteady, free convective heat and mass transfer were considered as test cases to examine the applicability of the proposed SLLM. The accuracy of the SLLM was confirmed against known accurate results reported in the published literature and against numerical results generated using Matlab’sbvp4c routine for solving boundary value problems. The performance of the SLLM was measured against the spectral quasi-linearization method (SQLM) which uses the Chebyshev collocation method to discretize the quasilinearization method (a Newton-Raphson based iteration method for solving boundary value problem).From this preliminary study, the following observations were made about the proposed SLLM.(1)
The algorithm of the SLLM is very easy to develop and implement, as it is based on a simple univariate linearization of nonlinear functions.(2)
The method is numerically efficient, since it results in a series of equations which are solved in a sequential manner by reusing the information from the solution of one equation in the next equation.(3)
The method was found to converge rapidly to the expected solutions for all the input parameters considered in this study.(4)
The convergence speed of the method can readily be improved by using successive overrelaxation (SRM) techniques. In the case of the Blasius boundary layer, overrelaxation improved convergence, and in the unsteady heat and mass transfer problem, the convergence was improved by underrelaxation.(5)
The convergence speed of the SLLM improves with an increase in the values of the parametersA,fw,  and  K and slows down with an increase in Gr.(6)
The method is more computationally efficient than the SQLM in the sense that it gives accurate results even with a few collocation points (smallN-). The accuracy of the SQLM was seen to deteriorate when the number of collocation points became very large. This is due to the numerical difficulties associated with solving large matrix systems. Unlike the SQLM, the SLLM does not suffer from the loss of accuracy when the number of collocation points is large.Because of its demonstrated accuracy, efficiency, and simplicity, it is envisaged that the proposed SLLM method could be used as a viable method for solving some classes of similarity variable boundary layer problems.



---

*Source: 423628-2013-09-17.xml*