# On the Comparison between Compact Finite Difference and Pseudospectral Approaches for Solving Similarity Boundary Layer Problems

Authors: P. G. Dlamini; S. S. Motsa; M. Khumalo
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We introduce two methods based on higher order compact finite
differences for solving boundary layer problems. The methods called compact
finite difference relaxation method (CFD-RM) and compact finite difference
quasilinearization method (CFD-QLM) are an alternative form of the
spectral relaxation method (SRM) and spectral quasilinearization method
(SQLM). The SRM and SQLM are Chebyshev pseudospectral-based methods
which have been successfully used to solve boundary layer problems. The
main objective of this paper is to give a comparison of the compact finite difference
approach against the pseudo-spectral approach in solving similarity
boundary layer problems. In particular, we seek to identify the most accurate
and computationally efficient method for solving systems of boundary
layer equations in fluid mechanics. The results of the two approaches are
comparable in terms of accuracy for small systems of equations. For larger
systems of equations, the proposed compact finite difference approaches are
more accurate than the spectral-method-based approaches.


---

## Body

## 1. Introduction

Low-order (second or lower) finite difference schemes are not accurate enough for solving many problems in fluid dynamics and other areas. Recently the focus has shifted to higher-order compact finite difference (CFD) schemes. Researchers have found significant improvement to the accuracy of numerical solutions by using fourth or sixth-order CFD schemes compared to the convectional second order finite central difference scheme [1]. Various CFD schemes used for applications such as interpolation, filtering, and evaluating high-order derivatives were discussed in detail by Lele [2]. CFD schemes have largely been applied to solve partial differential equations, for example, Burger’s equation [3, 4], Navier-Stokes equation [5], Korteweg-de Vries equation [6], Black-Scholes equation [7], and many more [8–10]. A limited number of researchers have utilised the CFD schemes in ordinary differential equations. For example in [11], Zhao solved two-point boundary value problems. In [12], Zhao and Corless used the CFD schemes to solve integro-differential equations.The advantage of the higher-order CFD schemes is that they give high accuracy on coarser grids with greater computational efficiency [13]. The difficulty that comes with the higher order CFD schemes especially near the boundaries makes many researchers shy away from using them in practical computations. To retain the accuracy of the schemes at the boundaries, the schemes are adjusted for the boundary points.When compared to spectral methods, compact schemes are more flexible in terms of application to complex geometries and boundary conditions. Lele [2] pointed out that the use of spectral methods in turbulent fluid flows is limited to flows in simple domains and simple boundary conditions. Spectral methods become less accurate for problems with complex geometries. Rai [14] gave a comprehensive comparison between results obtained using finite differences and spectral methods for direct simulation of turbulent flows. They used high-order accurate upwind schemes. They concluded that the spectral method is extremely accurate but it has restrictions on the type of geometry and grids that can be efficiently handled. They further state that unlike the spectral methods, the finite difference method can be used efficiently with curvilinear grids.In recent years Motsa and his coworkers [15–18] have developed successful methods based on the spectral method to solve nonlinear similarity boundary layer problems. The methods include, among others, the spectral relaxation method (SRM) [15, 19, 20], the spectral successive linearisation method [21–23], spectral homotopy analysis method [24–27], and the spectral quasilinearisation method (SQLM) [16]. The SRM is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the Chebyshev spectral collocation method. The SQLM combines the quasilinearization technique developed by Bellman and Kalaba [28] to linearize nonlinear differential equations and solve the resulting equations using the spectral method. In this work we present alternative approach to the implementation of the SRM and SQLM. Instead of using the spectral method in these methods, we use the higher order CFD schemes, and hence the resulting methods are the compact finite difference relaxation method (CFD-RM) and compact finite difference quasilinearization method (CFD-QLM), respectively.The main objective of this work is to compare the spectral-method-based and the CFD-based methods discussed above. We compare the performance of the methods in terms of accuracy and computational speed when solving nonlinear boundary layer problems in one dimension and three dimensions. We first consider the flow of a viscous incompressible electrically conducting fluid over a continuously shrinking sheet which is governed by a third-order nonlinear differential equation. The known exact analytical solution [29, 30] of this problem is used as a benchmark to validate the accuracy of the proposed algorithm discussed in this work. We also consider a 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 [23, 31].


## 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 Relaxation Method

The spectral relaxation method (SRM) is a new method that has been introduced recently by Motsa et al. [15, 19, 20] to solve initial and boundary value problems. The method is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the Chebyshev spectral collocation method. The algorithm for the method is summarized as follows.(1)
Arrange the governing nonlinear equations in a particular order, placing the equations with the least unknowns at the top of the equations list.(2)
Assign the labelsZ1,Z2,Z3,… to the ordered equations obtained in the above step, where each Zi(i=1,2,3,…) is an unknown function which, in the ith equation, is identified as the unknown function associated with the highest-order derivative.(3)
In the equation forZ1 (1st equation), the iteration scheme is developed by assuming that only linear terms in Z1 are to be evaluated at the current iteration level (denoted by r+1) and all other terms (linear and nonlinear) in Z2,Z3,… are assumed to be known from the previous iteration (denoted by r). In addition nonlinear terms in Z1 are also evaluated at the previous iteration. Furthermore, all derivative terms in Z1 are assumed to be known from the previous iteration.(4)
Similarly, in developing the iteration scheme in the equation forZ2 (2nd equation), only linear terms in Z2 are evaluated at the current iteration level (r+1) with all other terms evaluated at the previous level, except Z1 which is now known from the solution of the first equation.(5)
This process is repeated in theith equation (i=3,4,…) using the updated solutions for Zi-1 obtained from the previous i-1 equations.The resulting iteration scheme is integrated using the Chebyshev spectral method. The region of integration is discretized using the Gauss-Lobatto points defined by
(1)τj=cos(πjN),j=1,2,3,…,N,
where N is the number of collocation points used. The Chebyshev spectral collocation method is based on the idea of introducing a differentiation matrix D which is used to approximate the derivatives of the unknown variables Zi(x)(i=1,2,3,…) at the collocation points as the matrix vector product
(2)dZidx=∑k=0N‍DjkZi=DZi,j=1,2,…,N,
where Zi=[Zi(τ0),Zi(τ1),…,Zi(τN)] is the vector function at the collocation points τj.


### 2.2. Compact Finite Difference Relaxation Method

The compact finite difference relaxation method (CFD-RM) uses the same procedure followed in the SRM. The difference is that instead of using spectral methods to solve the resulting iteration schemes, higher order compact finite difference schemes are used. In this work we use sixth-order CFD schemes.In the derivation of the CFD schemes we consider a one-dimensional uniform mesh on the region[a,b] with nodes xi (i=1,2,…,N) where
(3)a=x1<x2<⋯<xN=b
and a corresponding function yi=y(xi) at the nodes. The distance between any two successive nodes is a constant h=xi-xi-1. Sixth-order approximations of the first, second, and third derivatives at interior nodes can be obtained using the following schemes (see [2] for details):
(4)13yi-1′+yi′+13yi+1′=149yi+1-yi-12h+19yi+2-yi-24h,(5)211yi-1′′+yi′′+211yi+1′′=1211yi+1-2yi+yi-1h2+311yi+2-2yi+yi-24h2,(6)716yi-1′′′+yi′′′+716yi+1′′′=2yi+2-2yi+1-2yi-1-yi-22h3-18yi+3-3yi+1-3yi-1-yi-38h3.For illustrative purposes we describe the application of the CFD schemes to second-order differential equations fory(x) with known boundary conditions at y(a) and y(b). Consider the nonlinear differential equations
(7)y′′+p(x)y′+q(x)y(x)+f(x,y,y′)=0,y(a)=ya,y(b)=yb,
where f(x,y,y′) is a nonlinear function, p(x) and q(x) are known functions of x, and ya and yb are known constants. In solving (7), we apply the CFD approximation for the first and second derivatives given by (4) and (5), respectively, at the interior nodes (i=2,…,N-1). Since we know boundary conditions at i=1 and i=N, the CFD schemes must be adjusted for the nodes near the boundary points. In order to maintain the order O(h6) accuracy at the boundary points as in the interior points and to maintain the same tridiagonal format, we use the following one-sided scheme at i=2:
(8)y2′+13y3′=1h(a1y1+a2y2+a3y3+a4y4+a5y5+a6y6+a7y7),
and when i=N-1, we use
(9)13yN-2′+yN-1′=1h(b1yN+b2yN-1+b3yN-3+b4yN-4+b5yN-5+b6yN-6+b7yN-7),
where ai,bi(i=1,…,7) are constants to be determined. To obtain a sixth-order accurate scheme, we use Taylor series expansion about x1 and xN up to O(h7), in (8) and (9), respectively, and equate terms of order h. In each case, we obtain a system of seven linear algebraic equations in seven unknowns which are solved to give
(10)a1=-745,a2=-1712,a3=8336,a4=-119,a5=23,a6=-37180,a7=136,bk=-ak,k=1,2,…,7.Similarly, for the second derivatives, we use(11)y2′′+211y3′′=1h2(c1y1+c2y2+c3y3+c4y4+c5y5+c6y6+c7y7+c8y8),
at i=2 and
(12)211yN-2′′+yN-1′′=1h2(d1yN+d2yN-1+d3yN-3+d4yN-4+d5yN-5+d6yN-6+d7yN-7+d8yN-8)
at i=N-1. The parameters ci,di can be determined by expanding (11) and (12) using Taylor series and equating powers of h and subsequently solving the resulting equations. This gives
(13)c1=3145,c2=-19110,c3=-339110,c4=1933396,c5=-4011,c6=9655,c7=-479990,c8=13220,di=ci,i=1,2,…,8.Using the above equations, the equations for approximating the first- and second-order derivatives can be expressed as(14)A2Y′′=B2Y+K2,A1Y′=B1Y+K1,
where
(15)A1=[1131311313113⋱⋱⋱13113131](N-2)×(N-2),K1=1h[a1y1-y1360⋮⋮0yN36-a1yN](N-2)×1,(16)B1=1h[a2a3a4a5a6a7-79079136-136-79079136⋱⋱⋱⋱⋱-136-79079136-136-79079-a7-a6-a5-a4-a3-a2](N-2)×(N-2),(17)A2=[121121112112111211⋱⋱⋱21112112111](N-2)×(N-2),K2=1h2[c1y13y1440⋮⋮03yN44c1yN](N-2)×1,B2=1h2[c2c3c4c5c6c7c81211-512212113443441211-51221211344⋱⋱⋱⋱⋱3441211-512212113443441211-51221211c8c7c6c5c4c3c2](N-2)×(N-2)Y′=[y2′,y3′,…,yN-2′,yN-1′]T,Y′′=[y2′′,y3′′,…,yN-2′′,yN-1′′]T.For the purposes of the examples given in this work, we also illustrate the application of the CFD schemes for the following set of boundary conditions:(18)y(a)=ya,y′(a)=yda,y′(b)=ydb.We adjust the schemes near the boundary points as discussed above. In this case, the schemes ati=N-2 are also adjusted. At i=2, we use the following one-sided scheme:
(19)y2′+13y3′=a1y1′+1h(a2y1+a3y2+a4y3+a5y4+a6y5+a7y6).
At i=N-1, we use
(20)13yN-2′+yN-1′=b1yN′+1h(b2yN-1+b3yN-2+b4yN-3+b5yN-4+b6yN-5+b7yN-6),
and at i=N-2, we use
(21)13yN-3′+yN-2′+13yN-1′=c1yN′+1h(c2yN-1+c3yN-2+c4yN-3+c5yN-4+c6yN-5+c7yN-6)
with the constants given by
(22)a1=-16,a2=-203360,a3=-512,a4=1918,a5=-19,a6=124,a7=-180,b1=463,b2=15184,b3=-701252,b4=311189,b5=-1921,b6=71252,b7=-29756,c1=5441,c2=373441,c3=-25294,c4=-9291323,c5=-31441,c6=-2147,c7=-52646.Similarly for the second derivative, we use(23)y2′′+211y3′′=d1y1′h+1h2(d2y1+d3y2+d4y3+d5y4+d6y5+d7y6+d8y7),
at i=2. At i=N-1, we use
(24)211yN-2′′+yN-1′′=e1yN′h+1h2(e2yN-1+e3yN-2+e4yN-3+e5yN-4+e6yN-5+e7yN-6+e8yN-7),
and at N-3,
(25)211yN-3′′+yN-2′′+211yN-1′′=f1yN′h+1h2(f2yN-1+f3yN-2+f4yN-3+f5yN-4+f6yN-5+f7yN-6+f8yN-7).The constants are found to be(26)d1=91220,d2=63413600,d3=-13544,d4=11188,d5=118,d6=-3176,d7=91100,d8=-1792,e1=1759124352,e2=5511732670,e3=-6394110890,e4=31288739204,e5=-194753267,e6=31161089,e7=-7780199010,e8=634165340,f1=19256088,f2=16971331,f3=-34531331,f4=55813993,f5=-4312662,f6=1471331,f7=-2457986,f8=51331.Combining the schemes for approximating the first and second derivatives at interior points with (19)–(23), the equations for approximating the first and second derivatives are given by
(27)A1Y′=B1Y+K1+K1′,A2Y′′=B2Y+K2+K2′,
where A1 and A2 are the same as in (14),
(28)B1=1h[a3a4a5a6a7-79079136-136-79079136⋱⋱⋱⋱⋱-136-79079136c7c6c5c4c3c1b7b6b5b4b3b2](N-2)×(N-2),K1=1h[a2y1-y1360⋮0](N-2)×1,K1′=[a1y1′0⋮c1yN′b1yN′](N-2)×1,B2=1h2[d3d4d5d6d7d81211-512212113443441211-51221211344⋱⋱⋱⋱⋱3441211-51221211344f8f7f6f5f4f3f2e8e7e6e5e4e3e2](N-2)×(N-2),K2=1h2[c2y13y1440⋮0](N-2)×1,K2′=1h[c1y1′0⋮f1yN′e1yN′](N-2)×1.Similarly, third derivatives can be expressed in the form(29)A3Y′′′=B3Y+K3+K3′.


### 2.3. Spectral Quasilinearization Method

In this section we give a brief description of the spectral quasilinearization method (SQLM). In this method we make use of the quasilinearization method (QLM) and spectral method. The QLM was initially proposed by Bellman and Kalaba [28] to solve nonlinear boundary layer problems. To develop the SQLM, we consider a system of m nonlinear differential equations in m unknowns zi(η)(i=1,2,…,m) where η is the independent variable. The system can be written as a sum of its linear L and nonlinear components N as
(30)L[z1(η),z2(η),…,zm(η)]+N[z1(η),z2(η),…,zm(η)]=H(η).
Define the vector Zi to be the vector of the derivatives of the variable zi with respect to the independent variable η; that is,
(31)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 (30) can be written as
(32)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=∑j=1m‍∑p=0ni‍αij[p]Zj(p)+Ni[Z1,Z2,…,Zm]=Hi,
where αi,j[p] are the constant coefficients of zj(p), the derivative of zj(j=1,2,…,m) that appears in the ith equation for i=1,2,…,m. Assume that the solution zi(η) of (32) 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 (32) can be linearised using one-term Taylor series for multiple variables, so that (32) can be approximated as
(33)Li[Z1,r+1,…,Zm,r+1]+Ni[⋯]+∑j=1m‍∑p=0ni‍(Zj,r+1(p)-Zj,r(p))∂Ni∂Zj(p)[⋯]=Hi,
where
(34)[⋯]=[Z1,r,Z2,r,…,Zm,r].Equation (33) can be rewritten as
(35)Li[Z1,r+1,…,Zm,r+1]+∑j=1m‍∑p=0ni‍Zj,r+1(p)∂Ni∂Zj(p)[⋯]=Hi+∑j=1m‍∑p=0ni‍Zj,r(p)∂Ni∂Zj(p)[⋯]-Ni[⋯].To solve the iteration scheme (35), we use the Chebyshev spectral collocation method.


### 2.4. Compact Finite Difference Quasilinearization Method

Using compact finite difference schemes to solve the QLM scheme (35) results in the compact finite difference quasilinearization method (CFD-QLM). We will use the sixth-order finite difference schemes the same way we used them for the CFD-RM.


## 2.1. Spectral Relaxation Method

The spectral relaxation method (SRM) is a new method that has been introduced recently by Motsa et al. [15, 19, 20] to solve initial and boundary value problems. The method is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the Chebyshev spectral collocation method. The algorithm for the method is summarized as follows.(1)
Arrange the governing nonlinear equations in a particular order, placing the equations with the least unknowns at the top of the equations list.(2)
Assign the labelsZ1,Z2,Z3,… to the ordered equations obtained in the above step, where each Zi(i=1,2,3,…) is an unknown function which, in the ith equation, is identified as the unknown function associated with the highest-order derivative.(3)
In the equation forZ1 (1st equation), the iteration scheme is developed by assuming that only linear terms in Z1 are to be evaluated at the current iteration level (denoted by r+1) and all other terms (linear and nonlinear) in Z2,Z3,… are assumed to be known from the previous iteration (denoted by r). In addition nonlinear terms in Z1 are also evaluated at the previous iteration. Furthermore, all derivative terms in Z1 are assumed to be known from the previous iteration.(4)
Similarly, in developing the iteration scheme in the equation forZ2 (2nd equation), only linear terms in Z2 are evaluated at the current iteration level (r+1) with all other terms evaluated at the previous level, except Z1 which is now known from the solution of the first equation.(5)
This process is repeated in theith equation (i=3,4,…) using the updated solutions for Zi-1 obtained from the previous i-1 equations.The resulting iteration scheme is integrated using the Chebyshev spectral method. The region of integration is discretized using the Gauss-Lobatto points defined by
(1)τj=cos(πjN),j=1,2,3,…,N,
where N is the number of collocation points used. The Chebyshev spectral collocation method is based on the idea of introducing a differentiation matrix D which is used to approximate the derivatives of the unknown variables Zi(x)(i=1,2,3,…) at the collocation points as the matrix vector product
(2)dZidx=∑k=0N‍DjkZi=DZi,j=1,2,…,N,
where Zi=[Zi(τ0),Zi(τ1),…,Zi(τN)] is the vector function at the collocation points τj.


## 2.2. Compact Finite Difference Relaxation Method

The compact finite difference relaxation method (CFD-RM) uses the same procedure followed in the SRM. The difference is that instead of using spectral methods to solve the resulting iteration schemes, higher order compact finite difference schemes are used. In this work we use sixth-order CFD schemes.In the derivation of the CFD schemes we consider a one-dimensional uniform mesh on the region[a,b] with nodes xi (i=1,2,…,N) where
(3)a=x1<x2<⋯<xN=b
and a corresponding function yi=y(xi) at the nodes. The distance between any two successive nodes is a constant h=xi-xi-1. Sixth-order approximations of the first, second, and third derivatives at interior nodes can be obtained using the following schemes (see [2] for details):
(4)13yi-1′+yi′+13yi+1′=149yi+1-yi-12h+19yi+2-yi-24h,(5)211yi-1′′+yi′′+211yi+1′′=1211yi+1-2yi+yi-1h2+311yi+2-2yi+yi-24h2,(6)716yi-1′′′+yi′′′+716yi+1′′′=2yi+2-2yi+1-2yi-1-yi-22h3-18yi+3-3yi+1-3yi-1-yi-38h3.For illustrative purposes we describe the application of the CFD schemes to second-order differential equations fory(x) with known boundary conditions at y(a) and y(b). Consider the nonlinear differential equations
(7)y′′+p(x)y′+q(x)y(x)+f(x,y,y′)=0,y(a)=ya,y(b)=yb,
where f(x,y,y′) is a nonlinear function, p(x) and q(x) are known functions of x, and ya and yb are known constants. In solving (7), we apply the CFD approximation for the first and second derivatives given by (4) and (5), respectively, at the interior nodes (i=2,…,N-1). Since we know boundary conditions at i=1 and i=N, the CFD schemes must be adjusted for the nodes near the boundary points. In order to maintain the order O(h6) accuracy at the boundary points as in the interior points and to maintain the same tridiagonal format, we use the following one-sided scheme at i=2:
(8)y2′+13y3′=1h(a1y1+a2y2+a3y3+a4y4+a5y5+a6y6+a7y7),
and when i=N-1, we use
(9)13yN-2′+yN-1′=1h(b1yN+b2yN-1+b3yN-3+b4yN-4+b5yN-5+b6yN-6+b7yN-7),
where ai,bi(i=1,…,7) are constants to be determined. To obtain a sixth-order accurate scheme, we use Taylor series expansion about x1 and xN up to O(h7), in (8) and (9), respectively, and equate terms of order h. In each case, we obtain a system of seven linear algebraic equations in seven unknowns which are solved to give
(10)a1=-745,a2=-1712,a3=8336,a4=-119,a5=23,a6=-37180,a7=136,bk=-ak,k=1,2,…,7.Similarly, for the second derivatives, we use(11)y2′′+211y3′′=1h2(c1y1+c2y2+c3y3+c4y4+c5y5+c6y6+c7y7+c8y8),
at i=2 and
(12)211yN-2′′+yN-1′′=1h2(d1yN+d2yN-1+d3yN-3+d4yN-4+d5yN-5+d6yN-6+d7yN-7+d8yN-8)
at i=N-1. The parameters ci,di can be determined by expanding (11) and (12) using Taylor series and equating powers of h and subsequently solving the resulting equations. This gives
(13)c1=3145,c2=-19110,c3=-339110,c4=1933396,c5=-4011,c6=9655,c7=-479990,c8=13220,di=ci,i=1,2,…,8.Using the above equations, the equations for approximating the first- and second-order derivatives can be expressed as(14)A2Y′′=B2Y+K2,A1Y′=B1Y+K1,
where
(15)A1=[1131311313113⋱⋱⋱13113131](N-2)×(N-2),K1=1h[a1y1-y1360⋮⋮0yN36-a1yN](N-2)×1,(16)B1=1h[a2a3a4a5a6a7-79079136-136-79079136⋱⋱⋱⋱⋱-136-79079136-136-79079-a7-a6-a5-a4-a3-a2](N-2)×(N-2),(17)A2=[121121112112111211⋱⋱⋱21112112111](N-2)×(N-2),K2=1h2[c1y13y1440⋮⋮03yN44c1yN](N-2)×1,B2=1h2[c2c3c4c5c6c7c81211-512212113443441211-51221211344⋱⋱⋱⋱⋱3441211-512212113443441211-51221211c8c7c6c5c4c3c2](N-2)×(N-2)Y′=[y2′,y3′,…,yN-2′,yN-1′]T,Y′′=[y2′′,y3′′,…,yN-2′′,yN-1′′]T.For the purposes of the examples given in this work, we also illustrate the application of the CFD schemes for the following set of boundary conditions:(18)y(a)=ya,y′(a)=yda,y′(b)=ydb.We adjust the schemes near the boundary points as discussed above. In this case, the schemes ati=N-2 are also adjusted. At i=2, we use the following one-sided scheme:
(19)y2′+13y3′=a1y1′+1h(a2y1+a3y2+a4y3+a5y4+a6y5+a7y6).
At i=N-1, we use
(20)13yN-2′+yN-1′=b1yN′+1h(b2yN-1+b3yN-2+b4yN-3+b5yN-4+b6yN-5+b7yN-6),
and at i=N-2, we use
(21)13yN-3′+yN-2′+13yN-1′=c1yN′+1h(c2yN-1+c3yN-2+c4yN-3+c5yN-4+c6yN-5+c7yN-6)
with the constants given by
(22)a1=-16,a2=-203360,a3=-512,a4=1918,a5=-19,a6=124,a7=-180,b1=463,b2=15184,b3=-701252,b4=311189,b5=-1921,b6=71252,b7=-29756,c1=5441,c2=373441,c3=-25294,c4=-9291323,c5=-31441,c6=-2147,c7=-52646.Similarly for the second derivative, we use(23)y2′′+211y3′′=d1y1′h+1h2(d2y1+d3y2+d4y3+d5y4+d6y5+d7y6+d8y7),
at i=2. At i=N-1, we use
(24)211yN-2′′+yN-1′′=e1yN′h+1h2(e2yN-1+e3yN-2+e4yN-3+e5yN-4+e6yN-5+e7yN-6+e8yN-7),
and at N-3,
(25)211yN-3′′+yN-2′′+211yN-1′′=f1yN′h+1h2(f2yN-1+f3yN-2+f4yN-3+f5yN-4+f6yN-5+f7yN-6+f8yN-7).The constants are found to be(26)d1=91220,d2=63413600,d3=-13544,d4=11188,d5=118,d6=-3176,d7=91100,d8=-1792,e1=1759124352,e2=5511732670,e3=-6394110890,e4=31288739204,e5=-194753267,e6=31161089,e7=-7780199010,e8=634165340,f1=19256088,f2=16971331,f3=-34531331,f4=55813993,f5=-4312662,f6=1471331,f7=-2457986,f8=51331.Combining the schemes for approximating the first and second derivatives at interior points with (19)–(23), the equations for approximating the first and second derivatives are given by
(27)A1Y′=B1Y+K1+K1′,A2Y′′=B2Y+K2+K2′,
where A1 and A2 are the same as in (14),
(28)B1=1h[a3a4a5a6a7-79079136-136-79079136⋱⋱⋱⋱⋱-136-79079136c7c6c5c4c3c1b7b6b5b4b3b2](N-2)×(N-2),K1=1h[a2y1-y1360⋮0](N-2)×1,K1′=[a1y1′0⋮c1yN′b1yN′](N-2)×1,B2=1h2[d3d4d5d6d7d81211-512212113443441211-51221211344⋱⋱⋱⋱⋱3441211-51221211344f8f7f6f5f4f3f2e8e7e6e5e4e3e2](N-2)×(N-2),K2=1h2[c2y13y1440⋮0](N-2)×1,K2′=1h[c1y1′0⋮f1yN′e1yN′](N-2)×1.Similarly, third derivatives can be expressed in the form(29)A3Y′′′=B3Y+K3+K3′.


## 2.3. Spectral Quasilinearization Method

In this section we give a brief description of the spectral quasilinearization method (SQLM). In this method we make use of the quasilinearization method (QLM) and spectral method. The QLM was initially proposed by Bellman and Kalaba [28] to solve nonlinear boundary layer problems. To develop the SQLM, we consider a system of m nonlinear differential equations in m unknowns zi(η)(i=1,2,…,m) where η is the independent variable. The system can be written as a sum of its linear L and nonlinear components N as
(30)L[z1(η),z2(η),…,zm(η)]+N[z1(η),z2(η),…,zm(η)]=H(η).
Define the vector Zi to be the vector of the derivatives of the variable zi with respect to the independent variable η; that is,
(31)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 (30) can be written as
(32)Li[Z1,Z2,…,Zm]+Ni[Z1,Z2,…,Zm]=∑j=1m‍∑p=0ni‍αij[p]Zj(p)+Ni[Z1,Z2,…,Zm]=Hi,
where αi,j[p] are the constant coefficients of zj(p), the derivative of zj(j=1,2,…,m) that appears in the ith equation for i=1,2,…,m. Assume that the solution zi(η) of (32) 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 (32) can be linearised using one-term Taylor series for multiple variables, so that (32) can be approximated as
(33)Li[Z1,r+1,…,Zm,r+1]+Ni[⋯]+∑j=1m‍∑p=0ni‍(Zj,r+1(p)-Zj,r(p))∂Ni∂Zj(p)[⋯]=Hi,
where
(34)[⋯]=[Z1,r,Z2,r,…,Zm,r].Equation (33) can be rewritten as
(35)Li[Z1,r+1,…,Zm,r+1]+∑j=1m‍∑p=0ni‍Zj,r+1(p)∂Ni∂Zj(p)[⋯]=Hi+∑j=1m‍∑p=0ni‍Zj,r(p)∂Ni∂Zj(p)[⋯]-Ni[⋯].To solve the iteration scheme (35), we use the Chebyshev spectral collocation method.


## 2.4. Compact Finite Difference Quasilinearization Method

Using compact finite difference schemes to solve the QLM scheme (35) results in the compact finite difference quasilinearization method (CFD-QLM). We will use the sixth-order finite difference schemes the same way we used them for the CFD-RM.


## 3. Examples

In this section, two numerical examples are discussed to compare the CFD results against spectral method results. We consider one-dimensional and three-dimensional problems.


### 3.1. MHD Boundary Layer Flow over a Shrinking Sheet

We first consider a steady one-dimensional laminar flow of a viscous incompressible electrically conducting fluid over a continuously shrinking sheet. The governing equation is given (see, e.g., [29, 30]) in similarity form as
(36)f′′′+ff′′-f′2-M2f′=0
subject to
(37)f(0)=0,f′(0)=-1,f′(∞)=0,
where f=f(η) and M is the magnetic interaction parameter. The analytical solution for (36) is given by
(38)f(η)=1α(e-αη-1),α=M2-1.We solve (36) using the SRM and CFD-RM and the SQLM and CFD-QLM. In the next subsections we explain the development of each of the methods for the solution of (36).


#### 3.1.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (36), we first set f′(η)=g(η) and write the equation as the following system of equations:
(39)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(40)f(0)=0,g(0)=-1,g(∞)=0.Applying the SRM and CFD-RM on (36), we obtain the following iteration scheme:
(41)fr+1′=gr,fr+1(0)=0,(42)gr+1′′+fr+1gr+1′-M2gr+1=gr2,gr+1(0)=-1,gr+1(∞)=0.
The SRM and CFD-RM schemes that result from applying the spectral method and compact finite difference schemes on (41)-(42), respectively, are given by
(43)X1fr+1=Y1,fr+1(τN)=0,X2gr+1=Y2,gr+1(τN)=-1,gr+1(τ0)=0,
where for the SRM,
(44)X1=D,Y1=gr,X2=D2+diag[fr+1]D-M2I,Y2=gr2,
where I is an (N+1)×(N+1) identity matrix and diag is a function that places a vector [] on the main diagonal of an (N+1)×(N+1) matrix of zeros. For the CFD-RM, we obtain
(45)X1=E1,Y1=gr-A1-1K1,X2=E2+diag[fr+1]E1-M2I,Y2=gr2-A2-1K2-A1-1K1,
where E1=A1-1B1 and E2=A2-1B2. I is an (N-1)×(N-1) identity matrix.


#### 3.1.2. SQLM and CFD-QLM

To solve (36) using the SQLM and CFD-QLM, we first linearize the equation using the quasilinearization method. We first reduce the order of the equation by setting f′(η)=g(η) and write it as the following system of equations:
(46)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(47)f(0)=0,g(0)=-1,g(∞)=0.Applying the quasilinearization method on (46), we obtain
(48)fr+1′-gr+1=0,(49)gr′fr+1+gr+1′′+frgr+1′-(2gr+M2)gr+1=frgr′-gr2
subject to
(50)fr+1(0)=0,gr+1(0)=-1,gr+1(0)=0.The SQLM and CFD-QLM schemes that result from applying the spectral method and compact finite difference schemes on (48)-(49), respectively, are given by
(51)[Δ1,1Δ1,2Δ2,1Δ2,2][fr+1gr+1]=[Φ1,rΦ2,r],
where in the framework of the SQLM,
(52)Δ1,1=D,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=D2+diag[fr]D-diag[2gr+M2],Φ1,r=O,Φ2,r=frgr′-gr2,
where O is an (N+1)×1 zero vector. For the CFD-QLM, we have
(53)Δ1,1=E1,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=E2+diag[fr]E1-diag[2gr+M2],Φ1,r=-A1-1K1,Φ2,r=frgr′-gr2-A2-1K2-diag[fr](A1-1K1).


### 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 [23, 31] for this problem are given as the following dimensionless system of equations and boundary conditions:
(54)f′′′+ff′′-(f′)2-Kf′-A(f′+η2f′′)+Grθ+Gcϕ=0,1Prθ′′-f′θ+fθ′-A(θ+12ηθ′)=0,1Scϕ′′-f′ϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0
subject to
(55)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.


#### 3.2.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (54), we set f′(η)=g(η) and obtain the following iteration scheme:
(56)fr+1′=gr,g′′r+1+fr+1gr+1′-(A+K)gr+1+η2gr+1′=gr2-Grθr-Gcϕr,1Prθr+1-gr+1θr+1+fr+1θr+1′-A(θr+1+η2θr+1)=0,1Scϕr+1-gr+1ϕr+1+fr+1ϕr+1′-(A+γ)ϕr+1-Aη2ϕr+1′=0
subject to
(57)fr+1(0)=fw,gr+1(0)=1,gr+1(∞)=0,θr+1(0)=1,θr+1(∞)=0,ϕr+1(0)=1,ϕr+1(∞)=0.Applying the Chebyshev pseudo-spectral method and finite difference schemes, we obtain(58)X1fr+1=Y1,fr+1(0)=fw,X2gr+1=Y2,gr+1(τ0)=1,gr+1(τN)=0,X3θr+1=Y3,θr+1(τ0)=1,θ(τN)=0,X4ϕr+1=Y4,ϕr+1(τ0)=1,ϕ(τN)=0.For the SRM, we have(59)X1=D,Y1=gr,X2=D2+diag[fr+1-Aη2]D-(K+A)I,Y2=gr2-Grθr-Gcϕr,X3=1PrD2+diag[fr+1-Aη2]D-diag[gr+1]-AI,Y3=O,X4=1ScD2+diag[fr+1-Aη2]D-diag[gr+1]-(A+γ)I,Y4=O,
where O is an (N+1)×1 zero vector.For the CFD-RM, we have(60)X1=E1,X2=E2+diag[fr+1-Aη2]E1-(K+A)I,X3=1PrE2+diag[fr+1-Aη2]E1-diag[gr+1]-AI,X4=1ScE2+diag[fr+1-Aη2]E1-diag[gr+1]-(A+γ)I,Y1=gr-A1-1K1,Y2=gr2-Grθr-Gcϕr-A2-1K2-diag[fr+1-Aη2](A1-1K1),Y3=-1Pr(A2-1K2)-diag[fr+1-Aη2](A1-1K1),Y4=-1Sc(A2-1K2)-diag[fr+1-Aη2](A1-1K1).


#### 3.2.2. SQLM and CFD-QLM

To apply the SQLM and CFD-QLM on (54), we first apply the quasilinearization method and obtain the following iteration scheme:
(61)fr+1′′′+(fr-Aη2)fr+1′′-(2fr′+K+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=θr′fr-fr′θr,1Scϕr+1′′+(fr-Aη2)ϕr+1′-(fr′+A+γ)ϕr+1-ϕrfr+1′+ϕr′fr+1=ϕr′fr-fr′ϕr
subject to
(62)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 spectral method and compact finite difference schemes, we obtain the SQLM and CFD-QLM iteration schemes, respectively. In matrix form, the SQLM and CFD-QLM schemes are given by(63)[Δ1,1Δ1,2Δ1,3Δ2,1Δ2,2Δ2,3Δ3,1Δ3,2Δ3,3][fr+1θr+1ϕr+1]=[Φ1,rΦ2,rΦ3,r].For the SQLM, we have(64)Δ1,1=D3-diag[fr-Aη2]D2-diag[2fr′+K+A]D+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]D+diag[θr′],Δ2,2=1PrD2+diag[fr-Aη2]D-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]D+diag[ϕr′],Δ3,2=O,Δ3,3=1ScD2+diag[fr-Aη2]D-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2,Φ2,r=frθr′-fr′θr,Φ3,r=frϕr-fr′ϕr,
where O is an (N+1)×(N+1) zero matrix, and for the CFD-QLM, we have
(65)Δ1,1=E3-diag[fr-Aη2]E2-diag[2fr′+K+A]E1+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]E1+diag[θr′],Δ2,2=1PrE2+diag[fr-Aη2]E1-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]E1+diag[ϕr′],Δ3,2=O,Δ3,3=1ScE2+diag[fr-Aη2]E1-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2-A3-1K3-A3-1Kr′-diag[fr-Aη2](A2-1K2)-diag[fr-Aη2](A2-1K2′)+diag[2fr′](A1-1K1)-(A+K)(A1-1K1)+diag[2fr′](A1-1K1′)-(A+K)(A1-1K1′),Φ2,r=frθr′-fr′θr-1Pr(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[θr](A1-1K1)+diag[θr](A1-1K1′),Φ3,r=frϕr-fr′ϕr-1Sc(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[ϕr](A1-1K1)+diag[ϕr](A1-1K1′),
where O is an (N-1)×(N-1) zero matrix and E3=A3-1B3.


## 3.1. MHD Boundary Layer Flow over a Shrinking Sheet

We first consider a steady one-dimensional laminar flow of a viscous incompressible electrically conducting fluid over a continuously shrinking sheet. The governing equation is given (see, e.g., [29, 30]) in similarity form as
(36)f′′′+ff′′-f′2-M2f′=0
subject to
(37)f(0)=0,f′(0)=-1,f′(∞)=0,
where f=f(η) and M is the magnetic interaction parameter. The analytical solution for (36) is given by
(38)f(η)=1α(e-αη-1),α=M2-1.We solve (36) using the SRM and CFD-RM and the SQLM and CFD-QLM. In the next subsections we explain the development of each of the methods for the solution of (36).


### 3.1.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (36), we first set f′(η)=g(η) and write the equation as the following system of equations:
(39)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(40)f(0)=0,g(0)=-1,g(∞)=0.Applying the SRM and CFD-RM on (36), we obtain the following iteration scheme:
(41)fr+1′=gr,fr+1(0)=0,(42)gr+1′′+fr+1gr+1′-M2gr+1=gr2,gr+1(0)=-1,gr+1(∞)=0.
The SRM and CFD-RM schemes that result from applying the spectral method and compact finite difference schemes on (41)-(42), respectively, are given by
(43)X1fr+1=Y1,fr+1(τN)=0,X2gr+1=Y2,gr+1(τN)=-1,gr+1(τ0)=0,
where for the SRM,
(44)X1=D,Y1=gr,X2=D2+diag[fr+1]D-M2I,Y2=gr2,
where I is an (N+1)×(N+1) identity matrix and diag is a function that places a vector [] on the main diagonal of an (N+1)×(N+1) matrix of zeros. For the CFD-RM, we obtain
(45)X1=E1,Y1=gr-A1-1K1,X2=E2+diag[fr+1]E1-M2I,Y2=gr2-A2-1K2-A1-1K1,
where E1=A1-1B1 and E2=A2-1B2. I is an (N-1)×(N-1) identity matrix.


### 3.1.2. SQLM and CFD-QLM

To solve (36) using the SQLM and CFD-QLM, we first linearize the equation using the quasilinearization method. We first reduce the order of the equation by setting f′(η)=g(η) and write it as the following system of equations:
(46)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(47)f(0)=0,g(0)=-1,g(∞)=0.Applying the quasilinearization method on (46), we obtain
(48)fr+1′-gr+1=0,(49)gr′fr+1+gr+1′′+frgr+1′-(2gr+M2)gr+1=frgr′-gr2
subject to
(50)fr+1(0)=0,gr+1(0)=-1,gr+1(0)=0.The SQLM and CFD-QLM schemes that result from applying the spectral method and compact finite difference schemes on (48)-(49), respectively, are given by
(51)[Δ1,1Δ1,2Δ2,1Δ2,2][fr+1gr+1]=[Φ1,rΦ2,r],
where in the framework of the SQLM,
(52)Δ1,1=D,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=D2+diag[fr]D-diag[2gr+M2],Φ1,r=O,Φ2,r=frgr′-gr2,
where O is an (N+1)×1 zero vector. For the CFD-QLM, we have
(53)Δ1,1=E1,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=E2+diag[fr]E1-diag[2gr+M2],Φ1,r=-A1-1K1,Φ2,r=frgr′-gr2-A2-1K2-diag[fr](A1-1K1).


## 3.1.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (36), we first set f′(η)=g(η) and write the equation as the following system of equations:
(39)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(40)f(0)=0,g(0)=-1,g(∞)=0.Applying the SRM and CFD-RM on (36), we obtain the following iteration scheme:
(41)fr+1′=gr,fr+1(0)=0,(42)gr+1′′+fr+1gr+1′-M2gr+1=gr2,gr+1(0)=-1,gr+1(∞)=0.
The SRM and CFD-RM schemes that result from applying the spectral method and compact finite difference schemes on (41)-(42), respectively, are given by
(43)X1fr+1=Y1,fr+1(τN)=0,X2gr+1=Y2,gr+1(τN)=-1,gr+1(τ0)=0,
where for the SRM,
(44)X1=D,Y1=gr,X2=D2+diag[fr+1]D-M2I,Y2=gr2,
where I is an (N+1)×(N+1) identity matrix and diag is a function that places a vector [] on the main diagonal of an (N+1)×(N+1) matrix of zeros. For the CFD-RM, we obtain
(45)X1=E1,Y1=gr-A1-1K1,X2=E2+diag[fr+1]E1-M2I,Y2=gr2-A2-1K2-A1-1K1,
where E1=A1-1B1 and E2=A2-1B2. I is an (N-1)×(N-1) identity matrix.


## 3.1.2. SQLM and CFD-QLM

To solve (36) using the SQLM and CFD-QLM, we first linearize the equation using the quasilinearization method. We first reduce the order of the equation by setting f′(η)=g(η) and write it as the following system of equations:
(46)f′=g,g′′+fg′-g2-M2g=0
with the boundary conditions
(47)f(0)=0,g(0)=-1,g(∞)=0.Applying the quasilinearization method on (46), we obtain
(48)fr+1′-gr+1=0,(49)gr′fr+1+gr+1′′+frgr+1′-(2gr+M2)gr+1=frgr′-gr2
subject to
(50)fr+1(0)=0,gr+1(0)=-1,gr+1(0)=0.The SQLM and CFD-QLM schemes that result from applying the spectral method and compact finite difference schemes on (48)-(49), respectively, are given by
(51)[Δ1,1Δ1,2Δ2,1Δ2,2][fr+1gr+1]=[Φ1,rΦ2,r],
where in the framework of the SQLM,
(52)Δ1,1=D,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=D2+diag[fr]D-diag[2gr+M2],Φ1,r=O,Φ2,r=frgr′-gr2,
where O is an (N+1)×1 zero vector. For the CFD-QLM, we have
(53)Δ1,1=E1,Δ1,2=-I,Δ2,1=diag[gr′],Δ2,2=E2+diag[fr]E1-diag[2gr+M2],Φ1,r=-A1-1K1,Φ2,r=frgr′-gr2-A2-1K2-diag[fr](A1-1K1).


## 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 [23, 31] for this problem are given as the following dimensionless system of equations and boundary conditions:
(54)f′′′+ff′′-(f′)2-Kf′-A(f′+η2f′′)+Grθ+Gcϕ=0,1Prθ′′-f′θ+fθ′-A(θ+12ηθ′)=0,1Scϕ′′-f′ϕ+fϕ′-A(ϕ+12ηϕ′)-γϕ=0
subject to
(55)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.


### 3.2.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (54), we set f′(η)=g(η) and obtain the following iteration scheme:
(56)fr+1′=gr,g′′r+1+fr+1gr+1′-(A+K)gr+1+η2gr+1′=gr2-Grθr-Gcϕr,1Prθr+1-gr+1θr+1+fr+1θr+1′-A(θr+1+η2θr+1)=0,1Scϕr+1-gr+1ϕr+1+fr+1ϕr+1′-(A+γ)ϕr+1-Aη2ϕr+1′=0
subject to
(57)fr+1(0)=fw,gr+1(0)=1,gr+1(∞)=0,θr+1(0)=1,θr+1(∞)=0,ϕr+1(0)=1,ϕr+1(∞)=0.Applying the Chebyshev pseudo-spectral method and finite difference schemes, we obtain(58)X1fr+1=Y1,fr+1(0)=fw,X2gr+1=Y2,gr+1(τ0)=1,gr+1(τN)=0,X3θr+1=Y3,θr+1(τ0)=1,θ(τN)=0,X4ϕr+1=Y4,ϕr+1(τ0)=1,ϕ(τN)=0.For the SRM, we have(59)X1=D,Y1=gr,X2=D2+diag[fr+1-Aη2]D-(K+A)I,Y2=gr2-Grθr-Gcϕr,X3=1PrD2+diag[fr+1-Aη2]D-diag[gr+1]-AI,Y3=O,X4=1ScD2+diag[fr+1-Aη2]D-diag[gr+1]-(A+γ)I,Y4=O,
where O is an (N+1)×1 zero vector.For the CFD-RM, we have(60)X1=E1,X2=E2+diag[fr+1-Aη2]E1-(K+A)I,X3=1PrE2+diag[fr+1-Aη2]E1-diag[gr+1]-AI,X4=1ScE2+diag[fr+1-Aη2]E1-diag[gr+1]-(A+γ)I,Y1=gr-A1-1K1,Y2=gr2-Grθr-Gcϕr-A2-1K2-diag[fr+1-Aη2](A1-1K1),Y3=-1Pr(A2-1K2)-diag[fr+1-Aη2](A1-1K1),Y4=-1Sc(A2-1K2)-diag[fr+1-Aη2](A1-1K1).


### 3.2.2. SQLM and CFD-QLM

To apply the SQLM and CFD-QLM on (54), we first apply the quasilinearization method and obtain the following iteration scheme:
(61)fr+1′′′+(fr-Aη2)fr+1′′-(2fr′+K+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=θr′fr-fr′θr,1Scϕr+1′′+(fr-Aη2)ϕr+1′-(fr′+A+γ)ϕr+1-ϕrfr+1′+ϕr′fr+1=ϕr′fr-fr′ϕr
subject to
(62)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 spectral method and compact finite difference schemes, we obtain the SQLM and CFD-QLM iteration schemes, respectively. In matrix form, the SQLM and CFD-QLM schemes are given by(63)[Δ1,1Δ1,2Δ1,3Δ2,1Δ2,2Δ2,3Δ3,1Δ3,2Δ3,3][fr+1θr+1ϕr+1]=[Φ1,rΦ2,rΦ3,r].For the SQLM, we have(64)Δ1,1=D3-diag[fr-Aη2]D2-diag[2fr′+K+A]D+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]D+diag[θr′],Δ2,2=1PrD2+diag[fr-Aη2]D-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]D+diag[ϕr′],Δ3,2=O,Δ3,3=1ScD2+diag[fr-Aη2]D-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2,Φ2,r=frθr′-fr′θr,Φ3,r=frϕr-fr′ϕr,
where O is an (N+1)×(N+1) zero matrix, and for the CFD-QLM, we have
(65)Δ1,1=E3-diag[fr-Aη2]E2-diag[2fr′+K+A]E1+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]E1+diag[θr′],Δ2,2=1PrE2+diag[fr-Aη2]E1-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]E1+diag[ϕr′],Δ3,2=O,Δ3,3=1ScE2+diag[fr-Aη2]E1-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2-A3-1K3-A3-1Kr′-diag[fr-Aη2](A2-1K2)-diag[fr-Aη2](A2-1K2′)+diag[2fr′](A1-1K1)-(A+K)(A1-1K1)+diag[2fr′](A1-1K1′)-(A+K)(A1-1K1′),Φ2,r=frθr′-fr′θr-1Pr(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[θr](A1-1K1)+diag[θr](A1-1K1′),Φ3,r=frϕr-fr′ϕr-1Sc(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[ϕr](A1-1K1)+diag[ϕr](A1-1K1′),
where O is an (N-1)×(N-1) zero matrix and E3=A3-1B3.


## 3.2.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (54), we set f′(η)=g(η) and obtain the following iteration scheme:
(56)fr+1′=gr,g′′r+1+fr+1gr+1′-(A+K)gr+1+η2gr+1′=gr2-Grθr-Gcϕr,1Prθr+1-gr+1θr+1+fr+1θr+1′-A(θr+1+η2θr+1)=0,1Scϕr+1-gr+1ϕr+1+fr+1ϕr+1′-(A+γ)ϕr+1-Aη2ϕr+1′=0
subject to
(57)fr+1(0)=fw,gr+1(0)=1,gr+1(∞)=0,θr+1(0)=1,θr+1(∞)=0,ϕr+1(0)=1,ϕr+1(∞)=0.Applying the Chebyshev pseudo-spectral method and finite difference schemes, we obtain(58)X1fr+1=Y1,fr+1(0)=fw,X2gr+1=Y2,gr+1(τ0)=1,gr+1(τN)=0,X3θr+1=Y3,θr+1(τ0)=1,θ(τN)=0,X4ϕr+1=Y4,ϕr+1(τ0)=1,ϕ(τN)=0.For the SRM, we have(59)X1=D,Y1=gr,X2=D2+diag[fr+1-Aη2]D-(K+A)I,Y2=gr2-Grθr-Gcϕr,X3=1PrD2+diag[fr+1-Aη2]D-diag[gr+1]-AI,Y3=O,X4=1ScD2+diag[fr+1-Aη2]D-diag[gr+1]-(A+γ)I,Y4=O,
where O is an (N+1)×1 zero vector.For the CFD-RM, we have(60)X1=E1,X2=E2+diag[fr+1-Aη2]E1-(K+A)I,X3=1PrE2+diag[fr+1-Aη2]E1-diag[gr+1]-AI,X4=1ScE2+diag[fr+1-Aη2]E1-diag[gr+1]-(A+γ)I,Y1=gr-A1-1K1,Y2=gr2-Grθr-Gcϕr-A2-1K2-diag[fr+1-Aη2](A1-1K1),Y3=-1Pr(A2-1K2)-diag[fr+1-Aη2](A1-1K1),Y4=-1Sc(A2-1K2)-diag[fr+1-Aη2](A1-1K1).


## 3.2.2. SQLM and CFD-QLM

To apply the SQLM and CFD-QLM on (54), we first apply the quasilinearization method and obtain the following iteration scheme:
(61)fr+1′′′+(fr-Aη2)fr+1′′-(2fr′+K+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=θr′fr-fr′θr,1Scϕr+1′′+(fr-Aη2)ϕr+1′-(fr′+A+γ)ϕr+1-ϕrfr+1′+ϕr′fr+1=ϕr′fr-fr′ϕr
subject to
(62)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 spectral method and compact finite difference schemes, we obtain the SQLM and CFD-QLM iteration schemes, respectively. In matrix form, the SQLM and CFD-QLM schemes are given by(63)[Δ1,1Δ1,2Δ1,3Δ2,1Δ2,2Δ2,3Δ3,1Δ3,2Δ3,3][fr+1θr+1ϕr+1]=[Φ1,rΦ2,rΦ3,r].For the SQLM, we have(64)Δ1,1=D3-diag[fr-Aη2]D2-diag[2fr′+K+A]D+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]D+diag[θr′],Δ2,2=1PrD2+diag[fr-Aη2]D-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]D+diag[ϕr′],Δ3,2=O,Δ3,3=1ScD2+diag[fr-Aη2]D-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2,Φ2,r=frθr′-fr′θr,Φ3,r=frϕr-fr′ϕr,
where O is an (N+1)×(N+1) zero matrix, and for the CFD-QLM, we have
(65)Δ1,1=E3-diag[fr-Aη2]E2-diag[2fr′+K+A]E1+diag[fr′′],Δ1,2=GrI,Δ1,3=GcI,Δ2,1=-diag[θr]E1+diag[θr′],Δ2,2=1PrE2+diag[fr-Aη2]E1-AI-diag[fr′],Δ2,3=O,Δ3,1=-diag[ϕr]E1+diag[ϕr′],Δ3,2=O,Δ3,3=1ScE2+diag[fr-Aη2]E1-(A+γ)I-diag[fr′],Φ1,r=frfr′′-fr′2-A3-1K3-A3-1Kr′-diag[fr-Aη2](A2-1K2)-diag[fr-Aη2](A2-1K2′)+diag[2fr′](A1-1K1)-(A+K)(A1-1K1)+diag[2fr′](A1-1K1′)-(A+K)(A1-1K1′),Φ2,r=frθr′-fr′θr-1Pr(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[θr](A1-1K1)+diag[θr](A1-1K1′),Φ3,r=frϕr-fr′ϕr-1Sc(A2-1K2)-diag[fr-Aη2](A1-1K1)+diag[ϕr](A1-1K1)+diag[ϕr](A1-1K1′),
where O is an (N-1)×(N-1) zero matrix and E3=A3-1B3.


## 4. Results and Discussion

In this section we present results for Examples 1 and 2 and give a comparison between the compact-finite-difference-based methods and the spectral-method-based methods, that is, the SRM and CFD-RM and the SQLM and CFD-QLM. The solutionf(η) and the velocity profiles f′(η) for varying values of M for Example 1 are shown by Figure 1. The numerical results are compared against the exact analytical solution and good agreement is observed in all cases of varying parameter M. Tables 1 and 2 show the computed skin friction f′′(0) to an accuracy of 10-9 obtained using the four methods.Table 1
Comparison of the SRM and CFD-RM results for the solution off′′(0) in Example 1.
M
Exact
SRM
CFD6-RM
N
f
′
′
(
0
)
CPU time
N
f
′
′
(
0
)
CPU time
2
1.732050808
40
1.732050808
0.04
270
1.732050808
0.23
5
4.898979486
40
4.898979486
0.01
550
4.898979486
1.11
10
9.949874371
40
9.949874371
0.01
700
9.949874371
2.19Table 2
Comparison of the SQLM and CFD-QLM results for the solution off′′(0) in Example 1.
M
Exact
SQLM
CFD6-QLM
N
f
′
′
(
0
)
CPU time
N
f
′
′
(
0
)
CPU time
2
1.732050808
40
1.732050808
0.06
270
1.732050808
0.31
5
4.898979486
40
4.898979486
0.03
550
4.898979486
1.87
10
9.949874371
40
9.949874371
0.01
700
9.949874371
2.84f
(
η
) and the velocity profile f′(η) for M=2,5,10.
(a)
(b)It can be seen from Table1 that in terms of computational speed, the SRM is efficient compared to the CFD-RM. This is because quite a large number of grid points are needed for the CFD-RM to give an accuracy of 10-9 as compared to the SRM. From Table 2 we observe similar results with the SQLM being computationally faster than the CFD-QLM. Again the number of grid points is the main reason for the difference in speed. The spectral method based approaches only require few grid points to give highly accurate results.Figure2 shows the errors at different iterations for the CFD-RM, and CFD-QLM respectively. In all cases the error decreases with each iteration which shows convergence of the methods. Increasing the value of the magnetic interaction parameter, M increases the convergence of both the CFD-RM and CFD-QLM. This may be explained by the observation that when M is very large, the dominant terms in the equation from which the SRM scheme is derived are g′′-M2g. Solving this equation gives a solution of the form g=ce-Mη (where c is a constant). For large values of M, the M2-1 appearing in the given exact solution (38) is approximately equal to M2. Consequently, the exact solution reduces to an exponential equation that has the same form as the approximate equation obtained from the SRM scheme at large M. The CFD-QLM shows a faster convergence than the CFD-RM as shown in Figure 2.Effect of varyingM on the error for the CFD-RM and CFD-QLM.
(a)
(b)Figures3 and 4 show the velocity f′(η), temperature θ(η), and concentration ϕ(η) for Example 2 for varying values of Gr and fw, respectively. We observe that the graphs are qualitatively similar to those reported in [23, 31]. The comparison of the performance of the methods is given by Tables 3 and 4. In Example 2 we compute the solution to an accuracy of 10-8. The SRM is computationally faster than the CFD-RM since it requires fewer grid points to give the required accuracy. We observe different results for the SQLM in Example 2 compared to Example 1. When comparing the results with the SRM, CFD-RM, and CFD-QLM, we observe that the SQLM fails to give results to the accuracy of 10-8 as it can be seen in Table 4. In this case the CFD-QLM is much better than the SQLM in terms of accuracy. This is one instance where the advantage of using CFD over spectral approach to integrate linearised equations is highlighted. The spectral approach is less accurate than the CFD approach when large systems of equations are considered.Table 3
Comparison of the SRM and CFD-RM results for the solution off′′(0) in Example 2.
A
f
w
Gr
K
SRM
CFD-RM
N
f
′
′
(
0
)
CPU time
N
f
′
′
(
0
)
CPU time
0
0
1
1
40
−0.63371481
0.43
250
−0.63371481
0.86
2
0
1
1
40
−1.30748109
0.50
200
−1.30748109
0.26
4
0
1
1
40
−1.79277127
0.55
230
−1.79277127
0.26
1
1
1
1
40
−1.55880093
0.12
210
−1.55880093
0.34
1
2
1
1
40
−2.31812942
0.23
340
−2.31812942
0.88
1
4
1
1
40
−4.14430512
0.29
410
−4.14430512
1.02
1
0
0
1
40
−1.32052206
0.11
180
−1.32052206
0.19
1
0
5
1
40
0.17151463
0.28
280
0.17151463
1.09
1
0
10
1
40
1.49371495
0.52
500
1.49371495
8.47
1
0
1
2
40
−0.61644148
0.18
200
−0.61644148
0.44
1
0
1
4
40
−1.84589956
0.27
310
−1.84589956
0.74
1
0
1
8
40
−2.65350852
0.36
620
−2.65350852
6.12Table 4
Comparison of the SQLM and CFD-QLM results for the solution off′′(0) in Example 2.
A
f
w
Gr
K
SQLM
CFD-QLM
N
f
′
′
(
0
)
CPU time
N
f
′
′
(
0
)
CPU time
0
1
1
1
200
−0.63371481
64.49
210
−0.63371481
0.27
2
1
1
1
200
−1.44059181
23.20
490
−1.30748109
1.68
4
1
1
1
200
−3.62469942
20.89
650
−1.79277127
3.60
1
1
1
1
200
−1.55886443
64.09
230
−1.55880093
0.35
1
2
1
1
200
−2.31813091
85.57
310
−2.31812942
0.83
1
4
1
1
200
−4.14430509
105.76
440
−4.14430512
2.03
1
1
0
1
200
−1.32106621
26.51
300
−1.32052206
0.69
1
1
5
1
200
0.17105418
33.92
300
0.17151463
0.67
1
1
10
1
200
1.49345960
36.09
420
1.49371495
1.77
1
1
1
2
200
−0.61872790
31.20
360
−0.61644148
1.24
1
1
1
4
200
−2.07158002
33.19
380
−1.84589956
1.20
1
1
1
8
200
−2.98816061
27.64
450
−2.65350852
2.01Velocityf(η), temperature θ(η), and concentration ϕ(η) profiles in Example 2 for various values of Gr(Gr=0,5,20) when Pr=0.72,Sc=0.6,A=0.5,γ=1,Gc=1,fw=1, and K=1.
(a)
(b)
(c)Velocityf(η), temperature θ(η), and concentration ϕ(η) profiles in Example 2 for various values of fw(fw=0,2,4) when Pr=0.72,Sc=0.6,A=0.5,γ=1,Gc=1,Gr=1, and K=1.
(a)
(b)In terms of convergence between the CFD-RM and CFD-QLM, we observe a faster convergence with the CFD-QLM than the CFD-RM. This is depicted in Figures5 and 6.Effect of varyingGr on the errors for the CFD-RM and the CFD-QLM.
(a)
(b)Effect of varyingfw on the errors for the CFD-RM and the CFD-QLM.
(a)
(b)


## 5. Conclusion

In this work we have introduced two new methods for solving systems of nonlinear boundary value equations. These methods called the compact finite difference relaxation method (CFD-RM) and compact finite difference quasilinearization method (CFD-QLM) are a modification of the spectral relaxation method (SRM) and the spectral quasilinearization method (SQLM), respectively. We compared the CFD schemes and spectral methods in solving nonlinear similarity boundary layer problems by comparing the CFD-based CFD-RM and CFD-QLM against the spectral-method-based SRM and SQLM. We did the comparison on one dimensional and  three-dimensional problems. After comparing these methods we can conclude that they are all highly accurate with the spectral method outperforming the CFD in terms of computational speed. The SQLM though was less accurate for the three-dimensional problems. The CFD-QLM was able to handle the three-dimensional problem. We also observed that the CFD-QLM converges faster than the CFD-RM. When solving nonlinear boundary value problems, we recommend that the spectral-method-based algorithms must be used in problems involving small systems of equations and when computational speed is of importance because they are easy to code and computationally faster. CFD-based algorithms are ideal for solving large systems of equations when high accuracy is required.



---

*Source: 746489-2013-08-27.xml*