Solution of the Falkner–Skan Equation Using the Chebyshev Series in Matrix Form

A numerical method for the solution of the Falkner–Skan equation, which is a nonlinear differential equation, is presented. )e method has been derived by truncating the semi-infinite domain of the problem to a finite domain and then expanding the required approximate solution as the elements of the Chebyshev series. Using matrix representation of a function and their derivatives, the problem is reduced to a system of algebraic equations in a simple way. From the computational point of view, the results are in excellent agreement with those presented in published works.

e Falkner-Skan equation arises in the study of laminar boundary layers exhibiting similarity. e similarity solutions of the two-dimensional incompressible laminar boundary layer equations are well-known as the Falkner-Skan equation.
e F-S equation is a one-dimensional third-order nonlinear two-point boundaryvalue problem which has no closed-form solution. e problem is given by where β 0 and β are constants and primes denote differentiation with respect to η. e associated boundary conditions are given by e solution of (1) and (2) is characterized by f ″ (0) � α. e numerical treatment of this problem was addressed by many authors, namely, Lakestani [8], Parand et al. [9], El-Nady and Abd Rabbo [10,11], Cebeci and Keller [12], Na [13], Asaithambi [14], Asaithambi [15], Elgazery [16], and Ganapol [17]. ese techniques have mainly used shooting algorithms or invariant imbedding. e Chebshev collocation matrix method [18] has been presented the numerical solution of nonlinear differential equations. e method in [18] transforms the nonlinear differential equation into the matrix equation, which corresponds to a system of nonlinear algebraic equations with unknown Chebyshev coefficients, via Chebyshev collocation points. e method in [19] transforms the nonlinear differential equation into the system of nonlinear algebraic equations with unknown shifted Chebyshev coefficients, via Chebyshev-Gauss collocation points. e solution of this system yields the Chebyshev coefficients of the solution function. e method is valid for both initial-value and boundary-value problems. e purpose of this paper is to develop an efficient method based on the Chebyshev series which is much more straightforward and simpler than the other existing algorithms. Using matrix notation of the Chebyshev series, the nonlinear differential equation converts into a system of algebraic equations which can be easily solved. e importance of the present method arises from its simplicity and the fact that it does not require to guess the value of f ″ (0).

Chebyshev Series
Any continuous function f(ξ) in the interval 0 ≤ ξ ≤ 1 can be written in a Chebyshev series as follows [20]: e series expansion (3) is fast-converging, and a good approximation is obtained by taking a few terms. erefore, equation (3) is approximated by where + sign means that the 1 st term must be halved and a r are constants to be determined so as to obtain the best possible fit.
T r (ξ) � Cos r t , e shifted Chebyshev polynomials satisfy the recurrence relations and the orthogonality conditions where for a known function f(ξ), the coefficients are given by Also, the derivatives of the function f(ξ) can be expanded in a Chebyshev series keeping a recurrence relation between the coefficients of the function and its derivatives.
e first derivative f ′ (ξ) is expressed in a Chebyshev series as [20] f ′ (ξ) � e coefficients a (1) r satisfy the recursive relation Similarly, the m th derivative is written as where a (m)

Matrix Representation of Function and Function
Derivatives. Any continuous function f(ξ) in the interval 0 ≤ ξ ≤ 1 and its derivatives can be written in a matrix form of the Chebyshev series as follows: e first-order-derivative coefficients a (1) r in equation (10) can be written in terms of the original function coefficients {a i } using matrix notation as follows: From equation (13), it is noted that the first-order-derivative coefficients are written in terms of the N coefficients e second-order-derivative coefficients a (2) r can be written in terms of the function coefficients {a i } using matrix notation as follows: after deleting the last row and last column. − 1 [A01] is the matrix [A01] after deleting the first row.
[A02] � 1 0 3 0 e general form of the m th derivative coefficients a (m) r can be written in terms of the function coefficients {a i } using matrix notation as follows: where m the order of derivative, where [T r ] is a row matrix whose elements are T r (ξ). Note that in the matrix form the first term of [T r ] must be halved. Now, consider the general nonhomogenous differential equation of m th order.
After expanding each term in the Chebyshev series, the abovementioned differential equation can be written as (23) e forcing-function coefficients p r can be evaluated using equation (8). Equating the coefficients of like Chebyshev polynomial terms on either side, the resulting N + 1 − m algebraic equations can be written in a matrix form using equations (16) and (21) as where all matrices in equation (24) are of the same order (N + 1 − m × N + 1). For all derivatives lower than the highest derivative, the first N + 1 − m rows are chosen so as to satisfy equation (24). In order to be able to solve equation (24), m additional equations are needed. ese additional equations are supplied by the problem boundary conditions.

Matrix Representation of Function
Products. If f(ξ) and g(ξ) are two continuous functions represented by truncated Chebyshev series as then the product of these functions can be written in a Chebyshev series as where e {c r } coefficients can be written in terms of the {a i } coefficients only using matrix notation as follows: erefore, equation (28) takes the form e Falkner-Skan equation (1) can be written in the nondimensional form as

Application to the Falkner-Skan Equation
e associated boundary conditions (2) are written in the nondimensional form as Expanding f(ξ) in (N + 1)-Chebyshev terms, we get N + 1 unknown coefficients. Now, using matrix notation for the functions, function derivatives, and functions products and applying the rule of matrix multiplication, equation (31) can be written as a system of algebraic equations in the following matrix form: where [Ha] is the matrix of coefficients of the first derivative, [Hb] is the matrix of coefficients of the second derivative, {f i } are Chebyshev coefficients of the function f, and {β r } are Chebyshev coefficients of β and obtained from equation (8). e highest derivative in equation (33) is of order 3, so the number of algebraic equations is N − 2 along with 3 boundary conditions at ξ � 0 and ξ � 1, leading to N + 1 equations in N + 1 unknowns, which can be easily solved. It is to be noted that all matrices in equation (33) are of order (N − 2 × N + 1).

Boundary Conditions.
Fortunately, it is relatively easy to represent any boundary conditions for the functions expanded in the Chebyshev series. e boundary conditions given in (32) can be written as where [TR0] is a row matrix of N + 1 Chebeychev terms at ξ � 0, [TR01] is a row matrix of N Chebeychev terms at ξ � 0, and [TR11] is a row matrix of N Chebeychev terms at ξ � 1. Note: the first term of Chebyshev terms must be halved.

Procedure of Solution.
e present algorithm consists of the following procedure: (a) It is started with a relatively small value of L as the initial value for semi-infinite domain.
(b) e system of equations (33) and (34) are solved simultaneously to find the velocity profile f ′ (η). Figure 1 shows the results for β 0 � 1 and β � 2. e Falkner-Skan equation was solved using the proposed Chebyshev series in a matrix form for different values of β 0 and β. It was noted that the solution is stable and convergent even for small number of terms less than ten. Finally, in all solutions, the number of Chebyshev terms is taken N � 20. Table 1 compares the initial slope f ″ (0) obtained using the present method with the values obtained in [17] by using highly accurate algorithm based on a Maclaurin series representation. It is seen that the present method delivers the same order of accuracy.

Results and Discussion
It is noted that L has a great effect on the solution. e results show that the system equations (33) and (34) have a solution with high level of accuracy for different values of the semi-infinite interval L as presented in Table 2. e velocity profile for different flows mentioned above is obtained using the present technique and Runge-Kutta 4 th order method and presented in the following figures. Figure 2 shows the velocity profile of the Blasius flow for β 0 � ½ and β � 0. Figure 3 shows the velocity profile of the Homann axisymmetric stagnation flow for β 0 � 1 and β � ½. Figure 4 shows the velocity profile of the Hiemenz flow for β 0 � 1 and β � 1. It is clear that the present solution has an excellent agreement with that obtained by the Runge-Kutta 4 th order method. e solutions shown are not a single solution, but it is a set of solutions for different lengths and gives the same value of f ″ (0) as presented in Table 2. Figure 5 shows the velocity profile of the Pohlhausen flow for β 0 � 0 and β � 1. It is clear in this figure that the solution obtained by the Runge-Kutta 4 th order method is divergent at L nearly equal 3.25 while the present solution is still convergent until L � 6.98476 which means that the present solution is powerful. is note is valid for all values of β and β 0 � 0. Figure 6 shows the velocity profile of the problem of Homann describing steady flow in the boundary layer along a surface of revolution near the stagnation point in which β 0 � 2 and β � 1. As mentioned previously, the present solution has an excellent agreement with that obtained by the Runge-Kutta 4 th order method.
e values of f ″ (0) and the largest L for each case are presented.
For β 0 � 1 and β < 0, two families of solutions are obtained. One family of solutions corresponds to forward flow and the other corresponds to reverse flow as shown in Figures 7 and 8. e step of increasing L is 0.05, and the error between two consecutive values of f ″ (0) is 1E − 2.
Refining ΔL � 0.001 and adopting the error between two consecutive values of f ″ (0) to be 1E − 6, the results are shown in Figures 9 and 10.

Conclusions
In the present study, a technique based on the Chebyshev series expansion is applied to determine an approximate solution for the Falkner-Skan (F-S) third-order nonlinear differential equation. e nonlinear differential equation has been transformed into a system of algebraic equations which was presented in a matrix form. is system of equations is easily solved by matrix inversion, by assuming initial values for chebyshev coefficients and iterating this process until the correct values are within acceptable error. e new proposed technique is straightforward and well adapted to the computer implementation. One of the advantages of this method is that the solution is expressed as a system of algebraic equations which is directly solved using the computer program without any computational effort. A comparison of the present results with the published results indicates excellent accuracy of the solution. It is concluded that the present technique is an accurate tool in handling the Falkner-Skan (F-S) equation with a high level of accuracy.
is technique can be considered as a powerful tool to solve the linear and nonlinear differential equations defined on a finite range.

Data Availability
All data and analysis are allowed to be shared by researchers throughout the published paper.

Conflicts of Interest
e authors declare that they have no conflicts of interest.