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

Authors: Abdelrady Okasha Elnady; M. Fayek Abd Rabbo; Hani M. Negm
Journal: Journal of Engineering
 (2020)

---

## Abstract

A numerical method for the solution of the Falkner–Skan equation, which is a nonlinear differential equation, is presented. The 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.


---

## Body

## 1. Introduction

Ordinary differential equations are important tools in solving real-world problems. Several natural phenomena are modelled by ODEs that have been used in several fields, such as physics, engineering, and biology [1–7].The Falkner–Skan equation arises in the study of laminar boundary layers exhibiting similarity. The similarity solutions of the two-dimensional incompressible laminar boundary layer equations are well-known as the Falkner–Skan equation. The F–S equation is a one-dimensional third-order nonlinear two-point boundary-value problem which has no closed-form solution. The problem is given by(1)f‴+β0ff″+β1−f′2=0,0≤η<∞,where β0 and β are constants and primes denote differentiation with respect to η. The associated boundary conditions are given by(2)f0=0,f′0=0,limη⟶∞f′η=1.The solution of (1) and (2) is characterized by f″0=α.The 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]. These techniques have mainly used shooting algorithms or invariant imbedding. The Chebshev collocation matrix method [18] has been presented the numerical solution of nonlinear differential equations. The 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. The 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. The solution of this system yields the Chebyshev coefficients of the solution function. The method is valid for both initial-value and boundary-value problems.The 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. The importance of the present method arises from its simplicity and the fact that it does not require to guess the value off″0.


## 2. Chebyshev Series

Any continuous functionf(ξ) in the interval 0 ≤ ξ ≤ 1 can be written in a Chebyshev series as follows [20]:(3)fξ=∑r=0∞a+Trξ.The series expansion (3) is fast-converging, and a good approximation is obtained by taking a few terms. Therefore, equation (3) is approximated by(4)fξ=∑r=0Na+rTrξ,where + sign means that the 1st term must be halved and ar are constants to be determined so as to obtain the best possible fit.(5)Trξ=Cosrt,Cost=2ξ‐1,0≤ξ≤1.The shifted Chebyshev polynomials satisfy the recurrence relations(6)Tr+1ξ=22ξ−1TrξTr−1ξ,2≤r≤N,T0=1,T1=2ξ‐1,and the orthogonality conditions(7)∫01TmξTnξξ1−ξdξ=0,form≠n,π2,form=n≠0,π,form=n=0,where for a known function fξ, the coefficients are given by(8)ar=2π∫01fξTrξξ1−ξdξ,0≤r≤N.Also, the derivatives of the functionfξ can be expanded in a Chebyshev series keeping a recurrence relation between the coefficients of the function and its derivatives. The first derivative f′ξ is expressed in a Chebyshev series as [20](9)f′ξ=∑r=0N−1a+r1Trξ.The coefficientsar1 satisfy the recursive relation(10)ar−11−ar+11=4rar,1≤r≤N.Similarly, themth derivative is written as(11)fmξ=∑r=0N−ma+rmTrξ,where ar−1m−ar+1m=4rarm−1, 1 ≤ r ≤ N-(m - 1).


### 2.1. Matrix Representation of Function and Function Derivatives

Any continuous functionf(ξ) in the interval 0 ≤ ξ ≤ 1 and its derivatives can be written in a matrix form of the Chebyshev series as follows:(12)fξ=TrIai,r=0,1,2,…,N,i=0,1,2,…,N.The first-order-derivative coefficientsar1 in equation (10) can be written in terms of the original function coefficients {ai} using matrix notation as follows:(13)ar1=4Aai,r=0,1,2,…,N−1,i=1,2,3,…,N,where [A] is an upper triangular matrix of order N × N.The elements of the matrixaij are defined as(14)aij=0,i>jori+jodd,j,i≤jandi+jeven.The form of matrix [A] for N = 5 for example is(15)A=1030502040003050004000005.From equation (13), it is noted that the first-order-derivative coefficients are written in terms of the N coefficients {ai}, (i = 1, 2, …, N) of the function f(ξ). To represent ar1 in terms of all function coefficients {ai}, i = 0, 1, 2, …, N, we add a new left column with zero entries in the matrix [A], and the new matrix is termed [A01]. Thus,(16)ar1=4A01ai,r=0,1,2,…,N−1,i=1,2,3,…,N,where [A01] is of order N × N + 1 and takes; for N = 5, the form(17)A01=010305002040000305000040000005.The second-order-derivative coefficientsar2 can be written in terms of the function coefficients {ai} using matrix notation as follows:(18)ar2=16A−1,−1A01ai,r=0,1,2,…,N−2,ar2=16A02ai,i=0,1,2,…,N,where [A02] is of order N − 1 × N + 1. [A]−1,−1 is the matrix [A] after deleting the last row and last column. −1[A01] is the matrix [A01] after deleting the first row.(19)A02=1030020400300004∗002040000305000040000005=0020160000603000001200000020.The general form of themth derivative coefficients arm can be written in terms of the function coefficients {ai} using matrix notation as follows:(20)arm=4mA1−m,1−mA0m−1−1ai,r=0,1,2,…,N−m,arm=4mA0mai;i=0,1,2,…,N,where m the order of derivative, [A]1−m, 1−m matrix [A] after deleting the last (m − 1) rows and (m − 1) columns. −1[ ] matrix [ ] after deleting the first row.(21)fmξ=4mTrA0mai,r=0,1,2,…,N−m,i=0,1,2,…,N,where [Tr] is a row matrix whose elements are Tr(ξ). Note that in the matrix form the first term of [Tr] must be halved.Now, consider the general nonhomogenous differential equation of mth order.(22)fm+fm−1+fm−2+⋯+f′+f=P.After expanding each term in the Chebyshev series, the abovementioned differential equation can be written as(23)∑r=0N−ma+rmTrξ+∑r=0N−m−1a+rm−1Trξ+⋯+∑r=0N−1a+r1Trξ+∑r=0Na+rTrξ=∑r=0Np+rTrξ.The forcing-function coefficientspr 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(24)4mA0m+4m−1A0m−1+⋯+4A01+Iai=pr;r=0,1,2,…,N+1−m,i=0,1,2,…,N,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. These additional equations are supplied by the problem boundary conditions.


### 2.2. Matrix Representation of Function Products

Iff(ξ) and g(ξ) are two continuous functions represented by truncated Chebyshev series as(25)fξ=∑r=0Na+rTrξ,gξ=∑r=0Mb+rTrξ,then the product of these functions can be written in a Chebyshev series as(26)gξfξ=∑r=0M+Nc+rTrξ,where(27)co=∑i=0N+Ma+ibi,cr=12∑i=0N+Ma+ibi+r+bi−r,1≤r≤M+N,i+r≤M,i−r≤M∀i.The {cr} coefficients can be written in terms of the {ai} coefficients only using matrix notation as follows:(28)cr=Hai,where r = 0, 1, 2, …, N + M, i = 0, 1, 2, …, N, {cr} is a column matrix of order N + M + 1 × 1, [H] is rectangular matrix of order N + M + 1 × N + 1, and {ai} is a column matrix of order N + 1 × 1.The elementshij of the matrix [H] can be written as(29)hij=14bi−j+bi+j−2,fori=1,2,…,N+M+1,j=1,12bi−j+bi+j−2,fori=1,2,…,N+M+1,j=2,3,…,N+1,i−j≤M,i+j−2≤M,∀i,j.Therefore, equation (28) takes the form(30)c0c1c2c3⋅⋅⋅cN+M=12b02b12b22b3⋅⋅⋅b1b0+b2b1+b3b2+b4⋅⋅⋅b2b1+b3b0+b4b1+b5⋅⋅⋅b3b2+b4b1+b5b0+b6⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅a0a1a2a3⋅⋅⋅aN.


## 2.1. Matrix Representation of Function and Function Derivatives

Any continuous functionf(ξ) in the interval 0 ≤ ξ ≤ 1 and its derivatives can be written in a matrix form of the Chebyshev series as follows:(12)fξ=TrIai,r=0,1,2,…,N,i=0,1,2,…,N.The first-order-derivative coefficientsar1 in equation (10) can be written in terms of the original function coefficients {ai} using matrix notation as follows:(13)ar1=4Aai,r=0,1,2,…,N−1,i=1,2,3,…,N,where [A] is an upper triangular matrix of order N × N.The elements of the matrixaij are defined as(14)aij=0,i>jori+jodd,j,i≤jandi+jeven.The form of matrix [A] for N = 5 for example is(15)A=1030502040003050004000005.From equation (13), it is noted that the first-order-derivative coefficients are written in terms of the N coefficients {ai}, (i = 1, 2, …, N) of the function f(ξ). To represent ar1 in terms of all function coefficients {ai}, i = 0, 1, 2, …, N, we add a new left column with zero entries in the matrix [A], and the new matrix is termed [A01]. Thus,(16)ar1=4A01ai,r=0,1,2,…,N−1,i=1,2,3,…,N,where [A01] is of order N × N + 1 and takes; for N = 5, the form(17)A01=010305002040000305000040000005.The second-order-derivative coefficientsar2 can be written in terms of the function coefficients {ai} using matrix notation as follows:(18)ar2=16A−1,−1A01ai,r=0,1,2,…,N−2,ar2=16A02ai,i=0,1,2,…,N,where [A02] is of order N − 1 × N + 1. [A]−1,−1 is the matrix [A] after deleting the last row and last column. −1[A01] is the matrix [A01] after deleting the first row.(19)A02=1030020400300004∗002040000305000040000005=0020160000603000001200000020.The general form of themth derivative coefficients arm can be written in terms of the function coefficients {ai} using matrix notation as follows:(20)arm=4mA1−m,1−mA0m−1−1ai,r=0,1,2,…,N−m,arm=4mA0mai;i=0,1,2,…,N,where m the order of derivative, [A]1−m, 1−m matrix [A] after deleting the last (m − 1) rows and (m − 1) columns. −1[ ] matrix [ ] after deleting the first row.(21)fmξ=4mTrA0mai,r=0,1,2,…,N−m,i=0,1,2,…,N,where [Tr] is a row matrix whose elements are Tr(ξ). Note that in the matrix form the first term of [Tr] must be halved.Now, consider the general nonhomogenous differential equation of mth order.(22)fm+fm−1+fm−2+⋯+f′+f=P.After expanding each term in the Chebyshev series, the abovementioned differential equation can be written as(23)∑r=0N−ma+rmTrξ+∑r=0N−m−1a+rm−1Trξ+⋯+∑r=0N−1a+r1Trξ+∑r=0Na+rTrξ=∑r=0Np+rTrξ.The forcing-function coefficientspr 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(24)4mA0m+4m−1A0m−1+⋯+4A01+Iai=pr;r=0,1,2,…,N+1−m,i=0,1,2,…,N,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. These additional equations are supplied by the problem boundary conditions.


## 2.2. Matrix Representation of Function Products

Iff(ξ) and g(ξ) are two continuous functions represented by truncated Chebyshev series as(25)fξ=∑r=0Na+rTrξ,gξ=∑r=0Mb+rTrξ,then the product of these functions can be written in a Chebyshev series as(26)gξfξ=∑r=0M+Nc+rTrξ,where(27)co=∑i=0N+Ma+ibi,cr=12∑i=0N+Ma+ibi+r+bi−r,1≤r≤M+N,i+r≤M,i−r≤M∀i.The {cr} coefficients can be written in terms of the {ai} coefficients only using matrix notation as follows:(28)cr=Hai,where r = 0, 1, 2, …, N + M, i = 0, 1, 2, …, N, {cr} is a column matrix of order N + M + 1 × 1, [H] is rectangular matrix of order N + M + 1 × N + 1, and {ai} is a column matrix of order N + 1 × 1.The elementshij of the matrix [H] can be written as(29)hij=14bi−j+bi+j−2,fori=1,2,…,N+M+1,j=1,12bi−j+bi+j−2,fori=1,2,…,N+M+1,j=2,3,…,N+1,i−j≤M,i+j−2≤M,∀i,j.Therefore, equation (28) takes the form(30)c0c1c2c3⋅⋅⋅cN+M=12b02b12b22b3⋅⋅⋅b1b0+b2b1+b3b2+b4⋅⋅⋅b2b1+b3b0+b4b1+b5⋅⋅⋅b3b2+b4b1+b5b0+b6⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅a0a1a2a3⋅⋅⋅aN.


## 3. Application to the Falkner–Skan Equation

The purpose of the present section is to solve the Falkner–Skan equations (1) and (2) numerically by using the Chebyshev series in a matrix form. Again, to solve this problem, the semi-infinite physical domain [0, ∞ of the problem is truncated to a finite domain [0, L], where L is a sufficiently large number. The Falkner–Skan equation (1) can be written in the nondimensional form as(31)1L3d3fξdξ3+β0fξ1L2d2fξdξ2−β1Ldfξdξ1Ldfξdξ=β,ξ=ηL.The associated boundary conditions (2) are written in the nondimensional form as(32)f0=0,1Ldf0dξ=0,1Ldf1dξ=1.Expandingf(ξ) 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:(33)64L3A03+β016L2Hb−β4LHa4LA01fi=βr,where [Ha] is the matrix of coefficients of the first derivative, [Hb] is the matrix of coefficients of the second derivative, {fi} are Chebyshev coefficients of the function f, and {βr} are Chebyshev coefficients of β and obtained from equation (8).The 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).


### 3.1. Boundary Conditions

Fortunately, it is relatively easy to represent any boundary conditions for the functions expanded in the Chebyshev series. The boundary conditions given in (32) can be written as(34)TR0fi=0,4LTR01A01fi=4LT0A01fi=0,4LTR11A01fi=4LT1A01fi=1,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.


### 3.2. Procedure of Solution

The present algorithm consists of the following procedure:(a)
It is started with a relatively small value ofL as the initial value for semi-infinite domain.(b)
The system of equations (33) and (34) are solved simultaneously to find the velocity profile f′η.(c)
L is continuously increased up to limη⟶∞f′η=1. This is achieved if two successive values f″0 differ by about 1E − 6. Figure 1 shows the results for β0=1 and β=2.Figure 1
Velocity profile for (β0 = 1 and β = 2).


## 3.1. Boundary Conditions

Fortunately, it is relatively easy to represent any boundary conditions for the functions expanded in the Chebyshev series. The boundary conditions given in (32) can be written as(34)TR0fi=0,4LTR01A01fi=4LT0A01fi=0,4LTR11A01fi=4LT1A01fi=1,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.


## 3.2. Procedure of Solution

The present algorithm consists of the following procedure:(a)
It is started with a relatively small value ofL as the initial value for semi-infinite domain.(b)
The system of equations (33) and (34) are solved simultaneously to find the velocity profile f′η.(c)
L is continuously increased up to limη⟶∞f′η=1. This is achieved if two successive values f″0 differ by about 1E − 6. Figure 1 shows the results for β0=1 and β=2.Figure 1
Velocity profile for (β0 = 1 and β = 2).


## 4. Results and Discussion

The Falkner–Skan equation has two coefficientsβ0 and β. The solutions corresponding to β = 0 have been known as constant flow, those corresponding to β > 0 are called accelerating flow, and those corresponding to β < 0 are known as decelerating flow. Physically relevant solutions exist only for −0.1988 < β ≤ 2. If β0 = 1/2 and β = 0, it is called Blasius flow; if β0 = 1 and β = ½, it describes the Homann axisymmetric stagnation flow, if β0 = 1 and β = 1, it describes the Hiemenz flow, if β0 = 0 and β = 1, it is called the Pohlhausen flow, and if β0 = 2 and β = 1, it represents the problem of Homann, describing the steady flow in the boundary layer along a surface of revolution near the stagnation point. Sometimes, the Falkner–Skan equation specifically refers to β0 = 1.The 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.Table 1
Comparison off″0 for β0 = 1 and different values of β.
β
f
″
0 [17]
Present method
Absolute error
L
40
7.314784974
7.314784949
2.5E − 08
1.19
30
6.338208628
6.338208608
2.0E − 08
1.35
20
5.180718025
5.180718025
3.9E − 10
1.61
15
4.491486898
4.491486992
−9.4E − 08
1.78
10
3.675234101
3.675234108
−6.8E − 09
2.18
2
1.687218169
1.687218169
3.1E − 11
4.37
1
1.232587657
1.232588439
−7.8E − 07
4.30
0.5
0.92768004
0.927686117
−6.1E − 06
4.41
0
0.469599988
0.4704412
−8.4E − 04
4.03It is noted thatL has a great effect on the solution. The 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.Table 2
Effect ofη on f″0 for β = 40 and β0 = 1 (error = 1E − 7).
η
1.01
1.02
1.03
1.04
1.05
1.06
f
″
0
7.3147854
7.3147853
7.3147852
7.3147852
7.3147851
7.3147851
η
1.07
1.08
1.09
1.1
1.11
1.12
f
″
0
7.3147851
7.3147850
7.3147850
7.3147849
7.3147849
7.3147849
η
1.13
1.14
1.15
1.16
1.19
f
″
0
7.3147848
7.3147848
7.3147847
7.3147846
7.3147857The velocity profile for different flows mentioned above is obtained using the present technique and Runge–Kutta 4th order method and presented in the following figures.Figure2 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 4th order method. The 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 2
Velocity profile of the Blasius flow (β0 = ½ and β = 0). f″0 = 0.332057336151 for L = 5.7205339053.Figure 3
Velocity profile of the Homann axisymmetric stagnation flow (β0 = 1 and β = ½). f″0 = 0.927686124423 for L = 4.404.Figure 4
Velocity profile of the Hiemenz flow (β0 = 1 and β = 1). f″0 = 1.232588438646 for L = 4.3.Figure5 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 4th 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. This note is valid for all values of β and β0 = 0.Figure 5
Velocity profile of the Pohlhausen flow (β0 = 0 and β = 1). f″0 = 1.154700538379 for L = 6.98476.Figure6 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 4th order method.Figure 6
Velocity profile of the problem of Homann, describing steady flow in the boundary layer along a surface of revolution near the stagnation point (β0 = 2 and β = 1). f″0 = 1.311937693752 for L = 3.18407676.The values off″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.Figure 7
Velocity profile forβ=−0.1.Figure 8
Velocity profile forβ=−0.15.The step of increasingL 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.Figure 9
Velocity profile forβ=−0.1, ΔL = 0.001; f″0 =0.326218125788279 (Forward flow), L = 3.913; f″0 = −0.140741363248159 (reverse flow), L = 8.Figure 10
Velocity profile forβ=−0.15, ΔL = 0.001; f″0 = 0.233877086171824 (Forward flow), L = 3.88; f″0 = −0.133696366858524 (reverse flow), L = 7.135.


## 5. 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. The nonlinear differential equation has been transformed into a system of algebraic equations which was presented in a matrix form. This 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. The 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.This technique can be considered as a powerful tool to solve the linear and nonlinear differential equations defined on a finite range.



---

*Source: 3972573-2020-05-01.xml*