On Extending the Quasilinearization Method to Higher Order Convergent Hybrid Schemes Using the Spectral Homotopy Analysis Method

Abstract

We propose a sequence of highly accurate higher order convergent iterative schemes by embedding the quasilinearization algorithm within a spectral collocation method. The iterative schemes are simple to use and significantly reduce the time and number of iterations required to find solutions of highly nonlinear boundary value problems to any arbitrary level of accuracy. The accuracy and convergence properties of the proposed algorithms are tested numerically by solving three Falkner-Skan type boundary layer flow problems and comparing the results to the most accurate results currently available in the literature. We show, for instance, that precision of up to 29 significant figures can be attained with no more than 5 iterations of each algorithm.

1. Introduction

The quasilinearization method (QLM) was originally developed by Bellman and Kalaba [1] as a generalization of the Newton-Raphson method to provide lower and upper bound solutions of nonlinear differential equations. The attraction of quasilinearization is that the algorithm is easy to understand and the method generally converges rapidly if the initial guess is close to the true solution.

Bellman and Kalaba [1] established that the method converges quadratically. However, the original proof of quadratic convergence was subject to restrictive conditions of small step size and convexity or concavity of nonlinear functions, Maleknejad and Najafi [2]. These conditions were subsequently relaxed and the method generalized to be applicable to a wider class of problems; see, for instance, papers by Mandelzweig and his coworkers [3–6] and Lakshmikantham [7, 8]. Parand et al. [9] used the quasilinearization method to solve Volterra's model for population growth in a closed system. Other uses of the quasilinearization method include application to reaction diffusion equations, Jiang and Vatsala [10], and to Volterra integro-differential equations, Ahmad [11], Pandit [12], and Ramos [13].

An often noted disadvantage of quasilinearization is the instability of the method whenever a poor initial guess is chosen, Tuffuor and Labadie [14]. To improve the accuracy and convergence of the quasilinearization method for all initial guesses, we propose in this paper to embed the QLM algorithm within the spectral homotopy analysis method (SHAM) to obtain a sequence of integration schemes with arbitrary higher order convergence.

The spectral homotopy analysis method was introduced by Motsa et al. [15, 16] to address some limitations of the standard homotopy analysis method of Liao [17, 18] by, for example, improving the rate of convergence and extending the region of validity of solutions. The SHAM has been used to solve nonlinear equations that arise in the study of fluid flow problems and other areas of science and engineering, Sibanda et al. [19] and Motsa and Sibanda [20].

Abbasbandy [21] and Chun [22] proposed and studied several methods for nonlinear equations with higher order convergence by using the Adomian decomposition technique [23–25]. Higher order Newton-like iteration formulae for the computation of the solutions of nonlinear equations were also derived in Chun [26] using the homotopy analysis method and by [27] using the homotopy perturbation method. In this paper we extend the ideas used for the solution of nonlinear equations to obtain higher order iteration schemes for solving nonlinear boundary value problems. We propose an extension of the quasilinearization method by using the spectral homotopy analysis method within the QLM algorithm to obtain a sequence of highly accurate and convergent higher order iterative schemes for solving boundary value problems. For illustration purposes we have presented three QLM-SHAM hybrid iteration schemes that are used to solve the Blasius and Falkner-Skan equations. The results are compared to the most accurate skin friction coefficients currently available in the literature by Boyd [28] and Ganapol [29]. Ganapol [29] used an algorithm based on a Maclaurin series with Wynn-epsilon convergence acceleration and analytical continuation to obtain highly accurate skin friction coefficients for the Blasius and Falkner-Skan boundary layer flows. The schemes derived in this paper however require neither convergence acceleration nor analytical continuation to remain steady and accurate for up to 29 digits of precision. In addition, the present schemes are highly efficient with 29-digit precision achieved with five or fewer iterations as compared with at least 104 iterations of Ganapol's algorithm.

The structure of this paper is as follows. Section 2 gives a general framework for the derivation of the hybrid quasilinearization-SHAM schemes for the solution of nonlinear differential equations. Section 3 illustrates the application of the three schemes derived in this paper to the solution of Blasius and Falkner-Skan equations. In Section 4, the results are presented and comparison made with the most accurate skin friction results for Blasius and Falkner-Skan equations currently available in the literature.

2. Derivation of the Iterative Schemes

In this section we present a framework for the derivation of general QLM-SHAM iterative schemes for solving one-dimensional nonlinear differential equations. We consider a general 𝑛-order nonlinear ordinary differential equation of the form 𝐿𝑦(𝑥),𝑦⁽¹⁾(𝑥),𝑦⁽²⁾(𝑥),…,𝑦^(𝑛)+𝐹𝑦(𝑥),𝑦⁽¹⁾(𝑥),𝑦⁽²⁾(𝑥),…,𝑦^(𝑛)=𝜓(𝑥), (1) where 𝜓(𝑥) is a known function of the independent variable 𝑥 and 𝑦(𝑥) is an unknown function. The functions 𝐿 and 𝐹 represent the linear and nonlinear components of the governing equation, respectively. We assume that (1) is to be solved for 𝑥∈[𝑎,𝑏] subject to the boundary conditions 𝐵_𝑎𝑦(𝑎),𝑦⁽¹⁾(𝑎),…,𝑦^(𝑛−1)(𝑎)𝐵=0,_𝑏𝑦(𝑏),𝑦⁽¹⁾(𝑏),…,𝑦^(𝑛−1)(𝑏)=0, (2) where 𝐵_𝑎 and 𝐵_𝑏 are linear operators.

Following [22, 27], we assume that the true solution of (1) is 𝑦_𝛼(𝑥) and that 𝑦_𝛾(𝑥) is an initial approximation that is sufficiently close to 𝑦_𝛼(𝑥). After expanding 𝐹 using Taylor series up to first order about 𝑦_𝛾,𝑦^′_𝛾,…,𝑦_𝛾^(𝑛) we obtain the following coupled system: 𝐿𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)+𝐹𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)+^𝑛_𝑠=0𝑦^(𝑠)−𝑦_𝛾^(𝑠)𝜕𝐹𝜕𝑦^(𝑠)(⋯)+𝐺𝑦,…,𝑦^(𝑛)(𝑥)𝐺=𝜓𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)=𝐹𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)−𝐹𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)−^𝑛_𝑠=0𝑦^(𝑠)−𝑦_𝛾^(𝑠)𝜕𝐹𝜕𝑦^(𝑠)(⋯), (3)(4) where, for compactness, (⋯) denotes (𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)). Note that adding (3) and (4) gives (1). Equation (3) can be rewritten in the form ℒ₁𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)+𝒢₁𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)=Φ𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛), (5) where ℒ₁𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)=𝐿𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)+^𝑛_𝑠=0𝑦^(𝑠)𝜕𝐹𝜕𝑦^(𝑠)𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛),Φ𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)=^𝑛_𝑠=0𝑦_𝛾^(𝑠)𝜕𝐹𝜕𝑦^(𝑠)𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)−𝐹𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛)(𝑥),𝒢+𝜓₁𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛)=𝐺𝑦,𝑦⁽¹⁾,…,𝑦^(𝑛). (6)(7)(8)

Here 𝒢₁ is a nonlinear function that is decomposed using the spectral homotopy analysis method [15, 16]. We define the following zeroth-order deformation equations: 1−𝑞ℒ₁𝑌𝑥;𝑞−𝑦₀(𝑥)=𝑞ℏ𝒩₁𝑌𝑥;𝑞−Φ𝑦_𝛾, (9) where 𝑞∈[0,1] is an embedding parameter, 𝑌(𝑥;𝑞) are unknown functions, ℏ is the convergence controlling parameter, and the “bar” has been introduced for convenience to denote the associated function and its 𝑛 derivatives. For example, 𝑦≡𝑦,𝑦⁽¹⁾,𝑦⁽²⁾,…,𝑦^(𝑛). (10)

The nonlinear operator 𝒩₁ is defined by 𝒩₁𝑌𝑥;𝑞=ℒ₁𝑌_𝑟𝑡;𝑞+𝒢₁𝑌𝑥;𝑞. (11)

By differentiating the zeroth-order equations (11) 𝑚 times with respect to 𝑞, setting 𝑞=0, and finally dividing the resulting equations by 𝑚! (see, e.g., [17, 30–32]), we obtain the following 𝑚th order deformation equations:ℒ₁𝑦_𝑚(𝑥)−𝜒_𝑚+ℏ𝑦_𝑚−1(𝑥)=ℏ𝑅_𝑚−1𝑦₀,𝑦₁,…,𝑦_𝑚−1, (12) where 𝑅_𝑚−1𝑦₀,𝑦₁,…,𝑦_𝑚−1=1()!𝜕𝑚−1^𝑚−1𝒢₁𝑌𝑥;𝑞−Φ𝑦_𝛾𝜕𝑞^𝑚−1_𝑞=0,𝜒_𝑚=0,𝑚⩽1,1,𝑚>1. (13)

After obtaining solutions for (12), the approximate solution for 𝑦(𝑥) is determined as the series solution 𝑦(𝑥)=^∞_𝑚=0𝑦_𝑚(𝑥). (14)

The SHAM solution is said to be of order 𝑀 if the previous series is truncated at 𝑚=𝑀, that is, if 𝑦(𝑥)=^𝑀_𝑚=0𝑦_𝑚(𝑥). (15)

The initial approximation 𝑦₀ required for solving the sequence of linear higher order deformation equations (12) is chosen as the solution that results from solving the linear part of (5) subject to the given boundary conditions (2). That is, we solve ℒ₁𝑦₀,𝑦₀⁽¹⁾,…,𝑦₀^(𝑛)=Φ𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛). (16)

We note that with ℒ₁ as defined in (6), (16) cannot be solved exactly by means of analytical techniques. Numerical methods such as finite differences, finite element method, and spectral method can be used to solve equations of the form (16). Thus, if only the initial approximation is used to approximate the solution 𝑦(𝑥) of the governing nonlinear differential equation (1), that is, if 𝑦(𝑥)≈𝑦₀(𝑥), the (𝑟+1)th approximation of (1) is a solution of ℒ₁𝑦_𝑟+1,𝑦⁽¹⁾_𝑟+1,…,𝑦^(𝑛)_𝑟+1=Φ𝑦_𝛾,𝑦_𝛾⁽¹⁾,…,𝑦_𝛾^(𝑛) (17) which, on using the definitions (6) and (7), can be written as 𝐿𝑦_𝑟+1,𝑦⁽¹⁾_𝑟+1,…,𝑦^(𝑛)_𝑟+1+^𝑛_𝑠=0𝑦^(𝑠)_𝑟+1−𝑦_𝑟^(𝑠)𝜕𝐹𝜕𝑦^(𝑠)𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)+𝐹𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)=0. (18)

We note that the iterative scheme (18) is, in fact, the quasilinearization method of Bellman and Kalaba [1]. For 𝑀=1, we have 𝑦(𝑥)≈𝑦₀(𝑥)+𝑦₁(𝑥), (19) where 𝑦₁ is obtained as a solution of ℒ₁𝑦₁,𝑦₁⁽¹⁾,…,𝑦₁^(𝑛)=ℏℒ₁𝑦₀,𝑦₀⁽¹⁾,…,𝑦₀^(𝑛)+ℏ𝑅₀𝑦₀,𝑦₀⁽¹⁾,…,𝑦₀^(𝑛). (20)

This produces the iteration scheme ℒ₁𝑦_𝑟+1,…,𝑦^(𝑛)_𝑟+1=Φ𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)+ℏℒ₁𝑦_0,𝑟+1,…,𝑦^(𝑛)_0,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1,…,𝑦^(𝑛)_0,𝑟+1, (21) where ℒ₁𝑦_𝑟+1,𝑦⁽¹⁾_𝑟+1,…,𝑦^(𝑛)_𝑟+1=𝐿𝑦_𝑟+1,𝑦⁽¹⁾_𝑟+1,…,𝑦^(𝑛)_𝑟+1+^𝑛_𝑠=0𝑦^(𝑠)_𝑟+1𝜕𝐹𝜕𝑦^(𝑠)𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛) (22) and 𝑦_0,𝑟+1 is the solution of 𝐿𝑦_0,𝑟+1,𝑦⁽¹⁾_0,𝑟+1,…,𝑦^(𝑛)_0,𝑟+1+^𝑛_𝑠=0𝑦^(𝑠)_0,𝑟+1𝜕𝐹𝜕𝑦^(𝑠)𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)=Φ𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛). (23)

For 𝑀=2, we have 𝑦(𝑥)≈𝑦₀(𝑥)+𝑦₁(𝑥)+𝑦₂(𝑥), (24) where 𝑦₂ is obtained as a solution of ℒ₁𝑦₂,𝑦₂⁽¹⁾,…,𝑦₂^(𝑛)=(1+ℏ)ℒ₁𝑦₁,𝑦₁⁽¹⁾,…,𝑦₁^(𝑛)+ℏ𝑅₁𝑦₀,𝑦₀⁽¹⁾,…,𝑦₀^(𝑛);𝑦₁,𝑦₁⁽¹⁾,…,𝑦₁^(𝑛). (25)

This produces the iteration scheme ℒ𝑦_𝑟+1=Φ𝑦_𝑟+ℏℒ₁𝑦_0,𝑟+1+(1+ℏ)ℒ₁𝑦_1,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1+𝑅₁𝑦_0,𝑟+1,𝑦_1,𝑟+1, (26) where 𝑦_1,𝑟+1 is obtained as the solution of ℒ₁𝑦_1,𝑟+1=ℏℒ₁𝑦_0,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1. (27)

In general, for any 𝑀>1, we have 𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥)+⋯+𝑦_𝑀(𝑥), (28) where 𝑦_𝑀(𝑥) is obtained as a solution of ℒ₁𝑦_𝑀=(1+ℏ)ℒ₁𝑦_𝑀−1+ℏ𝑅_𝑀−1𝑦₀,𝑦₁,…,𝑦_𝑀−1. (29)

Thus, a general scheme when the SHAM is truncated at order 𝑀 (where 𝑀≥1), hereinafter referred to as scheme-𝑀, can be obtained as ℒ₁𝑦_𝑟+1,𝑦⁽¹⁾_𝑟+1,…,𝑦^(𝑛)_𝑟+1=Φ𝑦_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)+^𝑀−1_𝑝=0𝜒_𝑝+1+ℏℒ₁𝑦_𝑝,𝑟+1,𝑦⁽¹⁾_𝑝,𝑟+1,…,𝑦^(𝑛)_𝑝,𝑟+1+ℏ^𝑀−1_𝑝=0𝑅_𝑝𝑦_0,𝑟+1,𝑦_1,𝑟+1,…,𝑦_𝑝,𝑟+1, (30) where each 𝑦_𝑝,𝑟+1 is obtained as the solution of ℒ₁𝑦_𝑝,𝑟+1=𝜒_𝑝+ℏℒ₁𝑦_{𝑝−1,𝑟+1}+ℏ𝑅_𝑝−1𝑦_0,𝑟+1,𝑦_1,𝑟+1,…,𝑦_{𝑝−1,𝑟+1}Φ𝑦when𝑝≥1,_𝑟,𝑦_𝑟⁽¹⁾,…,𝑦_𝑟^(𝑛)when𝑝=0. (31)

3. Solution of the Falkner-Skan Equation

In this section we demonstrate how the numerical schemes derived in the previous section may be used to solve the Falkner-Skan equation: 𝑓^′′′𝜂+𝛽𝑓𝜂𝑓^′′𝜂+𝛽₁1−𝑓^′𝜂²=0,𝜂∈[0,∞), (32) subject to the boundary conditions 𝑓(0)=𝑓^′(0)=0,lim_𝜂→∞𝑓𝜂=1. (33) It is convenient to first define 𝐹𝑦,𝑦^′,𝑦^′′=𝛽𝑦𝑦^′′−𝛽₁𝑦^′2,Ψ=−𝛽₁, (34) so that Φ𝑦,𝑦^′,𝑦^′′=𝛽𝑦𝑦^′′−𝛽₁𝑦^′2−𝛽₁,ℒ₁𝑦,𝑦^′,𝑦^′′,𝑦^′′′=𝑦^′′′+𝑎_0,𝛾𝑦^′′+𝑎_1,𝛾𝑦^′+𝑎_2,𝛾𝑅𝑦,_𝑚−1𝑦₀,𝑦₁,…,𝑦_𝑚−1=𝛽^𝑚−1_𝑗=0𝑦_𝑗𝑦^′′_{𝑚−1−𝑗}−𝛽₁^𝑚−1_𝑗=0𝑦^′_𝑗𝑦^′_{𝑚−1−𝑗}+𝛽₁1−𝜒_𝑚+2𝛽₁𝑦^′_𝛾𝑦^′_𝑚−1−𝛽𝑦_𝛾𝑦^′′_𝑚−1+𝑦_𝛾^′′𝑦_𝑚−1, (35)(36)(37) where 𝑎_0,𝛾=𝛽𝑦_𝛾,𝑎_1,𝛾=−2𝛽₁𝑦^′_𝛾,𝑎_2,𝛾=𝛽𝑦_𝛾^′′. (38)Using (35)–(37) the first three iterative schemes corresponding to 𝑀=0,1,2 may now be defined as follows.

3.1. Scheme-0

In this scheme we set ℒ₁𝑦_𝑟+1,𝑦^′_𝑟+1,𝑦^′′_𝑟+1,𝑦^′′′_𝑟+1=Φ𝑦_𝑟,𝑦^′_𝑟,𝑦_𝑟^′′ (39) subject to 𝑦_𝑟+1(0)=𝑦^′_𝑟+1(0)=0,𝑦^′_𝑟+1(∞)=1. (40) It is worth noting that Scheme-0 is, in fact, equivalent to the original QLM algorithm; see Mandelzweig and Tabakin [5] and Mandelzweig [6].

3.2. Scheme-1

For this scheme we set ℒ₁𝑦_𝑟+1,𝑦^′_𝑟+1,𝑦^′′_𝑟+1,𝑦^′′′_𝑟+1=Φ𝑦_𝑟,𝑦^′_𝑟,𝑦_𝑟^′′+ℏℒ₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1,𝑦^′′′_0,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1, (41) subject to 𝑦_𝑟+1(0)=𝑦^′_𝑟+1(0)=0,𝑦^′_𝑟+1(∞)=1, (42) where 𝑅₀𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1=𝛽𝑦_0,𝑟+1𝑦^′′_0,𝑟+1−𝛽₁𝑦^′2_0,𝑟+1+𝛽₁+2𝛽₁𝑦^′_𝑟𝑦^′_0,𝑟+1−𝛽𝑦_𝑟𝑦^′′_0,𝑟+1+𝑦_𝑟^′′𝑦_0,𝑟+1, (43) and 𝑦_0,𝑟+1 is the solution of ℒ₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1,𝑦^′′′_0,𝑟+1=Φ𝑦_𝑟,𝑦^′_𝑟,𝑦_𝑟^′′ (44) subject to 𝑦_0,𝑟+1(0)=𝑦^′_0,𝑟+1(0)=0,𝑦^′_0,𝑟+1(∞)=1. (45)

3.3. Scheme-2

The complexity of the defining equations increases with the order of the scheme. For Scheme-2 we have ℒ₁𝑦_𝑟+1,𝑦^′_𝑟+1,𝑦^′′_𝑟+1,𝑦^′′′_𝑟+1=Φ𝑦_𝑟,𝑦^′_𝑟,𝑦_𝑟^′′+ℏℒ₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1,𝑦^′′′_0,𝑟+1+(1+ℏ)ℒ₁𝑦_1,𝑟+1,𝑦^′_1,𝑟+1,𝑦^′′_1,𝑟+1,𝑦^′′′_1,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1+ℏ𝑅₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1;𝑦_1,𝑟+1,𝑦^′_1,𝑟+1,𝑦^′′_1,𝑟+1, (46) subject to 𝑦_𝑟+1(0)=𝑦^′_𝑟+1(0)=0,𝑦^′_𝑟+1(∞)=1, (47) where 𝑅₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1;𝑦_1,𝑟+1,𝑦^′_1,𝑟+1,𝑦^′′_1,𝑟+1=𝛽𝑦_0,𝑟+1𝑦^′′_1,𝑟+1+𝑦_1,𝑟+1𝑦^′′_0,𝑟+1−2𝛽₁𝑦^′_0,𝑟+1𝑦_1,𝑟+1+2𝛽₁𝑦^′_𝑟𝑦^′_1,𝑟+1−𝛽𝑦_𝑟𝑦^′′_1,𝑟+1+𝑦_𝑟^′′𝑦_1,𝑟+1, (48) and 𝑦_1,𝑟+1 is the solution of ℒ₁𝑦_1,𝑟+1,𝑦^′_1,𝑟+1,𝑦^′′_1,𝑟+1,𝑦^′′′_1,𝑟+1=ℏℒ₁𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1,𝑦^′′′_0,𝑟+1+ℏ𝑅₀𝑦_0,𝑟+1,𝑦^′_0,𝑟+1,𝑦^′′_0,𝑟+1, (49) with 𝑦_1,𝑟+1(0)=𝑦^′_1,𝑟+1(0)=0,𝑦^′_1,𝑟+1(∞)=0. (50) Equations (39), (41), and (46) describing the three solution schemes can be solved numerically using standard methods such as finite difference, finite elements, and spline collocation methods. In this study we use the Chebyshev spectral collocation method to solve the iteration schemes, (see [33–36]). To allow for numerical implementation of the pseudospectral method, the physical region [0,∞) is truncated to [0,𝐿] where 𝐿 is chosen to be sufficiently large. The truncated region is further transformed to the space [−1,1] using the transformation 2𝜉=𝐿_∞𝜂−1. (51) As with any other numerical approximation method, some sort of discretization is introduced in the interval [−1,1]. We choose the Gauss-Lobatto collocation points to define the nodes in [−1,1] as 𝜉_𝑗=cos𝜋𝑗𝑁,𝑗=0,1,…,𝑁, (52) where (𝑁+1) is the number of collocation points. The essence of the Chebyshev spectral collocation method is the idea of introducing a differentiation matrix 𝐷. The differentiation matrix maps a vector of the function values 𝐘=[𝑦(𝜉₀),…,𝑦(𝜉_𝑁)]^𝑇 at the collocation points to a vector 𝐘^′ defined as 𝐘^′=^𝑁_𝑘=0𝐷_𝑗𝑘𝑓𝜉_𝑘=𝐷𝐘, (53) In general, a derivative of order 𝑝 for the function 𝑦(𝜂) can be expressed as 𝑦^(𝑝)𝜂=𝐃^𝑝𝐘, (54) where 𝐃=2𝐷/𝐿_∞. The matrix 𝐷 is of size (𝑁+1)×(𝑁+1) and its entries are defined as 𝐷_𝑗𝑘=𝑐_𝑗𝑐_𝑘(−1)^𝑗+𝑘𝜏_𝑗−𝜏_𝑘𝐷,𝑗≠𝑘;𝑗,𝑘=0,1,…,𝑁,_𝑘𝑘𝜏=−^𝑘21−𝜏²_𝑘𝐷,𝑘=1,2,…,𝑁−1,₀₀=2𝑁²+16=−𝐷_𝑁𝑁, (55) with 𝑐_𝑘=2𝑘=0,𝑁1−1≤𝑘≤𝑁−1. (56) Thus, applying the spectral method to the iteration Scheme-0 (39) and the corresponding boundary conditions gives the following matrix system: 𝐀_𝑟𝐘_𝑟+1=Φ_𝑟 (57) with boundary conditions 𝑦_𝑟+1𝜉_𝑁=0,^𝑁_𝑘=0𝐃_𝑁𝑘𝑦_𝑟+1𝜉_𝑘=0,_𝑁𝑘=0𝐃_0𝑘𝑦_𝑟+1𝜉_𝑘=1, (58) where 𝐀_𝑟=𝐃³+𝐚_0,𝑟𝐃²+𝐚_1,𝑟𝐃+𝐚_2,𝑟, (59)where Φ_𝑟 corresponds to the function Φ(𝑦,𝑦,𝑦^′′) when evaluated at the collocation points and 𝐚_𝑖,𝑟  (𝑖=0,1,2) is a diagonal matrix corresponding to the vector of 𝑎_𝑖,𝑟.

The boundary conditions (58) are imposed on the first, 𝑁th, and (𝑁+1)th rows of 𝐴_𝑟 and Φ_𝑟 to obtain a system of the form 𝐃_0,0𝐃_0,1⋯𝐃_0,𝑁−1𝐃_0,𝑁𝐀_𝐫𝐃_𝑁,0𝐃_𝑁,1⋯𝐃_{𝑁,𝑁−1}𝐃_𝑁,𝑁𝑦00⋯01^𝑟+1𝜉₀𝑦_𝑟+1𝜉₁⋮𝑦_𝑟+1𝜉_𝑁−2𝑦_𝑟+1𝜉_𝑁−1𝑦_𝑟+1𝜉_𝑁=1Φ_𝑟𝜉₁⋮Φ_𝑟𝜉_𝑁−200. (60) Starting from a suitable initial guess 𝑦₀(𝜂), the iteration scheme (60) can be used to iteratively give approximate solutions of the governing equation (32) for Scheme-0. The application of the pseudospectral method for Scheme-1 and Scheme-2 can be done in a similar manner. The initial approximation used in all the algorithms is 𝑦₀𝜂=𝜂+𝑒^−𝜂+1. (61) The number of collocation points used in all the results presented here is 𝑁=200 with 𝐿_∞=20.

4. Results and Discussion

In this section we present solutions of the Falkner-Skan equation (32) using the QLM-SHAM hybrid iteration schemes. Numerical simulations were conducted for the following special classes of the F-S equations:(i)Blasius flow: 𝛽=1/2,𝛽₁=0, (ii)Pohlhausen flow: 𝛽=0,𝛽₁=1, (iii)Homann flow: 𝛽=2,𝛽₁=1.

To assess the accuracy and performance of our schemes, the numerical results were compared to the recently reported results of Ganapol [29]. To date, these results are the most accurate results for the Blasius and Falkner-Skan class of equations. Ganapol [29] reported highly accurate results between 10 and 30 decimal places using a robust algorithm based on Maclaurin series with convergence acceleration and analytical continuation techniques.

The comparison between the present findings and the results in the literature is made for the skin friction which is proportional to 𝑓^′′(0). Table 1 shows a comparison between the computed skin friction values of the Blasius equation using the three QLM-SHAM iteration schemes. The results are compared with the results reported in Ganapol [29] which are accurate to 29 decimal places. We observe that all the iteration schemes rapidly converge to the results of [29] to all 29 decimal places. Full convergence is achieved after 6 iterations when using Scheme-0, 4 iterations when using Scheme-1, and after 3 iterations when using Scheme-2. It is worth noting that the results of [29] were achieved after 104 decimal places. Prior to Ganapol [29], the most accurate Blasius skin friction results had been published to 17 decimal places by Boyd [28] as 𝑓^′′(0)=0.33205733621519630. This result was obtained after 5 iterations using Scheme-0 and 3 iterations for both Schemes-1 and -2. The value reported after 52 iterations in [29] is 𝑓^′′(0)=0.3320573362151965. It is clear that the proposed iteration schemes converge significantly faster than the method of [29]. That the results of Boyd [28] and Ganapol [29] were obtained only after a few iterations validates both the higher order convergence and the accuracy of the present solution methods.

Table 1 

Comparison between the computed values of the skin friction 𝑓^′′(0) for the Blasius flow (𝛽=1/2,𝛽₁=0).

Table 2 shows the computed values of 𝑓^′′(0) using Schemes-0 and -1 for the Pohlhausen flow (𝛽=0,𝛽₁=1). For this particular flow, the exact value of 𝑓^′′(0) is known to be 2/√3. The iteration Scheme-0 matches the exact result after only 5 iterations and Scheme-1 converges after only 3 iterations. The method used in Ganapol [29] converged to the exact result after 104 iterations. This again demonstrates the superior convergence of the present method.

Table 2 

Comparison between the computed values of the skin friction 𝑓^′′(0) for the Pohlhausen flow (𝛽=0,𝛽₁=1).

Table 3 gives the numerical simulations of the skin friction results for the Homann flow. We observe that the 29-digit results reported in [29] are achieved in 5 iterations, 4 iterations, and 3 iterations for Scheme-0, -1, and -2, respectively. This result indicates that adding an additional level in the QLM-SHAM scheme would further significantly increase the convergence of the iteration scheme. We further note from Tables 1–3 that all three schemes converge significantly much faster than the spectral homotopy analysis method on its own.

Table 3 

Comparison between the computed values of the skin friction 𝑓′′(0) for the Homann flow (𝛽=2,𝛽₁=1).

5. Conclusion

In this study we presented three hybrid QLM-SHAM iteration schemes for the solution of Falkner-Skan type boundary layer equations. We have shown through numerical experimentation that the proposed numerical schemes significantly enhance the convergence rate of the quasilinearization method. By comparison with the most accurate solutions of the Falkner-Skan equations currently available in the literature, we have shown that the schemes are highly accurate and efficient in terms of the number of iterations required to determine the solution to the required level of accuracy. The schemes presented provide robust tools for the efficient solution of nonlinear equations by offering superior accuracy to many existing methods. In addition, the approach used in deriving these schemes provides a suitable framework for extension to higher level schemes by adding more terms of the SHAM component of the method.