- Email: [email protected]
Numerical solution of the Falkner–Skan equation using piecewise linear functions
Applied Mathematics and Computation 159 (2004) 267–273 www.elsevier.com/locate/amc
Numerical solution of the Falkner–Skan equation using piecewise linear functions Asai Asaithambi Department of Computer Science, Parks College of Engineering and Aviation, Saint Louis University, 3450 Lindell Blvd., Saint Louis, MO 63040/63103, USA
Abstract We present a finite-element method for the solution of the Falkner–Skan equation. The method uses a coordinate transformation to map the semi-infinite domain of the problem to the unit interval ½0; 1. By means of a suitable change of variables, the transformed third-order boundary value problem is further rewritten as a system of differential equations consisting of a second-order equation and a first-order equation. The second-order equation is approximated using a Galerkin formulation with piecewise linear elements, and the first-order equation is approximated using a centereddifference approximation. In the process of the initial transformation, a finite computational boundary is introduced. This boundary is obtained as part of the computational solution by imposing an additional asymptotic boundary condition. The solutions obtained thus are in excellent agreement with those obtained by previous authors. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Falkner–Skan equation; Coordinate transformation; Galerkin formulation; Piecewise linear elements; Centered-difference; Asymptotic boundary condition
1. Introduction The Falkner–Skan equation arises in the study of laminar boundary layers exhibiting similarity. The solutions of the one-dimensional third-order
E-mail address: [email protected] (A. Asaithambi). 0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.10.047
268
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
boundary-value problem described by the well-known Falkner–Skan equation are the similarity solutions of the two-dimensional incompressible laminar boundary layer equations. This is a nonlinear two-point boundary value problem for which no closed-form solutions are available. The problem is described as follows: " 2 # d3 f d2 f df þf 2 þb 1 ¼ 0; 0 < g < 1; ð1Þ dg3 dg dg along with the boundary conditions f ¼0
at g ¼ 0;
ð2Þ
df ¼ 0 at g ¼ 0; dg
ð3Þ
df ¼ 1 as g ! 1: dg
ð4Þ
In addition to the unknown function f , the solution of (1)–(4) is characterized by the value of a ¼ d2 f =dg2 at g ¼ 0. The condition (4) is usually replaced by the condition df ¼ 1 at g ¼ g1 ; dg
ð5Þ
for some ‘‘sufficiently large’’ g1 , which must be determined as part of the solution. The addition of the new unknown g1 to the problem warrants the introduction of the ‘‘asymptotic condition’’ d2 f ¼ 0 at g ¼ g1 : dg2
ð6Þ
Rosenhead [1] and Weyl [2] present mathematical treatments of this problem focusing on existence and uniqueness results. The first numerical treatment of this problem was presented by Hartree [3]. Other numerical treatments included those developed by Smith [4], Cebeci and Keller [5], and Na [6]. These previous approaches have mainly used shooting and invariant imbedding. It is also important to note that the methods of Cebeci and Keller [5] have experienced convergence difficulties, which were overcome by moving towards more complicated methods. Recent methods presented by Asaithambi [7,8] have improved the performance of the previous methods greatly by reducing the amount of computational effort while obtaining the same results as the previous authors.
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
269
2. Method of solution With the change of variable g ¼ f =g1 , and under the transformation n ¼ g=g1 , the Falkner–Skan equation (1) is transformed to 2
g000 þ g21 gg00 þ g21 b½1 ðg0 Þ ¼ 0;
ð7Þ
where the primes denote differentiation with respect to the variable n. The boundary conditions 2, 3 and 5 are transformed to g0 ð0Þ ¼ 0;
gð0Þ ¼ 0;
g0 ð1Þ ¼ 1:
ð8Þ
The asymptotic boundary condition (6) is transformed to g00 ð1Þ ¼ 0:
ð9Þ
Direct numerical treatment of (7)–(9) is somewhat inconvenient since we have a third-order two-point boundary value problem. In order to overcome this difficulty, we let u ¼ g0 ;
ð10Þ
so that (7) transforms to u00 þ g21 gu0 þ g21 bð1 u2 Þ ¼ 0:
ð11Þ
The conditions in (8) are transformed to gð0Þ ¼ 0;
uð0Þ ¼ 0;
uð1Þ ¼ 1
ð12Þ
and the asymptotic boundary condition (9) is transformed to u0 ð1Þ ¼ 0:
ð13Þ
We wish to determine g and u for n 2 ½0; 1. We divide the interval ½0; 1 into N subintervals of length h ¼ 1=N each, and define nj ¼ jh for j ¼ 0; 1; . . . ; N . Then, we let gj and uj denote the values of gðnj Þ and uðnj Þ, respectively, for j ¼ 0; 1; . . . ; N . The boundary conditions (12) will now correspond to g0 ¼ 0;
u0 ¼ 0;
uN ¼ 1:
ð14Þ
The method of this paper is based on a finite-element formulation of (11). The method will use standard finite-difference approximations for the remaining parts of the problem. For the purpose of the finite-element formulation, we represent uðnÞ and gðnÞ in the forms uðnÞ ¼
n X j¼0
uj /j ðnÞ;
gðnÞ ¼
n X j¼0
gj /j ðnÞ;
ð15Þ
270
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
where uj uðnj Þ and gj 8 nn j1 > < nj nj1 /j ðnÞ ¼ njþ1 n > : njþ1 nj 0;
gðnj Þ as indicated previously, and for n 2 ½nj1 ; nj ; ð16Þ
for n 2 ½nj ; njþ1 ; otherwise:
Also, we represent (11) in its weak, or Galerkin form, as Z 1 ½u00 þ g21 gu0 þ g21 bð1 u2 Þ/j ðnÞ dn ¼ 0
ð17Þ
0
for j ¼ 1; 2; . . . ; N 1. For the first two integrands in (17), we use integration by parts and the fact that /j ð0Þ ¼ /j ð1Þ ¼ 0 for i ¼ 1; 2; . . . ; N 1, and obtain the equivalent form Z 1 Z 1 Z 1 u0 /0j dn g21 ug/0j dn g21 ug0 /j dn 0
þ
g21 b
Z
0
0
1 2
ð1 u Þ/j dn ¼ 0:
ð18Þ
0
Using (15) and (16) in (18), performing the integrations, collecting terms, and simplifying yields 1 uj1 2uj þ ujþ1 þ hg21 Aj þ bh2 g21 Bj ¼ 0; 6
ð19Þ
for i ¼ 1; 2; . . . ; N 1 with Aj ¼ gj1 ðuj uj1 Þ þ 2gj ðujþ1 uj1 Þ þ gjþ1 ðujþ1 uj Þ;
ð20Þ
Bj ¼ 1 121 ðu2j1 þ 2uj1 uj þ 6u2j þ 2uj ujþ1 þ u2jþ1 Þ:
ð21Þ
Next, we use the discretization gj gj1 h uj þ uj1 i ¼ 0; h 2
ð22Þ
corresponding to Eq. (10) applied at n ¼ nj1=2 for j ¼ 1; 2; . . . ; N . Note that the ð2N 1Þ equations described by (19) and (22) form a nonN N linear system of algebraic equations in the ð2N þ 3Þ variables fgj gj¼0 , fuj gj¼0 , and g1 . Among these ð2N þ 3Þ variables, the three quantities g0 ; u0 , and uN are available from (14). Assuming that g1 is available, (19) and (22) will be used to N N 1 determine the ð2N 1Þ unknown quantities fgj gj¼1 , and fuj gj¼1 . Then, a numerical representation of the asymptotic boundary condition (13) will be used to obtain an improved value for g1 . We now proceed to describe the iterative process for the solution of the nonlinear system (19) and (22). We begin by letting
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
uT ¼ ½ u1
uN ;
gT ¼ ½ g1
gN
271
ð23Þ
and 2
3 H1 ðu; g; g1 Þ 6 7 .. Hðu; g; g1 Þ ¼ 4 5; . HN 1 ðu; g; g1 Þ
ð24Þ
where Hj ðu; g; g1 Þ ¼ ujþ1 2uj þ uj1 þ 16hg21 Aj þ bh2 g21 Bj
ð25Þ
for j ¼ 1; 2; . . . ; N 1. If we consider (19) as a nonlinear system in u alone, then solving (19) is equivalent to solving the system described by Hðu; g; g1 Þ ¼ 0: Note that the Jacobian matrix JH ðuÞ defined by 2 oH 3 1 1 ouoHN1 ou1 6 . .. 7 .. 7 JH ðuÞ ¼ 6 4 .. . 5 . oHN 1 N 1 oH ou1 ouN1
ð26Þ
ð27Þ
is tridiagonal. The entries in the jth row of JH ðu; g; g1 Þ are given by oHj ¼ 1 16hg21 ðgj1 þ 2gj Þ 16bh2 g21 ðuj1 þ uj Þ; ouj1
ð28Þ
oHj ¼ 2 16hg21 ðgjþ1 gj1 Þ 16bh2 g21 ðuj1 þ 6uj þ ujþ1 Þ; ouj
ð29Þ
oHj ¼ 1 þ 16hg21 ð2gj þ gjþ1 Þ 16bh2 g21 ðuj þ ujþ1 Þ oujþ1
ð30Þ
and oHj =oum ¼ 0 for jj mj > 1. k;l Let gðkÞ and gk;l denote 1 denote the kth iterate corresponding to g1 , and u ðkÞ the lth iterates for u, g corresponding to g1 respectively. k;l We begin with uk;0 j ¼ j=N , and g0 0. Rewriting (22) in the form k;l k;l gjk;l ¼ gj1 þ 12hðuk;l j þ uj1 Þ;
ð31Þ
provides a procedure for obtaining recursively the components of gk;l . NewtonÕs method for the solution of (26) proceeds to obtain subsequent iterates for u as uk;lþ1 ¼ uk;l þ Duk;l ; where Duk;l satisfies the equation
ð32Þ
272
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
JH ðuk;l ÞDuk;l ¼ Hðuk;l ; gk;l ; gðkÞ 1 Þ:
ð33Þ
The iterations described by (32) and (33) may be repeated in succession until kDuk;l k1 < e1
ð34Þ
for some prescribed error tolerance e1 . Let uðkÞ denote the converged value of uk;l . Note that this is the converged result corresponding to a fixed approximate value gðkÞ 1 for g1 . Since the asymptotic condition (13) is to be imposed corresponding to g ¼ g1 , the discretization of (13) provides the additional equation needed in order to determine g1 . Let ðkÞ
wðgðkÞ 1Þ ¼
ðkÞ
ðkÞ
3uN 4uN 1 þ uN 2 : 2h
ð35Þ
We wish to determine g1 such that wðg1 Þ ¼ 0. Then, successive iterates for g1 may be obtained using the secant formula " # ðk1Þ gðkÞ ðkþ1Þ ðkÞ ðkÞ 1 g1 g1 ¼ g1 wðg1 Þ : ð36Þ ðkÞ ðk1Þ wðg1 Þ wðg1 Þ The iterations will continue until jwðgðkÞ 1 j < e2
ð37Þ
for some prescribed error tolerance e2 . We will assume that we begin our ð0Þ ð1Þ computations with initial guesses g1 and g1 for g1 . 3. Numerical results and discussion Solutions of the Falkner–Skan equation for various values of b have been reported in the literature. For b > 0, the solutions obtained are called accelerating flows; the solutions corresponding to b ¼ 0 are called constant flows; and, those corresponding to b < 0 are known as decelerating flows. Physically relevant solutions exist only for 0:19884 < b 6 2. In order to compare the results obtained by the present method with those of previously reported methods, the quantity a ¼ d2 f =dg2 for g ¼ 0 has been computed and tabulated. To express a in terms of g and n, we note that d2 f 1 du ¼ dg2 g¼0 g1 dn n¼0 and approximate it as ah ¼ ð3u0 4u1 þ u2 Þ=2h=g1 . Increased accuracy is accomplished in the computation of a by using an extrapolation formula to
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
273
Table 1 b
a (Present)
a (Previous) [5–8]
2.0000 1.0000 0.5000 0.0000 )0.1000 )0.1200 )0.1500 )0.1800 )0.1988
1.687222 1.232588 0.927680 0.469600 0.319270 0.281761 0.216361 0.128636 0.005220
1.687222 1.232589 0.927682 0.469601 0.319270 0.281762 0.216360 0.128637 0.005217
obtain aext ¼ ð4ah=2 ah Þ=3. The present numerical method has been able to obtain accelerating, constant, and decelerating flows. Table 1 compares the values of aext to the results previously reported. The converged values of u yield values of a that are in excellent agreement with those obtained in [5,6]. In general, the method behaved well for arbitrary ð0Þ ð1Þ starting values of g1 and g1 , and reasonably arbitrary choices of values for e1 and e2 . The method was able to obtain the results of the previous authors corresponding to 0:1988 6 b 6 2. We have used mesh sizes of h ¼ 0:002 and h ¼ 0:001 in order to obtain the extrapolated values reported in Table 1. Error tolerances of e1 ¼ e2 ¼ 108 have been used for all computations reported in this paper.
References [1] L. Rosenhead, Laminar Boundary Layers, Clarendon Press, Oxford, 1963. [2] H. Weyl, On the differential equations of the simplest boundary-layer problem, Ann. Math. 43 (1942) 381–407. [3] D.R. Hartree, On an equation occurring in Falkner and SkanÕs approximate treatment of the equations of the boundary layer, Proc. Camb. Philos. Soc. 33 (1937) 223–239. [4] A.M.O. Smith, Improved solutions of the Falkner and Skan boundary-layer equation, J. Aero. Sci. Sherman M. Fairchild, 1954 (Fund paper). [5] T. Cebeci, H.B. Kelle, Shooting and parallel shooting methods for solving the Falkner–Skan boundary-layer equation, J. Comp. Phys. 7 (1971) 289–300. [6] T.Y. Na, Computational Methods in Engineering Boundary Value Problems, Academic Press, New York, 1979. [7] N.S. Asaithambi, A numerical method for the solution of the Falkner–Skan equation, Appl. Math. Comput. 81 (1997) 259–264. [8] A. Asaithambi, A finite-difference method for the solution of the Falkner–Skan equation, Appl. Math. Comput. 92 (1998) 135–141.
Numerical solution of the Falkner–Skan equation using piecewise linear functions Asai Asaithambi Department of Computer Science, Parks College of Engineering and Aviation, Saint Louis University, 3450 Lindell Blvd., Saint Louis, MO 63040/63103, USA
Abstract We present a finite-element method for the solution of the Falkner–Skan equation. The method uses a coordinate transformation to map the semi-infinite domain of the problem to the unit interval ½0; 1. By means of a suitable change of variables, the transformed third-order boundary value problem is further rewritten as a system of differential equations consisting of a second-order equation and a first-order equation. The second-order equation is approximated using a Galerkin formulation with piecewise linear elements, and the first-order equation is approximated using a centereddifference approximation. In the process of the initial transformation, a finite computational boundary is introduced. This boundary is obtained as part of the computational solution by imposing an additional asymptotic boundary condition. The solutions obtained thus are in excellent agreement with those obtained by previous authors. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Falkner–Skan equation; Coordinate transformation; Galerkin formulation; Piecewise linear elements; Centered-difference; Asymptotic boundary condition
1. Introduction The Falkner–Skan equation arises in the study of laminar boundary layers exhibiting similarity. The solutions of the one-dimensional third-order
E-mail address: [email protected] (A. Asaithambi). 0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.10.047
268
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
boundary-value problem described by the well-known Falkner–Skan equation are the similarity solutions of the two-dimensional incompressible laminar boundary layer equations. This is a nonlinear two-point boundary value problem for which no closed-form solutions are available. The problem is described as follows: " 2 # d3 f d2 f df þf 2 þb 1 ¼ 0; 0 < g < 1; ð1Þ dg3 dg dg along with the boundary conditions f ¼0
at g ¼ 0;
ð2Þ
df ¼ 0 at g ¼ 0; dg
ð3Þ
df ¼ 1 as g ! 1: dg
ð4Þ
In addition to the unknown function f , the solution of (1)–(4) is characterized by the value of a ¼ d2 f =dg2 at g ¼ 0. The condition (4) is usually replaced by the condition df ¼ 1 at g ¼ g1 ; dg
ð5Þ
for some ‘‘sufficiently large’’ g1 , which must be determined as part of the solution. The addition of the new unknown g1 to the problem warrants the introduction of the ‘‘asymptotic condition’’ d2 f ¼ 0 at g ¼ g1 : dg2
ð6Þ
Rosenhead [1] and Weyl [2] present mathematical treatments of this problem focusing on existence and uniqueness results. The first numerical treatment of this problem was presented by Hartree [3]. Other numerical treatments included those developed by Smith [4], Cebeci and Keller [5], and Na [6]. These previous approaches have mainly used shooting and invariant imbedding. It is also important to note that the methods of Cebeci and Keller [5] have experienced convergence difficulties, which were overcome by moving towards more complicated methods. Recent methods presented by Asaithambi [7,8] have improved the performance of the previous methods greatly by reducing the amount of computational effort while obtaining the same results as the previous authors.
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
269
2. Method of solution With the change of variable g ¼ f =g1 , and under the transformation n ¼ g=g1 , the Falkner–Skan equation (1) is transformed to 2
g000 þ g21 gg00 þ g21 b½1 ðg0 Þ ¼ 0;
ð7Þ
where the primes denote differentiation with respect to the variable n. The boundary conditions 2, 3 and 5 are transformed to g0 ð0Þ ¼ 0;
gð0Þ ¼ 0;
g0 ð1Þ ¼ 1:
ð8Þ
The asymptotic boundary condition (6) is transformed to g00 ð1Þ ¼ 0:
ð9Þ
Direct numerical treatment of (7)–(9) is somewhat inconvenient since we have a third-order two-point boundary value problem. In order to overcome this difficulty, we let u ¼ g0 ;
ð10Þ
so that (7) transforms to u00 þ g21 gu0 þ g21 bð1 u2 Þ ¼ 0:
ð11Þ
The conditions in (8) are transformed to gð0Þ ¼ 0;
uð0Þ ¼ 0;
uð1Þ ¼ 1
ð12Þ
and the asymptotic boundary condition (9) is transformed to u0 ð1Þ ¼ 0:
ð13Þ
We wish to determine g and u for n 2 ½0; 1. We divide the interval ½0; 1 into N subintervals of length h ¼ 1=N each, and define nj ¼ jh for j ¼ 0; 1; . . . ; N . Then, we let gj and uj denote the values of gðnj Þ and uðnj Þ, respectively, for j ¼ 0; 1; . . . ; N . The boundary conditions (12) will now correspond to g0 ¼ 0;
u0 ¼ 0;
uN ¼ 1:
ð14Þ
The method of this paper is based on a finite-element formulation of (11). The method will use standard finite-difference approximations for the remaining parts of the problem. For the purpose of the finite-element formulation, we represent uðnÞ and gðnÞ in the forms uðnÞ ¼
n X j¼0
uj /j ðnÞ;
gðnÞ ¼
n X j¼0
gj /j ðnÞ;
ð15Þ
270
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
where uj uðnj Þ and gj 8 nn j1 > < nj nj1 /j ðnÞ ¼ njþ1 n > : njþ1 nj 0;
gðnj Þ as indicated previously, and for n 2 ½nj1 ; nj ; ð16Þ
for n 2 ½nj ; njþ1 ; otherwise:
Also, we represent (11) in its weak, or Galerkin form, as Z 1 ½u00 þ g21 gu0 þ g21 bð1 u2 Þ/j ðnÞ dn ¼ 0
ð17Þ
0
for j ¼ 1; 2; . . . ; N 1. For the first two integrands in (17), we use integration by parts and the fact that /j ð0Þ ¼ /j ð1Þ ¼ 0 for i ¼ 1; 2; . . . ; N 1, and obtain the equivalent form Z 1 Z 1 Z 1 u0 /0j dn g21 ug/0j dn g21 ug0 /j dn 0
þ
g21 b
Z
0
0
1 2
ð1 u Þ/j dn ¼ 0:
ð18Þ
0
Using (15) and (16) in (18), performing the integrations, collecting terms, and simplifying yields 1 uj1 2uj þ ujþ1 þ hg21 Aj þ bh2 g21 Bj ¼ 0; 6
ð19Þ
for i ¼ 1; 2; . . . ; N 1 with Aj ¼ gj1 ðuj uj1 Þ þ 2gj ðujþ1 uj1 Þ þ gjþ1 ðujþ1 uj Þ;
ð20Þ
Bj ¼ 1 121 ðu2j1 þ 2uj1 uj þ 6u2j þ 2uj ujþ1 þ u2jþ1 Þ:
ð21Þ
Next, we use the discretization gj gj1 h uj þ uj1 i ¼ 0; h 2
ð22Þ
corresponding to Eq. (10) applied at n ¼ nj1=2 for j ¼ 1; 2; . . . ; N . Note that the ð2N 1Þ equations described by (19) and (22) form a nonN N linear system of algebraic equations in the ð2N þ 3Þ variables fgj gj¼0 , fuj gj¼0 , and g1 . Among these ð2N þ 3Þ variables, the three quantities g0 ; u0 , and uN are available from (14). Assuming that g1 is available, (19) and (22) will be used to N N 1 determine the ð2N 1Þ unknown quantities fgj gj¼1 , and fuj gj¼1 . Then, a numerical representation of the asymptotic boundary condition (13) will be used to obtain an improved value for g1 . We now proceed to describe the iterative process for the solution of the nonlinear system (19) and (22). We begin by letting
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
uT ¼ ½ u1
uN ;
gT ¼ ½ g1
gN
271
ð23Þ
and 2
3 H1 ðu; g; g1 Þ 6 7 .. Hðu; g; g1 Þ ¼ 4 5; . HN 1 ðu; g; g1 Þ
ð24Þ
where Hj ðu; g; g1 Þ ¼ ujþ1 2uj þ uj1 þ 16hg21 Aj þ bh2 g21 Bj
ð25Þ
for j ¼ 1; 2; . . . ; N 1. If we consider (19) as a nonlinear system in u alone, then solving (19) is equivalent to solving the system described by Hðu; g; g1 Þ ¼ 0: Note that the Jacobian matrix JH ðuÞ defined by 2 oH 3 1 1 ouoHN1 ou1 6 . .. 7 .. 7 JH ðuÞ ¼ 6 4 .. . 5 . oHN 1 N 1 oH ou1 ouN1
ð26Þ
ð27Þ
is tridiagonal. The entries in the jth row of JH ðu; g; g1 Þ are given by oHj ¼ 1 16hg21 ðgj1 þ 2gj Þ 16bh2 g21 ðuj1 þ uj Þ; ouj1
ð28Þ
oHj ¼ 2 16hg21 ðgjþ1 gj1 Þ 16bh2 g21 ðuj1 þ 6uj þ ujþ1 Þ; ouj
ð29Þ
oHj ¼ 1 þ 16hg21 ð2gj þ gjþ1 Þ 16bh2 g21 ðuj þ ujþ1 Þ oujþ1
ð30Þ
and oHj =oum ¼ 0 for jj mj > 1. k;l Let gðkÞ and gk;l denote 1 denote the kth iterate corresponding to g1 , and u ðkÞ the lth iterates for u, g corresponding to g1 respectively. k;l We begin with uk;0 j ¼ j=N , and g0 0. Rewriting (22) in the form k;l k;l gjk;l ¼ gj1 þ 12hðuk;l j þ uj1 Þ;
ð31Þ
provides a procedure for obtaining recursively the components of gk;l . NewtonÕs method for the solution of (26) proceeds to obtain subsequent iterates for u as uk;lþ1 ¼ uk;l þ Duk;l ; where Duk;l satisfies the equation
ð32Þ
272
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
JH ðuk;l ÞDuk;l ¼ Hðuk;l ; gk;l ; gðkÞ 1 Þ:
ð33Þ
The iterations described by (32) and (33) may be repeated in succession until kDuk;l k1 < e1
ð34Þ
for some prescribed error tolerance e1 . Let uðkÞ denote the converged value of uk;l . Note that this is the converged result corresponding to a fixed approximate value gðkÞ 1 for g1 . Since the asymptotic condition (13) is to be imposed corresponding to g ¼ g1 , the discretization of (13) provides the additional equation needed in order to determine g1 . Let ðkÞ
wðgðkÞ 1Þ ¼
ðkÞ
ðkÞ
3uN 4uN 1 þ uN 2 : 2h
ð35Þ
We wish to determine g1 such that wðg1 Þ ¼ 0. Then, successive iterates for g1 may be obtained using the secant formula " # ðk1Þ gðkÞ ðkþ1Þ ðkÞ ðkÞ 1 g1 g1 ¼ g1 wðg1 Þ : ð36Þ ðkÞ ðk1Þ wðg1 Þ wðg1 Þ The iterations will continue until jwðgðkÞ 1 j < e2
ð37Þ
for some prescribed error tolerance e2 . We will assume that we begin our ð0Þ ð1Þ computations with initial guesses g1 and g1 for g1 . 3. Numerical results and discussion Solutions of the Falkner–Skan equation for various values of b have been reported in the literature. For b > 0, the solutions obtained are called accelerating flows; the solutions corresponding to b ¼ 0 are called constant flows; and, those corresponding to b < 0 are known as decelerating flows. Physically relevant solutions exist only for 0:19884 < b 6 2. In order to compare the results obtained by the present method with those of previously reported methods, the quantity a ¼ d2 f =dg2 for g ¼ 0 has been computed and tabulated. To express a in terms of g and n, we note that d2 f 1 du ¼ dg2 g¼0 g1 dn n¼0 and approximate it as ah ¼ ð3u0 4u1 þ u2 Þ=2h=g1 . Increased accuracy is accomplished in the computation of a by using an extrapolation formula to
A. Asaithambi / Appl. Math. Comput. 159 (2004) 267–273
273
Table 1 b
a (Present)
a (Previous) [5–8]
2.0000 1.0000 0.5000 0.0000 )0.1000 )0.1200 )0.1500 )0.1800 )0.1988
1.687222 1.232588 0.927680 0.469600 0.319270 0.281761 0.216361 0.128636 0.005220
1.687222 1.232589 0.927682 0.469601 0.319270 0.281762 0.216360 0.128637 0.005217
obtain aext ¼ ð4ah=2 ah Þ=3. The present numerical method has been able to obtain accelerating, constant, and decelerating flows. Table 1 compares the values of aext to the results previously reported. The converged values of u yield values of a that are in excellent agreement with those obtained in [5,6]. In general, the method behaved well for arbitrary ð0Þ ð1Þ starting values of g1 and g1 , and reasonably arbitrary choices of values for e1 and e2 . The method was able to obtain the results of the previous authors corresponding to 0:1988 6 b 6 2. We have used mesh sizes of h ¼ 0:002 and h ¼ 0:001 in order to obtain the extrapolated values reported in Table 1. Error tolerances of e1 ¼ e2 ¼ 108 have been used for all computations reported in this paper.
References [1] L. Rosenhead, Laminar Boundary Layers, Clarendon Press, Oxford, 1963. [2] H. Weyl, On the differential equations of the simplest boundary-layer problem, Ann. Math. 43 (1942) 381–407. [3] D.R. Hartree, On an equation occurring in Falkner and SkanÕs approximate treatment of the equations of the boundary layer, Proc. Camb. Philos. Soc. 33 (1937) 223–239. [4] A.M.O. Smith, Improved solutions of the Falkner and Skan boundary-layer equation, J. Aero. Sci. Sherman M. Fairchild, 1954 (Fund paper). [5] T. Cebeci, H.B. Kelle, Shooting and parallel shooting methods for solving the Falkner–Skan boundary-layer equation, J. Comp. Phys. 7 (1971) 289–300. [6] T.Y. Na, Computational Methods in Engineering Boundary Value Problems, Academic Press, New York, 1979. [7] N.S. Asaithambi, A numerical method for the solution of the Falkner–Skan equation, Appl. Math. Comput. 81 (1997) 259–264. [8] A. Asaithambi, A finite-difference method for the solution of the Falkner–Skan equation, Appl. Math. Comput. 92 (1998) 135–141.