# On an Interpolation Based Spectral Homotopy Analysis Method for PDE Based Unsteady Boundary Layer Flows

Authors: S. S. Motsa
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This work presents a new approach to the application of the spectral homotopy analysis method (SHAM) in solving non-linear partial differential equations (PDEs). The proposed approach is based on an innovative idea of seeking solutions that obey a rule of solution expression that is defined in terms of bivariate Lagrange interpolation polynomials. The applicability and effectiveness of the expanded SHAM approach are tested on a non-linear PDE that models the problem of unsteady boundary layer flow caused by an impulsively stretching plate. Numerical simulations are conducted to generate results for the important flow properties such as the local skin friction. The accuracy of the present results is validated against existing results from the literature and against results generated using the Keller-box method. The preliminary results from the proposed study indicate that the present method is more accurate and computationally efficient than more traditional methods used for solving PDEs that describe nonsimilar boundary layer flow.


---

## Body

## 1. Introduction

In this work, a new approach for applying the spectral homotopy analysis method (SHAM) is used to solve a partial differential equation (PDE) that models the problem of unsteady boundary layer flow caused by an impulsively stretching plate. The SHAM is a discrete numerical version of the homotopy analysis method (HAM) that has been widely applied to solve a wide variety of nonlinear ordinary and partial differential equations with applications in applied mathematics, physics, nonlinear mechanics, finance, and engineering. A detailed systematic description of the HAM and its applications can be found in two books (and a huge list of references cited therein) [1, 2] by Liao who is credited with developing the method.The SHAM was introduced in [3, 4] and it uses the principles of the traditional HAM and combines them with the Chebyshev spectral collocation method which is used to solve the so-called higher order deformation equations. One of the advantages of the SHAM is that it can accommodate very complex linear operators in its solution algorithm. In a recent application of the SHAM on solving PDE based problems [5], it was observed that the linear operator that gives the best results is one that is selected as the entire linear part of the governing PDEs. This is in sharp contrast to the traditional HAM which can only admit simple ordinary differential equation based linear differential operators with constant coefficient in its solution of PDEs describing the unsteady boundary layer problems of the type discussed in this paper (see, e.g., [6–17]). In the work of [5], the commonly used ODE based linear operator was identified as a limitation if the governing equation was a nonlinear PDE. Partial derivative linear operators were observed to significantly improve the accuracy and convergence speed of the SHAM approach when solving nonlinear partial differential equations. However, it was also observed that the use of PDE based linear operators leads to a complicated linear sequence of higher order deformation equations whose closed form solution was impossible to find. Consequently, the Chebyshev spectral collocation method [18, 19] coupled with a series form of solution constructed from monomial basis functions ({1,ξ,ξ2,…}) was used to solve the higher order deformation equations. However, as it will be pointed out in this work, this approach is only applicable to cases when ξ<1 and becomes less accurate when ξ is very close to 1. Furthermore, the approach of [5] cannot be used directly to solve evolution equations of the form ut=G(u), where u(x,t) is the solution and G(u) is a nonlinear operator which contains all the spatial derivatives of u. In this investigation, we propose the use of Lagrange interpolation polynomials in place of the monomial basis functions used in [5] for the solution of the PDEs that describe the unsteady boundary layer flow over an impulsively stretching sheet.The main objective of this work is to introduce an extended application of the SHAM that uses bivariate Lagrange interpolation polynomials as an alternative method for solving nonlinear PDEs. To validate accuracy of the present SHAM results, the governing PDEs are also solved using the Keller-box implicit finite difference approach. The proposed approach is found to be desirable because it is accurate, computationally efficient, and easy to use.


## 2. Governing Equations

We consider an unsteady boundary layer developed by an impulsively stretching plate in a constant pressure viscous flow. The governing partial differential equations can be obtained by using the standard stream function formulation in conjunction with the transformations suggested by Williams III and Rhyne [20]. The dimensionless governing equation is obtained (see [10, 21, 22] for details) as
(1)f′′′+η2(1-ξ)f′′+ξ(ff′′-(f′)2)-ξ(1-ξ)∂f′∂ξ=0,
subject to the boundary conditions
(2)f(0,ξ)=0,f′(0,ξ)=1,f′(∞,ξ)=0,
where the primes denote differentiation with respect to η and ξ∈[0,1] is the dimensionless time-scale defined as
(3)ξ=1-e-τ,τ=bt,
with b a positive constant and t is the time variable. In the analysis of boundary layer flow problems, a quantity that is of physical interest is the skin friction which is given in this model [10, 21, 22], in dimensionless form, as
(4)CfRex1/2=ξ-1/2f′′(ξ,0),
where Rex is the local Reynolds number.The governing equation (1) can be considered to have an initial solution at ξ=0 which is obtained as a solution of the equation
(5)f′′′+12ηf′′=0,f(0,0)=0,f′(0,0)=1,f′(∞,0)=0.
Solving (5) gives
(6)f(η,0)=ηerfc(η2)+2π[1-exp(-η24)],
where erfc(η) is the standard complementary error function defined by
(7)erfc(η)=2π∫η∞‍exp(-z2)dz.


## 3. Method of Solution

The solution method employed for solving the governing nonlinear PDE (1) is summarised in the following steps.(1)
Convert the nonlinear PDE (1) into sequence nonlinear ODEs using Lagrange polynomial interpolation in the ξ variable.(2)
Decompose the nonlinear sequence of ODEs into a sequence of linear ODEs using the Homotopy analysis method (HAM) approach.(3)
Solve the resulting sequence of linear ODEs using the Chebyshev collocation spectral method.We approximate the exact solution off(η,ξ) by a Lagrange form of polynomial fM(η,ξ) which interpolates f(η,ξ) at the selected points (called collocation points):
(8)0=ξ0<ξ1<ξ2<⋯<ξM=1.
The approximation fM(η,ξ) has the form
(9)fM(η,ξ)=∑j=0M‍gj(η)Lj(ξ),
where gj(η)=fM(η,ξj) and Lj(ξ) are the characteristic Lagrange cardinal polynomials defined as
(10)Lj(ξ)=∏j=0j≠kMξ-ξkξj-ξk
that obey the Kronecker delta equation
(11)Lj(ξk)=δjk={0ifj≠k1ifj=k.The equations for the solution ofgj(η) are obtained by substituting (9) in (1) and requiring that the equation be satisfied exactly at the collocation points ξi,  i=0,1,2,…,M. An important step in the substitution process is the evaluation of the derivatives of Lj(ξ) with respect to ξ. The derivatives can be evaluated analytically if the interpolating points are chosen as the Chebyshev-Gauss-Lobatto points ζi, such that
(12)ξM-i=ζi+12,ζi=cosiπM,i=0,1,2,…,M.
The values of the derivatives at the Chebyshev-Gauss-Lobatto points are computed as
(13)∂f′∂ξ|ξ=ξi=2∑j=0M‍di,jg′(η,ζj)=2∑j=0M‍di,jgj′(η),hhhhhhhhhhhhhhi=0,1,2,…,M,
where di,j  (i,j=0,1,…,M) are entries of the standard Chebyshev differentiation matrix d=[di,j] of size (M+1)×(M+1) (see, e.g., [18, 19]).Thus, by substituting (9) in (1) and making use of the spectral differentiation matrices, we obtain
(14)gi′′′+12(1-ξi)ηgi′′-2ξi(1-ξi)∑j=0M‍di,jgj′+ξi(gigi′′-(gi′)2)=0,
subject to the boundary conditions
(15)gi(0)=0,gi′(0)=1,gi′(∞)=0,i=0,1,…,M.
Equations (14) form a sequence of M+1 nonlinear ordinary differential equations that are to be solved for gi subject to the boundary conditions (15) for i=0,1,…,M. We remark that, since the solution for f(η,ξ) is known exactly when ξ=0 (corresponding to ξ0=ζM=0), we need only solve (14) up to M-1. The next step in the solution procedure is the linearisation of nonlinear ODE system (14) and (15). This is achieved by employing the Homotopy analysis method [1, 2].In the framework of the HAM, we begin by choosing a suitable linear operator which, as suggested in [3, 4], should be taken to be the entire linear component of the given differential equation. Thus, we write the governing equation as:
(16)𝒩(gi)=ℒ(gi)+ξi(gigi′′-(gi′)2)=2ξi(1-ξi)di,MgM′,
where gM′ is a known function that is obtained as the first derivative of (6) and the linear and nonlinear operators are defined as follows:
(17)ℒ(gi)=gi′′′+η2(1-ξi)gi′′-2ξi(1-ξi)∑j=0M-1‍dijgj′,𝒩(gi)=ℒ(gi)+ξi(gigi′′-(gi′)2).The next step of the HAM decomposition involves the construction of the so-called zeroth-order deformation equation which, for (14), is defined as
(18)(1-q)ℒ[Gi(η;q)-gi,0(η)]=qℏ{𝒩[Gi(η;q)]-2ξi(1-ξi)di,MgM′},
subject to the boundary conditions
(19)Gi(0;q)=0,Gi′(0;q)=0,Gi(∞;q)=0,
where q∈[0,1] is an embedding parameter, ℏ denotes a nonzero convergence controlling auxiliary parameter, and gi,0(η) is the initial approximation of the solution of gi for i=0,1,2,…,M-1. The initial approximation is chosen in such a way that
(20)ℒ[gi,0]=2ξi(1-ξi)di,MgM′,gi,0(0)=0,hhhhhhhhhhgi,0′=1,gi,0′(∞)=0.
We remark that, as it was proposed in the SHAM application in nonlinear ODEs [3, 4], it is convenient to obtain the initial guess for the SHAM algorithm as the solution of the equation formed by considering only the linear components of the given nonlinear equation. It can be noted from the zeroth-order deformation equation (18) that as q increases from 0 to 1, Gi(η;q) varies from the initial approximation gi,0(η) to the solution gi(η) of the original equation (14). Expanding Gi(η;q) using Taylor series about q gives
(21)Gi(η;q)=Gi(η;0)+∑n=1+∞‍gi,n(η)qn,
where
(22)gi,n(η)=1n!∂nGi(η;q)∂qn|q=0.
Thus, since Gi(η;1)=gi(η) and Gi(η;0)=gi,0(η,ξ) we obtain
(23)gi(η)=gi,0(η)+∑m=1+∞‍gi,m(η).
The series (23) converges when the auxiliary parameter ℏ is properly chosen.The functionsgi,m, m≥1 are obtained as a solution of higher order deformations equations. These higher order deformation equations are obtained by differentiating the zero-order deformation equations (18) m times with respect to q, then dividing by m!, and finally setting q=0. This gives
(24)ℒ[gi,m(η)-χmgi,m-1(η)]=ℏRi,m,
subject to the boundary conditions
(25)gi,m(0)=0,gi,m′(0)=0,gi,m′(∞)=0,
where
(26)Ri,m=ℒ[gi,m-1]-2(1-χm)ξi(1-ξi)di,MgM′=+ξi∑n=0m-1‍[gi,m-1-ngi,n′′-gi,m-1-n′gi,n′],χm={0,m≤1,1,m>1.The solution for the initial approximationgi,0 and higher order deformation equation giving gi,m are obtained using the Chebyshev spectral collocation method which is applied independently in the η direction using N+1 Chebyshev-Gauss-Lobatto points 0=η0<η1<⋯<ηN=ηe defined as
(27)ηN-r=ηeη¯r+12,η¯r=cosrπN,r=0,1,2,…,N,
where ηe is a finite value that is chosen to be sufficiently large to approximate the conditions at ∞. The derivatives with respect to η are defined in terms of the Chebyshev differentiation matrix as
(28)dpgi,mdηp|η=ηr=(2ηe)p∑s=0NDr,spgi,m(η¯s)=DpGi,m,hhhhhhhhhhhhhhr=0,1,2,…,N,
where p is the order of the derivative, D=(2/ηe)[Dr,s]  (r,s=0,1,2,…,N) with [Dr,s] being an (N+1)×(N+1) Chebyshev derivative matrix, and the vector Gi,m is defined as
(29)Gi,m=[gi,m(η¯0),gi,m(η¯1),…,gi,m(η¯N)]t.
Thus substituting (28) in the equation that gives the initial approximation (20) gives the following M(N+1)×M(N+1) matrix system:
(30)[A0,0A0,1⋯A0,M-1A1,0A1,1⋯A1,M-1⋮⋮⋱⋮AM-1,0AM-1,1⋯AM-1,M-1][G0,0G1,0⋮GM-1,0]=[R0,0R1,0⋮RM-1,0],
where
(31)Ai,i=D3+12(1-ξi)ηD2-2ξi(1-ξi)di,iD,hhhhhhhhhhhhhhhhhi=0,1,…,M-1,(32)Ai,j=-2ξi(1-ξi)di,jD,wheni≠j,(33)Ri,0=2ξi(1-ξi)di,MDGM,i=0,1,…,M-1,(34)η=diag([η0,η1,…,ηN]t).Solving equation (30) gives the initial approximation f0(η,ξ) which is computed as
(35)f0(η,ξ)≈f0M(η,ξ)=∑j=0M‍gj,0(η)Lj(ξ).
To obtain the approximate solution for fm(η,ξ), (for m≥1) the spectral collocation method, with discretisation in the η direction, is applied to (24). This gives the following M(N+1)×M(N+1) matrix system:
(36)[A0,0A0,1⋯A0,M-1A1,0A1,1⋯A1,M-1⋮⋮⋱⋮AM-1,0AM-1,1⋯AM-1,M-1][G0,mG1,m⋮GM-1,m]=[R0,mR1,m⋮RM-1,m],
where Ai,i and Ai,j are as defined in (31) and (32) and
(37)Ri,m=(ℏ+χm)(∑j=0M-1‍di,jDGj,m-1D3Gi,m-1+12(1-ξi)ηD2Gi,m-1=hhh(ℏ+χm)-2ξi(1-ξi)∑j=0M-1‍di,jDGj,m-1)=-2(1-χm)ℏξi(1-ξi)di,MDGM=+ξiℏ∑n=0m-1[Gi,m-1-n(D2Gi,n)hhhhhhhhhh-(DGi,m-1-n)(DGi,n)Gi,m-1-n(D2Gi,n)].The above procedure can easily be extended to any nonlinear partial differential equation system. In the context of the HAM approach the steps involved can be succinctly stated in terms of the HAM jargon as simplylooking for a solution that obeys the following rule of solution expression:
(38)∑m=0+∞‍∑i=0N‍∑j=0M‍fmM,N(η¯i,ξ¯j)Li(η¯)Lj(ξ¯),
where fmM,N(η¯,ξ¯) is a Lagrange form of polynomial which interpolates fm(η,ξ) independently at selected points in both the η and ξ directions defined as
(39){η¯i}={cos(πiN)}i=0N,{ξ¯j}={cos(πjM)}j=0M.
We remark that both η¯ and ξ¯ are defined between [-1,1] and suitable linear transformations are employed to transform the original domain of η and ξ to [-1,1].


## 4. Results and Discussion

In this section we present the numerical solutions of the governing partial differential equations (1) that were computed using the proposed bivariate spectral homotopy analysis method. Starting from the initial analytical solutions at ξ=0 (corresponding to τ=0), the SHAM was used to generate results up to solutions near the steady state values at ξ=1 (corresponding to τ→∞). The accuracy of the computed SHAM approximate results was confirmed against numerical results obtained by using the popular Keller-box implicit finite difference method as described by [23]. The Keller-box method is known to be accurate, fast, and easier to program for boundary layer flow problems of the type discussed in this work. The Keller-box method is still the preferred solution method for most researchers (see, e.g., [24, 25]). The algorithm of the Keller-box method begins with the reduction of the governing nonlinear PDEs into a system of first-order equations that are discretized using central differences. The nonlinear algebraic difference equations are linearised using Newton’s method and written in matrix-vector form. The linear matrix systems are solved in an efficient manner using a block-tridiagonal-elimination technique. The grid spacing in both the η- and ξ-direction is carefully selected to ensure that the Keller-box computations yield consistent results for the governing velocity and temperature distributions to a convergence level of at least 10-7. On the other hand, the SHAM was implemented with N=60 collocation points in the η variable. Furthermore, the finite value used to approximate the boundary conditions at infinity was set to be η∞=30.Using finite terms of the SHAM series we defineKth order approximation at the collocation points η¯i and ξ¯j (for i=0,1,2,…,N and j=0,1,2,…,M) as
(40)f(η,ξ)≈f0(η,ξ)+∑m=1K‍fm(η,ξ)=∑m=0K‍∑i=0N‍∑j=0M‍fmM,N(η¯i,ξ¯j).
Assuming that F(ηi,ξj) is the SHAM approximate solution obtained using (40) at the collocation (grid) points, the residual error is defined as
(41)Res(f)=𝒩[F(ηi,ξj;ℏ)],
where 𝒩 is defined as
(42)𝒩(f)=∂3f∂η3+12(1-ξ)η∂2f∂η2-ξ(1-ξ)∂2f∂ξ∂η-ξ[f∂2f∂η2-(∂f∂η)2].A crucial factor in obtaining accurate and converging SHAM series solutions is the selection of suitable convergence controlling parameterℏ. The infinity norm
(43)Er(ℏ)=∥𝒩[F(ηi,ξ;ℏ)]∥∞
of the residual error at particular values of ξ was used to identify the optimal value of ℏ that gives the best accuracy. Figure 1 is an illustration of a typical residual error curve that can be used to calculate the optimal value of ℏ when ξ=0.98. The optimal values of ℏ are chosen to be the clearly defined minimum of the residual curve. It can be seen from Figure 1 that the optimal ℏ value lies in the range -0.6<ℏ<-0.5. We also note that the residual error decreases with an increase in the number of collocation (grid) points (M) used in interpolating in the ξ direction.Figure 1
Residual error curve for locating optimalℏ when ξ=0.98 for M=5,10,15,20,25,30.Figure2 gives the residual error in the whole range of ξ using a fixed value of ℏ and different collocation points M in the ξ direction. We remark that similar curves are obtained for any fixed value of ℏ, but, for illustration purposes, ℏ=-0.5 was used in Figure 2. It can be seen from the figure that the error decreases with an increase in the number of collocation points used. Another interesting observation to emerge from Figure 2 is that the residual error appears to be nearly uniform for most parts of the domain of ξ. This result is what makes the present SHAM approach with interpolation better than the SHAM version proposed in [5] for the solution of PDEs. We remark that convergence of the SHAM version used in [5] was seen to significantly slow down when ξ was near 1. It can be clearly seen from Figure 2 that the present approach does not suffer from the same limitation.Figure 2
Residual error againstξ for M=5,10,15,20,25,30 when ℏ=-0.5.In Table1, the results for the skin friction f′′(0,ξ) are given for different values of time ξ. Table 1 also gives the number of collocation points (M) and the computational time required to obtain a solution that is consistent with at least six decimal places. The SHAM results match the Keller-box numerical results exactly for all values of ξ. It can be observed from this table that converged solutions are reached using small number of collocation points M for the selected values of ξ. In addition, from the comparison, of the computational times, it is clear that the proposed SHAM approach is more computationally efficient in terms of the amount of time and it takes the method to give the desired results. Consequently, it can be inferred from Table 1 that the SHAM is more efficient than the Keller-box in computing an accurate solution for (1). It is worth mentioning that the apparent computational speed of the proposed SHAM can be explained by the fact that, unlike the Keller-box and other numerical methods, very few grid points in both the η and ξ directions are required to give very accurate results.Table 1
Comparison between the SHAM and Keller-box numerical values of the skin frictionf′′(0,ξ) at different values of time ξ when ℏ=-0.6.
ξ
M
SHAM
Keller-box
f
′
′
(
0
,
ξ
)
Time (sec)
f
′
′
(
0
,
ξ
)
Time (sec)
0.1
2
−0.610468
0.02
−0.610468
3.96
0.3
4
−0.701267
0.05
−0.701267
11.45
0.5
4
−0.789828
0.05
−0.789828
18.78
0.7
4
−0.876267
0.05
−0.876266
26.17
0.8
4
−0.918701
0.05
−0.918701
29.88
0.9
8
−0.960538
0.2
−0.960538
34.05


## 5. Conclusion

In this paper, a new approach based on the spectral homotopy analysis method with bivariate Lagrange interpolating polynomials was introduced for the solution of partial differential equations. The applicability of the proposed method was tested on the problem of unsteady boundary layer flows caused by an impulsively stretching sheet. A residual error analysis was conducted in order to assess the accuracy of the present method. Computational efficiency of the method was demonstrated by comparing with results obtained using the Keller-Box implicit finite difference method. It was found that the proposed SHAM approach was significantly much faster than the Keller-box method. The numerical results presented in this study clearly demonstrate the potential of the proposed bivariate SHAM approach for the simulation of PDEs with high efficiency and accuracy. In future studies it would be interesting to explore the use of this method in other classes of nonlinear PDEs.



---

*Source: 848170-2014-01-02.xml*