- Email: [email protected]
Simple analytic approximations for the Blasius problem
Physica D 310 (2015) 72–78
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
Simple analytic approximations for the Blasius problem R. Iacono a,∗ , John P. Boyd b a
ENEA, CR Casaccia (SP 118), Via Anguillarese, 301, 00123 Roma, Italy
b
Department of Atmospheric, Oceanic & SpaceScience, University of Michigan, 2455 Hayward Avenue, Ann Arbor MI 48109, United States
highlights • • • •
We derive analytic approximations for the Blasius function f and its derivatives. We extend the integral iteration scheme for the Blasius problem devised by H. Weyl. We compute very accurate bounds for the second derivative of f at the origin. We discuss the new approximations in the context of generalized Padè theory.
article
info
Article history: Received 16 April 2015 Received in revised form 6 August 2015 Accepted 6 August 2015 Available online 13 August 2015 Communicated by M. Vergassola
abstract The classical boundary layer problem formulated by Heinrich Blasius more than a century ago is revisited, with the purpose of deriving simple and accurate analytical approximations to its solution. This is achieved through the combined use of a generalized Padé approach and of an integral iteration scheme devised by Hermann Weyl. The iteration scheme is also used to derive very accurate bounds for the value of the second derivative of the Blasius function at the origin, which plays a crucial role in this problem. © 2015 Elsevier B.V. All rights reserved.
Keywords: Blasius equation Boundary layer Generalized Padé approximants Integral iteration scheme
1. Introduction The Blasius problem is a prototypal example illustrating the development of a boundary layer as a result of the interaction between a weakly viscous flow and an obstacle, a phenomenon that is ubiquitous in nature. By now, it has more than a century of life, during which it has developed a vast bibliography; detailed information about this problem, together with a reminder of the essential literature, can be found in the recent review [1]. Here we only sketch the basic physical picture. Consider a steady, uniform flow impinging on a thin, wide plate parallel to it, and choose a (X , Y ) Cartesian coordinate frame such that the flow points in the positive X direction and the plate is placed along the half-plane Y = 0, X > 0. For X > 0, the stream is basically undisturbed at large Y , but the flow must slow down when approaching the plate, because the fluid velocity must be zero at
∗
Corresponding author. E-mail addresses: [email protected] (R. Iacono), [email protected] (J.P. Boyd). http://dx.doi.org/10.1016/j.physd.2015.08.003 0167-2789/© 2015 Elsevier B.V. All rights reserved.
the solid boundary. In a low viscosity fluid, such as air or water, this only happens in a narrow region near the plate, named ‘‘boundary layer’’, or ‘‘shear layer’’, in which viscosity effects become essential. The main problem then is to compute the flow in this layer, and to match it with the flow outside. Scale analysis of the steady-state Navier–Stokes equations shows that, in a first approximation, pressure gradients can be neglected within the boundary layer, so that the dominant balance, in the longitudinal (X ) direction, is between the nonlinear advection terms and the leading viscosity term,
∂u ∂u ∂ 2u +v =ν , ∂X ∂Y ∂Y 2 where (u, v) are the (X , Y ) components of the flow, and ν is the kinematic viscosity. Together with mass continuity, ∂X u +∂Y v = 0,
u
this equation can be solved to yield the steady flow components in the boundary layer. Alternatively, the equilibrium equations can be reduced to a single third-order, nonlinear partial differential equation (PDE) for the flow streamfunction (see, e.g., [1]). It turns out, however, that the problem can be further simplified, thanks to the existence of a continuous group of invariance,
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
Fig. 1. The Blasius function f (x) [solid black] is plotted with its first derivative [red dashed] and second derivative [blue dots]. The function (x/f (x))df /dx, which has a starring role in later derivations, is the thick green curve with markers.
which also has a physical meaning, being proportional to the shear stress on the plate. The parameter is a priori unknown, and not easy to compute analytically (we will come back to this issue later on). However, a second group invariance of the BE, discovered by Töpfer [3], allows for an efficient numerical computation of κ , which is approximately given by κ ≈ 0.332057336. High precision values of κ and of the other constants that enter the Blasius problem can be found in Table 4.1 of [1]. Although the existence of a complex pole had long been known [4], a detailed understanding of the structure of the singularities of the BE in the complex plane has been achieved much more recently [2]. Among other results, it was found that the three singularities closest to the origin, which determine the radius of the convergence of the series expansion, lie on a circle of radius S, where S ≈ 5.69. It was also shown in [2] that the series can be accelerated using Euler’s summation, and that the accelerated series appears to converge everywhere on the positive x axis, which is the physical domain of the flow. Treating the BE as a linear, first order ODE for fxx , yields the wellknown expression
fxx = κ exp − pointed out by Blasius himself. This implies that the streamfunction is not a function of X and Y separately, but is instead a function of a similarity variable x = Y /X 1/2 . As a consequence, the problem becomes one-dimensional, with the third-order PDE reducing to a third-order, nonlinear ordinary differential equation (ODE), 2fxxx + ffxx = 0,
(1)
which is referred to as the Blasius equation (BE hereafter). Here, f is related to the streamfunction of the two-dimensional problem, and the subscript denotes derivation with respect to the similarity variable x. Eq. (1) needs to be solved subject to boundary conditions both at the origin and at infinity (see [1] for further details): f (0) = fx (0) = 0,
fx (∞) = 1.
(2)
The boundary problem (1)–(2) was solved in [2] using the Chebyshev pseudospectral method. Fig. 1 shows the resulting profiles of f (x) and of its first two derivatives, together with that of the function xfx /f , which will be useful in the developments to follow. This very accurate numerical solution – errors for f and its first three derivatives are smaller than 10−9 over the whole real axis – will be used as a reference in the present work. All the functions plotted in Fig. 1 are smooth and monotonic. The Blasius function f has a simple structure: it grows like x2 at small x, and asymptotes to x plus a negative constant at large x. Despite the simple shape of the solution and the long history of the subject, analytic insight on the BE is still limited. While an accurate approximation for fxx is known (see the next section), to the best of our knowledge, nobody has yet been able to provide a simple analytic approximation for f (x) of comparable accuracy. Constructing such an approximation is the main objective of the present note. 2. Some basic results The series expansion of f about the origin, f (x) ≈
κ 2
x2 −
κ2 240
x5 +
11κ 3 161280
x8 −
5κ 4 4257792
x11 + · · · ,
(3)
was computed by Blasius himself, who derived a recurrence relation for the coefficients of the expansion. The coefficients depend on the parameter
κ ≡ fxx (0),
(4)
73
1 2
x
dξ f (ξ ) ,
(5)
0
which explains the fast decay of fxx with x shown by Fig. 1. Expression (5) also suggests that a guess for f that is good at small x should give an approximation for fxx accurate over a wider x-range, because the main contribution in the rhs of (5) comes from small x. One could place the series expansion (3) in (5), as Bairstow [5] did, but would then be confronted with the fact that the series has a finite radius of convergence. Bairstow noted, however, that taking f ≈ (κ/2)x2 yields
fxxCG ≡ fxx = κ exp −
κ 12
x3 ,
(6)
where CG stands for cubic-Gaussian. The approximation (6) is fairly accurate for x not too large: it is still within 1% from the numerical solution at x = 2.5, whereas the guess for f only preserves a similar accuracy up to x = 1.5. It turns out that Bairstow’s approach can be refined, to build a much better approximation for fxx , fxxPBS ≡ fxx = κ
exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3
,
(7)
which will be used here to test the accuracy of the approximations to be derived. The – non-trivial – derivation of this reference approximation, due to Parlange, Braddock and Sander [6], is sketched in the review [1]. x Finally, we note that integration by parts of the identity f (x) = dξ fξ yields an integral expression for f in terms of its second 0 derivative, f (x) =
x
dξ (x − ξ ) fξ ξ (ξ ).
(8)
0
This suggests an iterative approach based on (5) and (8): given a guess for f that satisfies the boundary conditions at the origin, one could compute fxx through (5), then get a new guess for f using (8), and so on. Would this approach prove convergent, it could be used to derive accurate numerical approximations to the solution of the BE. 3. Weyl’s integral approach The idea of iterating between f and fxx is not new. It goes back to a trio of articles on the BE and related boundary layer problems that Hermann Weyl published three quarters of a century ago.
74
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78 Table 1 Maximum errors in the iterates of Weyl’s integral for the second derivative of the Blasius function, W (z ) = fxx (z (4/κ)1/3 )/κ 1/3 . Chebyshev series with N = 200 on z ∈ [0, 8], but there is no significant change when N = 80 and the Chebyshev approximation interval is reduced to z ∈ [0, 6].
His second contribution [7] is a three-page note that completes the convergence proof for his iteration for the generalized Blasius equation, 2fxxx + ffxx + (k2 − fx2 ) = 0 where k is a parameter (a constant), but says nothing about the BE itself. His final paper [8] amplifies his earlier treatment of the Blasius problem as a prelude to expanding his methodology to similar nonlinear ODEs, but offers nothing new about the solutions to the BE. We therefore shall discuss only his first article [9] where Weyl repeats Töpfer’s group invariance and then continues ‘‘When H. Blasius carried out the integration . . . he made use of the power series for f around z = 0 and of a certain asymptotic expression for large values of z, adjusting the constant β so as to make both expressions dovetail in a middle region. It is a fact, although it seems to have escaped notice, that convergence of the power series of f stops somewhere between √ [converting from√ his convention to the usual definition of f ] 3 9/κ ≈ 3.004 and 3 30/κ ≈ 4.487. Thus, Blasius’ method . . . can yield results of limited accuracy only’’. Weyl was incorrect on a couple of counts. The radius of convergence is S ≈ 5.69, well outside the limits that he gives (see Section 8.2 for more details), and the Euler acceleration of the power series converges on the entire real axis. Nevertheless, his negative evaluation of prior work motivated him to set up an iterative scheme, based on an integral equation for the second derivative of f , which he proved to be convergent. Despite Weyl’s great reputation, his iteration has not been widely used. We extend his treatment in three ways. First, we use a modern computer algebra system (Maple) to extend the analytic formulas for iterates to one order higher than Weyl himself. Second, we show that with modern Chebyshev spectral methods, it is possible to numerically apply his method for arbitrary iteration order at a cost of only O(N log(N )) operations per iteration where N, the number of Chebyshev polynomials, is no larger than O(100) and where the error in the seventh iterate is a maximum of only about 10−14 . Third, we show that the numerical iterates can be used to compute very accurate bounds for the value of κ . Instead of explicitly iterating between (5) and (8), Weyl works with the nonlinear integral equation z
W (z ) = exp −
(z − y)2 W (y)dy ,
(9)
0
where W (z ) = fxx (z (4/κ)1/3 )/κ 1/3 and z ≡ (κ/4)1/3 x, and proves that this equation yields a convergent iteration scheme. The first few iterates are, with the maximum error of the iterate shown in square brackets, W0 = 0 [1.0]
(10)
W1 = 1 [1.0]
(11)
W2 = exp(−[1/3]z ) [0.23]
(12)
z 2 3 W3 = exp − (z − y) exp(−[1/3]y )
(13)
0
= exp(−w(z )) [0.00039]
(14)
where 3
3 3
2 2 1 3 Γ −Γ , z 3
∞
z
exp(−t ) t α−1 dt .
1 0.0235 0.000391 2.89 × 10−6 1.13 × 10−8 2.61 × 10−11 3.84 × 10−14
Taking into account the normalizations, it is easily seen that the second iterate, Eq. (12), is none other than Bairstow’s approximation (6). Although explicit forms for the iterates are not possible for order k > 3, a very fast numerical procedure is possible. First, the exponential of each iterate can be written without approximation as p(z )
z
(z − y)2 Wk−1 (y)dy
=
(17)
0
= z2
z
Wk−1 (y)dy
(18)
0 z
yWk−1 (y)dy +
− 2z 0
z
y2 Wk−1 (y)dy
(19)
0
Wk (z ) = exp {−p(z )} .
(20)
We can approximate each integrand on an arbitrary interval z ∈ [0, b] by a Chebyshev interpolant with N interpolation points. Because of the extremely rapid exponential decay of W (z ) with z, b = 6 gives an approximation accurate to near machine precision (2.2 × 10−16 ) in Matlab. We can then apply the usual recurrence to obtain the coefficients of the indefinite integral from those of the integrand, evaluate the integrand series at each interpolation point zj , and add the contributions of the three integrals to obtain Wk (zj ) = exp(−p(zj )). By using the Fast Fourier Transform for the interpolation and summation, each iterate can be calculated in O(N log2 (N )) operations. Convergence is exponential as shown in Table 1. b Elementary calculus requires that, for arbitrary f (x), a fxx dx = fx (b) − fx (a). The boundary conditions for the derivative of the Blasius function, fx (0) = 0 and fx (∞) = 1, then give the valuable identity ∞
fxx dx = 1.
(21)
3 3
This can be used to approximate κ , the second derivative of the Blasius function at the origin. First, though, we must sort the different normalization used by Weyl.1 Let f denote the Blasius function with Weyl’s scaling so that fzz = W (z ). Then the usual Blasius function f (x) with fxx (0) = κ is connected to f by (22)
It follows that (15)
fxx (x) = κ f′′ ([κ/4]1/3 x)
= κ W ([κ/4]1/3 x).
with Γ (α, z ) the usual incomplete gamma function,
Γ (α, z ) ≡
1 2 3 4 5 6 7
f (x) = 42/3 κ 1/3 f([κ/4]1/3 x).
1 1 3 1 3 −Γ , z w(z ) = 1 − e−1/3 z + z 2 3−2/3 Γ − 2z3
maxz ∈[0,∞] |Wk − W |
0
3
−1/3
Iteration number k
(16)
1 f + 2ff = 0, f (0) = 1 versus 2f + ff = 0, f (0) = κ . zzz zz zz xxx xx xx
(23)
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
75
Table 2 Approximations to κ and relative errors as obtained from the integral of Weyl’s iterates. No answer is given for W1 ≡ 1 because the infinite integral of this iterate is unbounded. Iteration number k
κk
|κ − κk |/κ
2 3 4 5 6 7 8
0.342095321717019 0.331898302513614 0.332058458263831 0.332057331979152 0.332057336224705 0.332057336215182 0.332057336215196
0.03023 0.00047893 3.3791e−06 1.2757e−08 2.8634e−11 4.1961e−14 1.6717e−16
The integral (21) then becomes ∞
κ W ([κ/4]1/3 x)dx = 1,
(24)
0
and, with the change of variable ζ ≡ (κ/4)1/3 x, ∞
W (ζ )dζ = 4−1/3 κ −2/3 .
(25)
0
This allows us to calculate κ from the integral of Weyl’s function as
κk =
1 2
∞ 0
Wk (ζ ) dζ
−3/2
(26)
Table 2 shows that the resulting approximations converge very rapidly. The eighth iterate is correct to all fifteen decimal places shown in the table. It should be noted that, since Weyl’s iteration produces alternating bounds for W (z ), (26) yields alternating bounds for κ . Only the first upper bound (iteration number 2) was computed in [9]. 4. New analytic approximations for f and fxx An approach often used to extend the accuracy of the series expansions beyond their convergence range is that of the Padé approximants, which are rational approximations constructed from the series (see, e.g., [10]). In the context of the Blasius problem, this approach was pursued in [11], in the attempt of constructing approximations for f that would also allow the determination of the parameter κ . Because of the structure of (3), a Padé approximant for f should be constructed as x2 times the ratio of two polynomials in x3 , but such an approximant cannot capture the large-x behavior, which is linear in x. High order, diagonal Padé approximants of this form were computed in [11] that proved very accurate at small and intermediate x, but no simple way was found to extend this accuracy to larger values of x. Here we tackle the problem from a different angle. Since it is difficult to construct good rational approximations for f , we may try to do so for some other variable that is best suited. A good candidate seems to be the variable y ≡ xfx /f that we plotted in Fig. 1. This is a smooth, decreasing function that equals 2 at the origin, tends towards 1 at large x, and is a function of x3 at small x. The simplest rational approximation for y is therefore of the form y ≈ (2 + α x3 )/(1 + α x3 ), with α a constant that can be determined by enforcing the x3 term in the series expansion of y about the origin. This yields y=
2 + (κ/40)x3 1 + (κ/40)x3
,
(27)
which is very accurate at small x, and is still within 1% from the numerical y(x) at x = 4.8. Upon integrating, we get a simple analytic approximation for f : fA =
κ x2 /2 . [1 + (κ/40)x3 ]1/3
(28)
Fig. 2. Absolute errors in approximations to fxx , the second derivative of the Blasius function. The thin solid red curve is the Parlange–Braddock–Sander approximation, fxxPBS ≡ κ exp(−0.0125152x3 )/(1 + 0.00505209x3 )3 . The dashed κ 3 2/3 black curve is fxxa ≡ κ exp −5 (1 + 40 x ) − 1 . The thick black curve is fxxb ≡
m A exp − 4p (1 + px3 )2/3 −
n 2q
(1 + qx3 )1/3 , with the coefficients given in Eqs. (43)
and (45).
This is already a pretty accurate approximation for small and intermediate values of x, with relative errors of 2.6 × 10−7 at x = 1, 4.4 × 10−4 at x = 4, and less than 3% at x = 8. Then it loses accuracy, because it does not have the correct asymptotic dependence, even though it captures the linear growth with x. Based on the consideration of the previous sections, we would expect that placing (28) in (5) should give a good approximation for fxx . This is indeed the case. Integration of (5) gives
fxx ≡ fxxa ≈ κ exp −5 (1 +
κ
40
x3 )2/3 − 1
,
(29)
which yields the absolute errors plotted as a dashed line in Fig. 2. The errors for (29) are smaller than those for the reference approximation (7) (red line) up to x ≃ 6, and the maximum error, 3.7816 × 10−5 , is smaller than the maximum error for (7), which is of 1.0176 × 10−4 . Thus, quite surprisingly, we have found a simple analytic approximation for fxx that is superior to the sophisticated approximation by Parlange et al. [6]. Placing (29) in (8) gives a new approximation for f : fB = κ
x
dξ (x − ξ ) exp −5
1+
0
κ 40
ξ3
2/3
−1 .
(30)
This time, the integral cannot be evaluated analytically, but f can be easily computed numerically. As shown in Fig. 3, this yields a very accurate approximation (dashed black curve), with relative errors smaller than 10−5 up to x ≃ 6, and still smaller than 10−4 up to x ≃ 100. Errors are much smaller than those for the approximation (28) (solid black curve), confirming that the iterative scheme based on (5) and (8) indeed works well. 5. Improving the small x behavior The form of (28) suggests a simple generalization. Define a function g so as to strip off the leading quadratic factor of f , enabling us to work with a function of z ≡ x3 only: f =
κ 2
x2 g (z )
↔ g (z ) =
2
κ x2
f (x).
(31)
76
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
Fig. 3. Relative errors in approximations to f , the Blasius function. Solid black: f A ≡ (κ/2)x2 /(1 + (κ/40)x3 )1/3 . Triangles are for the asymptotic approximation f asy,best , defined in Eq. (49).
Then, (28) can be considered as the lowest order [N /(N + 1)] Padé approximant to the cube of g, g 3 , i.e.,
g[0/1] =
1/3
1 1 + (κ/40)z
.
(32)
Higher order approximants are easily computed: 11200(396 + 43 κ z )
g[1/2] =
1/3
4435200 + 592480 κ z + 12073 κ 2 z 2
g[2/3] =
1663200
P2
1/3
Q3
,
(33)
,
(34)
with
Fig. 4. Errors in approximations to g (z ) which are each the cube root of a Padé approximation to g 3 (z ). The error in the first two orders of the large z-approximation is also shown (dashed line). The Padé orders label each curve: degree of the numerator polynomial in z on the left, denominator degree to the right of the slash. The horizontal dotted guideline is at 10−5 . The vertical guideline is at the radius of convergence of the power series.
illustrates the accuracy gain at small and intermediate values of x obtained by increasing the order of the approximant. 6. Improving the large x behavior The approximations derived in the previous sections are quite good, but could still be improved at large x, were we able to enforce the correct asymptotic behavior. Consider, for example, the lowest order approximation fA [Eq. (28)]. Since f ∼ B + x at large x, with B ≈ −1.7207876575, we can try to improve the large x behavior by adding a term with the same structure of the one in fA , but with a 2/3 power in the denominator, which would tend towards a constant in the limit x → ∞. That is, we may seek an approximation of the form fC =
P2 = 3964295171200 + 76608724296 κ z
+ 507276173 κ 2 z 2 ,
(35)
(1 + px3 )1/3
+ 4078150639904880 κ z + 16713554600203 κ z , (36) 3 3
g[3/4] =
9828 P3 5
Q4
,
1/3 (37)
nx2
(1 + qx3 )2/3
.
(40)
n/q2/3 = B,
(41)
which gives the correct large x behavior, and m + n = κ/2,
and
+
We determine the coefficients in (40) by asking that m/p1/3 = 1,
Q3 = 6593415728739840000 + 292251023467603200 κ z 2 2
mx2
mp + 2nq = κ 2 /80,
which enforces the first two terms of the series expansion about the origin. Solving the system (41)–(42) with Mathematica, we find a single real solution, given by
with
m = 0.198177262,
n = −0.0321485939,
P3 = 52740459296753088904184064000 + 1156660146147189506418534400 κ z
p = 0.0077832588,
q = 0.0025535952.
+ 9352330559888045666104080 κ 2 z 2 + 11843596516363041045577 κ 3 z 3 ,
(38)
Q4 = 103666646793697871550064196198400 + 4865197353109362482567876121600 κ z
+ 75992549747403848831924652000 κ 2 z 2 + 416014006359213902061551220 κ 3 z 3 + 401763679454602196561927 κ 4 z 4 .
Recalling that κ = 0.33205733621519630, we may compute the relative errors (absolute errors divided by g) for these approximations, which are displayed in Fig. 4. The figure clearly
(43)
The relative errors for this approximation are shown as red dots in Fig. 3. The approximation is slightly less accurate than fA [Eq. (28)] at small x, but much more accurate at large x, with the largest relative errors, at intermediate values of x, that remain smaller than 1%. Placing (40) in (5) gives a new analytic approximation for the second derivative of f ,
(39)
(42)
fxxb = A exp −
m 4p
(1 + px3 )2/3 −
n 2q
(1 + qx3 )1/3 ,
(44)
where A = κ exp
m 4p
+
n 2q
≈ 0.356393092.
(45)
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
The absolute errors in this approximation are graphed as a thick black curve in Fig. 2. The figure shows that (44) is definitely more accurate than the reference Parlange–Braddock–Sander approximation, with a maximum error that is about 1/8 of that for fxxPBS . Finally, we can place (44) in (8) to get a new approximation for f
8.2. Recurrence relation for power series In (4.2) of [1], the right-hand side should be divided by a factor of two. The correct power series recurrence is am = −
d ξ (x − ξ )
fD = A
×
0
2 (3m − 1)(3m − 2)(3m − 3)
(3j − 1) (3j − 2) aj am−j
(51)
j =1
m n × exp − (1 + pξ 3 )2/3 − (1 + qξ 3 )1/3 . 4p
2q
(46)
This yields the relative error curve marked by stars in Fig. 3, with a maximum error just above 10−5 at x = 100. 7. Improved large x asymptotics
fxx (x) ∼ Q exp(−(1/4)(x + B) ). 2
2 (47)
This is derived by substituting f ∼ B + x into 2fxxx + ffxx = 0 and solving as a first order ODE for fxx with the initial condition fxx (0) = κ . It turns out that Q = 0.2337276212850634. Integrating this twice, both with lower limits at zero, gives
− 2 exp(−(1/4)B2 ) + 2 exp(−(1/4)(x + B)2 ) .
1 120
f int ∼ −1.48925285928230415948146 (48)
We can therefore get a better asymptotic approximation by writing f asy,best ∼ (B + 1.48925285928230415948)
2
1
n (52)
36
exp(−(1 − 3β1 )y)
(53)
(1 + β1 y)3 exp(−0.452277442494834y)
(54)
(1 + 0.182574185835055y)3 exp(−(1 − 3β1 )(κ/12)x3 )
(55)
(1 + [β1 κ/12]x3 )3 exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3
.
(56)
Eq. (5.12) of [1] is in error, writing 0.4114 in the exponential instead of the correct 0.45227. However, (5.17) in the same article is correct.
− 2 exp(−(1/4)B2 ) (49)
The errors of (49) are graphed as triangles in Fig. 3. The approximation is superb for large x but poor for small x, and thus complementary to most of the other approximations derived here. 8. Corrections We are aware of several previously published typographical errors. In the spirit that repair is as important as new construction, we correct these errors in this section even though these errors are not directly related to any of our new approximations. 8.1. Large x approximation typo
9. Conclusion Padé approximants have proved so useful that they have been greatly generalized. Although we have not invoked generalized Padé theory in our derivations, we would be amiss in failing to note that many of our approximations fit in this generalized framework. A standard Padé approximant is the solution to the linear polynomial equation QN (z )fM .N − PM (z ) = 0 where the coefficient polynomials PM and QN are chosen so that the power series expansion of f[M .N ] agrees with the leading terms of f (z ). In his thesis, Padé – using ideas developed two decades earlier by his thesis advisor, Hermite – observed that there was no reason to restrict the polynomial equation to first degree. For example, one could define an approximation to f (x) by declaring that g[N ,K ,L,M ] is the solution to the cubic polynomial equation
Q(z )g[3N ,K ,L,M ] (z ) + R(z )g[2N ,K ,L,M ] + Sg[N ,K ,L,M ] − P(z ) = 0
The lowest order asymptotic approximation for large x is fxx ∼ exp(−B2 /4)Q exp(−(1/4)x(x + 2B)).
fxx = κ fxx = κ
+ (1 − .735875425275434309752)x √ π (B + x) [erf(x/2 + B/2) − erf(B/2)] +Q
x ≫ 1.
≤ cn ≤
1
8.3. PBS approximation
=κ
+ 0.735875425275434309752792 x.
n
where both the left and right hand sides are smaller by 4−n than the incorrect bounds in Weyl [9].
fxx = κ
At large x, this gives
Parlange, Braddock and Sanders give the approximation, using y = (κ/12)x3 ,
√ π (B + x) [erf(x/2 + B/2) − erf(B/2)]
+ 2 exp(−(1/4)(x + B)2 )
where the factor erroneously omitted in [1] is boxed. Weyl erroneously placed this same factor of two in the numerator instead of the denominator in [9], pg. 579. His power series coefficients cn [the coefficient of z 3n+2 ] are therefore too large by 4n . His error shrinks the radius of convergence (noting that the series proceeds in powers of z 3 ) by 41/3 (to 3.58 in our notation). His bounds on his coefficients cn should be 1
The asymptotic behavior of fxx is given by
f int ∼ Q
1
m−1
x
77
(50)
The exponential exp(−B2 /4) was incorrectly omitted from (3.2) of [2] and the error repeated in [1].
(57)
and where the coefficients in the cubic equation are polynomials in z only of degrees (N , K , L, M ) respectively. Standard Padé approximants to g (= f /(κ x2 /2)) would not be able to capture the decay of g as z 1/3 at large z. We therefore computed Padé approximants to the cube root of g. This is
78
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
equivalent to a cubic Hermite–Padé approximation for the special case K = L = 0 so that only the constant and the cubic term are nonzero. QN (z )g[3N ,0,0,M ] (z )] − PM (z ) = 0
↔ g[3N ,0,0,M ] (z )] =
PM QN
. (58)
Unfortunately, solutions to polynomial equations have only algebraic branch point singularities. More complex behavior can be obtained by replacing the polynomial equation by a simple differential equation. For example, one species of ‘‘Padédifferential approximant’’ can be defined as the solution of the first-order linear differential equation
A = κ exp
m
+
4p
n
≈ 0.356393092.
2q
(A.5)
f asy,best ∼ (B + 1.48925285928230415948) + (1√ − .735875425275434309752)x +Q π (B + x) [erf(x/2 + B/2) − erf(B/2)]
− 2 exp(−(1/4)B2 ) + 2 exp(−(1/4)(x + B)2 ) ,
x ≫ 1.
(A.6)
A.2. Approximations to the second derivative
x q(y) dy u = exp −
ux + q(x)u = 0
where
(59)
0
where q(x) is chosen to be a power series whose coefficients are such that the leading terms of the power series for u(x) match those of the function which is to be approximated by u(x). This is exactly what we have done with u ≡ fxx and q(x) ≡ 2f (x). Acknowledgments
fxxCG = κ exp − fxxPBS = κ
κ 12
x3 .
exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3 κ 3 2/3 fxxa = κ exp −5 (1 + x ) −1 . fxxb
Appendix. Summary of approximations
with A given by (A.5).
A.1. Approximations to the Blasius function
References
fA =
κ x /2 . [1 + (κ/40)x3 ]1/3
(A.1)
[fA is equivalent to (κ/2)x2 g[0/1] ] fB = κ
x
dξ (x − ξ ) exp −5 0
mx2
1+
κ 40
nx2
ξ3
2/3
−1 .
, (1 + px3 )1/3 (1 + qx3 )2/3 m = 0.198177262, n = −0.0321485939, p = 0.0077832588, q = 0.0025535952. x fD = A dξ (x − ξ ) 0 m n × exp − (1 + pξ 3 )2/3 − (1 + qξ 3 )1/3 fC =
+
4p
2q
(A.2)
(A.3)
(A.4)
.
40
We thank Tyler Nielson and his project group at Eastern Washington University for finding the error in the power series recurrence (4.2) of [1].
2
(A.7)
n m 3 2/3 3 1/3 = A exp − (1 + px ) − (1 + qx ) 4p
2q
(A.8) (A.9)
(A.10)
[1] J.P. Boyd, Themes illustrated by the blasius function: Numerical computations before computers, the value of narrow tricks, and interesting undergraduate projects and open research problems, SIAM Rev. 50 (2008) 791–804. [2] J.P. Boyd, The Blasius function in the complex plane, J. Exp. Math. 8 (1999) 381–394. [3] K. Töpfer, Bemerkung zu dem Aufsatz von H. Blasius Grenzschichten in Flüssigkeiten mit kleiner Reibung, Zeit. Math. Phys. 60 (1912) 397–398. [4] B. Punnis, Zur Differntialgleichung der Plattengrenzschicht von Blasius, Arch. d. Math. 7 (1956) 165–171. [5] L. Bairstow, J. R. Aero. Soc. 19 (1925) 3. [6] J. Parlange, R.D. Braddock, G. Sander, Analytical approximations to the solution of the Blasius equation, Acta Mech. 38 (1981) 119–125. [7] H. Weyl, Concerning the differential equations of some boundary layer problems, Proc. Natl. Acad. Sci. 28 (1942) 100–101. [8] H. Weyl, On the differential equations of the simplest boundary layer problems, Ann. of Math. 43 (1942) 381–407. [9] H. Weyl, Concerning the differential equations of some boundary layer problems, Proc. Natl. Acad. Sci. 27 (1941) 578–583. [10] G.A. Baker Jr., Padé approximants, in: Advances in Theoretical Physics, in: Advances in Theoretical Physics, vol. 1, Academic Press, New York, 1965, pp. 1–50. [11] J.P. Boyd, Padé approximant algorithm for solving nonlinear ODE boundary value problems on an unbounded domain, Comput. Phys. 11 (1997) 299–303.
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
Simple analytic approximations for the Blasius problem R. Iacono a,∗ , John P. Boyd b a
ENEA, CR Casaccia (SP 118), Via Anguillarese, 301, 00123 Roma, Italy
b
Department of Atmospheric, Oceanic & SpaceScience, University of Michigan, 2455 Hayward Avenue, Ann Arbor MI 48109, United States
highlights • • • •
We derive analytic approximations for the Blasius function f and its derivatives. We extend the integral iteration scheme for the Blasius problem devised by H. Weyl. We compute very accurate bounds for the second derivative of f at the origin. We discuss the new approximations in the context of generalized Padè theory.
article
info
Article history: Received 16 April 2015 Received in revised form 6 August 2015 Accepted 6 August 2015 Available online 13 August 2015 Communicated by M. Vergassola
abstract The classical boundary layer problem formulated by Heinrich Blasius more than a century ago is revisited, with the purpose of deriving simple and accurate analytical approximations to its solution. This is achieved through the combined use of a generalized Padé approach and of an integral iteration scheme devised by Hermann Weyl. The iteration scheme is also used to derive very accurate bounds for the value of the second derivative of the Blasius function at the origin, which plays a crucial role in this problem. © 2015 Elsevier B.V. All rights reserved.
Keywords: Blasius equation Boundary layer Generalized Padé approximants Integral iteration scheme
1. Introduction The Blasius problem is a prototypal example illustrating the development of a boundary layer as a result of the interaction between a weakly viscous flow and an obstacle, a phenomenon that is ubiquitous in nature. By now, it has more than a century of life, during which it has developed a vast bibliography; detailed information about this problem, together with a reminder of the essential literature, can be found in the recent review [1]. Here we only sketch the basic physical picture. Consider a steady, uniform flow impinging on a thin, wide plate parallel to it, and choose a (X , Y ) Cartesian coordinate frame such that the flow points in the positive X direction and the plate is placed along the half-plane Y = 0, X > 0. For X > 0, the stream is basically undisturbed at large Y , but the flow must slow down when approaching the plate, because the fluid velocity must be zero at
∗
Corresponding author. E-mail addresses: [email protected] (R. Iacono), [email protected] (J.P. Boyd). http://dx.doi.org/10.1016/j.physd.2015.08.003 0167-2789/© 2015 Elsevier B.V. All rights reserved.
the solid boundary. In a low viscosity fluid, such as air or water, this only happens in a narrow region near the plate, named ‘‘boundary layer’’, or ‘‘shear layer’’, in which viscosity effects become essential. The main problem then is to compute the flow in this layer, and to match it with the flow outside. Scale analysis of the steady-state Navier–Stokes equations shows that, in a first approximation, pressure gradients can be neglected within the boundary layer, so that the dominant balance, in the longitudinal (X ) direction, is between the nonlinear advection terms and the leading viscosity term,
∂u ∂u ∂ 2u +v =ν , ∂X ∂Y ∂Y 2 where (u, v) are the (X , Y ) components of the flow, and ν is the kinematic viscosity. Together with mass continuity, ∂X u +∂Y v = 0,
u
this equation can be solved to yield the steady flow components in the boundary layer. Alternatively, the equilibrium equations can be reduced to a single third-order, nonlinear partial differential equation (PDE) for the flow streamfunction (see, e.g., [1]). It turns out, however, that the problem can be further simplified, thanks to the existence of a continuous group of invariance,
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
Fig. 1. The Blasius function f (x) [solid black] is plotted with its first derivative [red dashed] and second derivative [blue dots]. The function (x/f (x))df /dx, which has a starring role in later derivations, is the thick green curve with markers.
which also has a physical meaning, being proportional to the shear stress on the plate. The parameter is a priori unknown, and not easy to compute analytically (we will come back to this issue later on). However, a second group invariance of the BE, discovered by Töpfer [3], allows for an efficient numerical computation of κ , which is approximately given by κ ≈ 0.332057336. High precision values of κ and of the other constants that enter the Blasius problem can be found in Table 4.1 of [1]. Although the existence of a complex pole had long been known [4], a detailed understanding of the structure of the singularities of the BE in the complex plane has been achieved much more recently [2]. Among other results, it was found that the three singularities closest to the origin, which determine the radius of the convergence of the series expansion, lie on a circle of radius S, where S ≈ 5.69. It was also shown in [2] that the series can be accelerated using Euler’s summation, and that the accelerated series appears to converge everywhere on the positive x axis, which is the physical domain of the flow. Treating the BE as a linear, first order ODE for fxx , yields the wellknown expression
fxx = κ exp − pointed out by Blasius himself. This implies that the streamfunction is not a function of X and Y separately, but is instead a function of a similarity variable x = Y /X 1/2 . As a consequence, the problem becomes one-dimensional, with the third-order PDE reducing to a third-order, nonlinear ordinary differential equation (ODE), 2fxxx + ffxx = 0,
(1)
which is referred to as the Blasius equation (BE hereafter). Here, f is related to the streamfunction of the two-dimensional problem, and the subscript denotes derivation with respect to the similarity variable x. Eq. (1) needs to be solved subject to boundary conditions both at the origin and at infinity (see [1] for further details): f (0) = fx (0) = 0,
fx (∞) = 1.
(2)
The boundary problem (1)–(2) was solved in [2] using the Chebyshev pseudospectral method. Fig. 1 shows the resulting profiles of f (x) and of its first two derivatives, together with that of the function xfx /f , which will be useful in the developments to follow. This very accurate numerical solution – errors for f and its first three derivatives are smaller than 10−9 over the whole real axis – will be used as a reference in the present work. All the functions plotted in Fig. 1 are smooth and monotonic. The Blasius function f has a simple structure: it grows like x2 at small x, and asymptotes to x plus a negative constant at large x. Despite the simple shape of the solution and the long history of the subject, analytic insight on the BE is still limited. While an accurate approximation for fxx is known (see the next section), to the best of our knowledge, nobody has yet been able to provide a simple analytic approximation for f (x) of comparable accuracy. Constructing such an approximation is the main objective of the present note. 2. Some basic results The series expansion of f about the origin, f (x) ≈
κ 2
x2 −
κ2 240
x5 +
11κ 3 161280
x8 −
5κ 4 4257792
x11 + · · · ,
(3)
was computed by Blasius himself, who derived a recurrence relation for the coefficients of the expansion. The coefficients depend on the parameter
κ ≡ fxx (0),
(4)
73
1 2
x
dξ f (ξ ) ,
(5)
0
which explains the fast decay of fxx with x shown by Fig. 1. Expression (5) also suggests that a guess for f that is good at small x should give an approximation for fxx accurate over a wider x-range, because the main contribution in the rhs of (5) comes from small x. One could place the series expansion (3) in (5), as Bairstow [5] did, but would then be confronted with the fact that the series has a finite radius of convergence. Bairstow noted, however, that taking f ≈ (κ/2)x2 yields
fxxCG ≡ fxx = κ exp −
κ 12
x3 ,
(6)
where CG stands for cubic-Gaussian. The approximation (6) is fairly accurate for x not too large: it is still within 1% from the numerical solution at x = 2.5, whereas the guess for f only preserves a similar accuracy up to x = 1.5. It turns out that Bairstow’s approach can be refined, to build a much better approximation for fxx , fxxPBS ≡ fxx = κ
exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3
,
(7)
which will be used here to test the accuracy of the approximations to be derived. The – non-trivial – derivation of this reference approximation, due to Parlange, Braddock and Sander [6], is sketched in the review [1]. x Finally, we note that integration by parts of the identity f (x) = dξ fξ yields an integral expression for f in terms of its second 0 derivative, f (x) =
x
dξ (x − ξ ) fξ ξ (ξ ).
(8)
0
This suggests an iterative approach based on (5) and (8): given a guess for f that satisfies the boundary conditions at the origin, one could compute fxx through (5), then get a new guess for f using (8), and so on. Would this approach prove convergent, it could be used to derive accurate numerical approximations to the solution of the BE. 3. Weyl’s integral approach The idea of iterating between f and fxx is not new. It goes back to a trio of articles on the BE and related boundary layer problems that Hermann Weyl published three quarters of a century ago.
74
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78 Table 1 Maximum errors in the iterates of Weyl’s integral for the second derivative of the Blasius function, W (z ) = fxx (z (4/κ)1/3 )/κ 1/3 . Chebyshev series with N = 200 on z ∈ [0, 8], but there is no significant change when N = 80 and the Chebyshev approximation interval is reduced to z ∈ [0, 6].
His second contribution [7] is a three-page note that completes the convergence proof for his iteration for the generalized Blasius equation, 2fxxx + ffxx + (k2 − fx2 ) = 0 where k is a parameter (a constant), but says nothing about the BE itself. His final paper [8] amplifies his earlier treatment of the Blasius problem as a prelude to expanding his methodology to similar nonlinear ODEs, but offers nothing new about the solutions to the BE. We therefore shall discuss only his first article [9] where Weyl repeats Töpfer’s group invariance and then continues ‘‘When H. Blasius carried out the integration . . . he made use of the power series for f around z = 0 and of a certain asymptotic expression for large values of z, adjusting the constant β so as to make both expressions dovetail in a middle region. It is a fact, although it seems to have escaped notice, that convergence of the power series of f stops somewhere between √ [converting from√ his convention to the usual definition of f ] 3 9/κ ≈ 3.004 and 3 30/κ ≈ 4.487. Thus, Blasius’ method . . . can yield results of limited accuracy only’’. Weyl was incorrect on a couple of counts. The radius of convergence is S ≈ 5.69, well outside the limits that he gives (see Section 8.2 for more details), and the Euler acceleration of the power series converges on the entire real axis. Nevertheless, his negative evaluation of prior work motivated him to set up an iterative scheme, based on an integral equation for the second derivative of f , which he proved to be convergent. Despite Weyl’s great reputation, his iteration has not been widely used. We extend his treatment in three ways. First, we use a modern computer algebra system (Maple) to extend the analytic formulas for iterates to one order higher than Weyl himself. Second, we show that with modern Chebyshev spectral methods, it is possible to numerically apply his method for arbitrary iteration order at a cost of only O(N log(N )) operations per iteration where N, the number of Chebyshev polynomials, is no larger than O(100) and where the error in the seventh iterate is a maximum of only about 10−14 . Third, we show that the numerical iterates can be used to compute very accurate bounds for the value of κ . Instead of explicitly iterating between (5) and (8), Weyl works with the nonlinear integral equation z
W (z ) = exp −
(z − y)2 W (y)dy ,
(9)
0
where W (z ) = fxx (z (4/κ)1/3 )/κ 1/3 and z ≡ (κ/4)1/3 x, and proves that this equation yields a convergent iteration scheme. The first few iterates are, with the maximum error of the iterate shown in square brackets, W0 = 0 [1.0]
(10)
W1 = 1 [1.0]
(11)
W2 = exp(−[1/3]z ) [0.23]
(12)
z 2 3 W3 = exp − (z − y) exp(−[1/3]y )
(13)
0
= exp(−w(z )) [0.00039]
(14)
where 3
3 3
2 2 1 3 Γ −Γ , z 3
∞
z
exp(−t ) t α−1 dt .
1 0.0235 0.000391 2.89 × 10−6 1.13 × 10−8 2.61 × 10−11 3.84 × 10−14
Taking into account the normalizations, it is easily seen that the second iterate, Eq. (12), is none other than Bairstow’s approximation (6). Although explicit forms for the iterates are not possible for order k > 3, a very fast numerical procedure is possible. First, the exponential of each iterate can be written without approximation as p(z )
z
(z − y)2 Wk−1 (y)dy
=
(17)
0
= z2
z
Wk−1 (y)dy
(18)
0 z
yWk−1 (y)dy +
− 2z 0
z
y2 Wk−1 (y)dy
(19)
0
Wk (z ) = exp {−p(z )} .
(20)
We can approximate each integrand on an arbitrary interval z ∈ [0, b] by a Chebyshev interpolant with N interpolation points. Because of the extremely rapid exponential decay of W (z ) with z, b = 6 gives an approximation accurate to near machine precision (2.2 × 10−16 ) in Matlab. We can then apply the usual recurrence to obtain the coefficients of the indefinite integral from those of the integrand, evaluate the integrand series at each interpolation point zj , and add the contributions of the three integrals to obtain Wk (zj ) = exp(−p(zj )). By using the Fast Fourier Transform for the interpolation and summation, each iterate can be calculated in O(N log2 (N )) operations. Convergence is exponential as shown in Table 1. b Elementary calculus requires that, for arbitrary f (x), a fxx dx = fx (b) − fx (a). The boundary conditions for the derivative of the Blasius function, fx (0) = 0 and fx (∞) = 1, then give the valuable identity ∞
fxx dx = 1.
(21)
3 3
This can be used to approximate κ , the second derivative of the Blasius function at the origin. First, though, we must sort the different normalization used by Weyl.1 Let f denote the Blasius function with Weyl’s scaling so that fzz = W (z ). Then the usual Blasius function f (x) with fxx (0) = κ is connected to f by (22)
It follows that (15)
fxx (x) = κ f′′ ([κ/4]1/3 x)
= κ W ([κ/4]1/3 x).
with Γ (α, z ) the usual incomplete gamma function,
Γ (α, z ) ≡
1 2 3 4 5 6 7
f (x) = 42/3 κ 1/3 f([κ/4]1/3 x).
1 1 3 1 3 −Γ , z w(z ) = 1 − e−1/3 z + z 2 3−2/3 Γ − 2z3
maxz ∈[0,∞] |Wk − W |
0
3
−1/3
Iteration number k
(16)
1 f + 2ff = 0, f (0) = 1 versus 2f + ff = 0, f (0) = κ . zzz zz zz xxx xx xx
(23)
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
75
Table 2 Approximations to κ and relative errors as obtained from the integral of Weyl’s iterates. No answer is given for W1 ≡ 1 because the infinite integral of this iterate is unbounded. Iteration number k
κk
|κ − κk |/κ
2 3 4 5 6 7 8
0.342095321717019 0.331898302513614 0.332058458263831 0.332057331979152 0.332057336224705 0.332057336215182 0.332057336215196
0.03023 0.00047893 3.3791e−06 1.2757e−08 2.8634e−11 4.1961e−14 1.6717e−16
The integral (21) then becomes ∞
κ W ([κ/4]1/3 x)dx = 1,
(24)
0
and, with the change of variable ζ ≡ (κ/4)1/3 x, ∞
W (ζ )dζ = 4−1/3 κ −2/3 .
(25)
0
This allows us to calculate κ from the integral of Weyl’s function as
κk =
1 2
∞ 0
Wk (ζ ) dζ
−3/2
(26)
Table 2 shows that the resulting approximations converge very rapidly. The eighth iterate is correct to all fifteen decimal places shown in the table. It should be noted that, since Weyl’s iteration produces alternating bounds for W (z ), (26) yields alternating bounds for κ . Only the first upper bound (iteration number 2) was computed in [9]. 4. New analytic approximations for f and fxx An approach often used to extend the accuracy of the series expansions beyond their convergence range is that of the Padé approximants, which are rational approximations constructed from the series (see, e.g., [10]). In the context of the Blasius problem, this approach was pursued in [11], in the attempt of constructing approximations for f that would also allow the determination of the parameter κ . Because of the structure of (3), a Padé approximant for f should be constructed as x2 times the ratio of two polynomials in x3 , but such an approximant cannot capture the large-x behavior, which is linear in x. High order, diagonal Padé approximants of this form were computed in [11] that proved very accurate at small and intermediate x, but no simple way was found to extend this accuracy to larger values of x. Here we tackle the problem from a different angle. Since it is difficult to construct good rational approximations for f , we may try to do so for some other variable that is best suited. A good candidate seems to be the variable y ≡ xfx /f that we plotted in Fig. 1. This is a smooth, decreasing function that equals 2 at the origin, tends towards 1 at large x, and is a function of x3 at small x. The simplest rational approximation for y is therefore of the form y ≈ (2 + α x3 )/(1 + α x3 ), with α a constant that can be determined by enforcing the x3 term in the series expansion of y about the origin. This yields y=
2 + (κ/40)x3 1 + (κ/40)x3
,
(27)
which is very accurate at small x, and is still within 1% from the numerical y(x) at x = 4.8. Upon integrating, we get a simple analytic approximation for f : fA =
κ x2 /2 . [1 + (κ/40)x3 ]1/3
(28)
Fig. 2. Absolute errors in approximations to fxx , the second derivative of the Blasius function. The thin solid red curve is the Parlange–Braddock–Sander approximation, fxxPBS ≡ κ exp(−0.0125152x3 )/(1 + 0.00505209x3 )3 . The dashed κ 3 2/3 black curve is fxxa ≡ κ exp −5 (1 + 40 x ) − 1 . The thick black curve is fxxb ≡
m A exp − 4p (1 + px3 )2/3 −
n 2q
(1 + qx3 )1/3 , with the coefficients given in Eqs. (43)
and (45).
This is already a pretty accurate approximation for small and intermediate values of x, with relative errors of 2.6 × 10−7 at x = 1, 4.4 × 10−4 at x = 4, and less than 3% at x = 8. Then it loses accuracy, because it does not have the correct asymptotic dependence, even though it captures the linear growth with x. Based on the consideration of the previous sections, we would expect that placing (28) in (5) should give a good approximation for fxx . This is indeed the case. Integration of (5) gives
fxx ≡ fxxa ≈ κ exp −5 (1 +
κ
40
x3 )2/3 − 1
,
(29)
which yields the absolute errors plotted as a dashed line in Fig. 2. The errors for (29) are smaller than those for the reference approximation (7) (red line) up to x ≃ 6, and the maximum error, 3.7816 × 10−5 , is smaller than the maximum error for (7), which is of 1.0176 × 10−4 . Thus, quite surprisingly, we have found a simple analytic approximation for fxx that is superior to the sophisticated approximation by Parlange et al. [6]. Placing (29) in (8) gives a new approximation for f : fB = κ
x
dξ (x − ξ ) exp −5
1+
0
κ 40
ξ3
2/3
−1 .
(30)
This time, the integral cannot be evaluated analytically, but f can be easily computed numerically. As shown in Fig. 3, this yields a very accurate approximation (dashed black curve), with relative errors smaller than 10−5 up to x ≃ 6, and still smaller than 10−4 up to x ≃ 100. Errors are much smaller than those for the approximation (28) (solid black curve), confirming that the iterative scheme based on (5) and (8) indeed works well. 5. Improving the small x behavior The form of (28) suggests a simple generalization. Define a function g so as to strip off the leading quadratic factor of f , enabling us to work with a function of z ≡ x3 only: f =
κ 2
x2 g (z )
↔ g (z ) =
2
κ x2
f (x).
(31)
76
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
Fig. 3. Relative errors in approximations to f , the Blasius function. Solid black: f A ≡ (κ/2)x2 /(1 + (κ/40)x3 )1/3 . Triangles are for the asymptotic approximation f asy,best , defined in Eq. (49).
Then, (28) can be considered as the lowest order [N /(N + 1)] Padé approximant to the cube of g, g 3 , i.e.,
g[0/1] =
1/3
1 1 + (κ/40)z
.
(32)
Higher order approximants are easily computed: 11200(396 + 43 κ z )
g[1/2] =
1/3
4435200 + 592480 κ z + 12073 κ 2 z 2
g[2/3] =
1663200
P2
1/3
Q3
,
(33)
,
(34)
with
Fig. 4. Errors in approximations to g (z ) which are each the cube root of a Padé approximation to g 3 (z ). The error in the first two orders of the large z-approximation is also shown (dashed line). The Padé orders label each curve: degree of the numerator polynomial in z on the left, denominator degree to the right of the slash. The horizontal dotted guideline is at 10−5 . The vertical guideline is at the radius of convergence of the power series.
illustrates the accuracy gain at small and intermediate values of x obtained by increasing the order of the approximant. 6. Improving the large x behavior The approximations derived in the previous sections are quite good, but could still be improved at large x, were we able to enforce the correct asymptotic behavior. Consider, for example, the lowest order approximation fA [Eq. (28)]. Since f ∼ B + x at large x, with B ≈ −1.7207876575, we can try to improve the large x behavior by adding a term with the same structure of the one in fA , but with a 2/3 power in the denominator, which would tend towards a constant in the limit x → ∞. That is, we may seek an approximation of the form fC =
P2 = 3964295171200 + 76608724296 κ z
+ 507276173 κ 2 z 2 ,
(35)
(1 + px3 )1/3
+ 4078150639904880 κ z + 16713554600203 κ z , (36) 3 3
g[3/4] =
9828 P3 5
Q4
,
1/3 (37)
nx2
(1 + qx3 )2/3
.
(40)
n/q2/3 = B,
(41)
which gives the correct large x behavior, and m + n = κ/2,
and
+
We determine the coefficients in (40) by asking that m/p1/3 = 1,
Q3 = 6593415728739840000 + 292251023467603200 κ z 2 2
mx2
mp + 2nq = κ 2 /80,
which enforces the first two terms of the series expansion about the origin. Solving the system (41)–(42) with Mathematica, we find a single real solution, given by
with
m = 0.198177262,
n = −0.0321485939,
P3 = 52740459296753088904184064000 + 1156660146147189506418534400 κ z
p = 0.0077832588,
q = 0.0025535952.
+ 9352330559888045666104080 κ 2 z 2 + 11843596516363041045577 κ 3 z 3 ,
(38)
Q4 = 103666646793697871550064196198400 + 4865197353109362482567876121600 κ z
+ 75992549747403848831924652000 κ 2 z 2 + 416014006359213902061551220 κ 3 z 3 + 401763679454602196561927 κ 4 z 4 .
Recalling that κ = 0.33205733621519630, we may compute the relative errors (absolute errors divided by g) for these approximations, which are displayed in Fig. 4. The figure clearly
(43)
The relative errors for this approximation are shown as red dots in Fig. 3. The approximation is slightly less accurate than fA [Eq. (28)] at small x, but much more accurate at large x, with the largest relative errors, at intermediate values of x, that remain smaller than 1%. Placing (40) in (5) gives a new analytic approximation for the second derivative of f ,
(39)
(42)
fxxb = A exp −
m 4p
(1 + px3 )2/3 −
n 2q
(1 + qx3 )1/3 ,
(44)
where A = κ exp
m 4p
+
n 2q
≈ 0.356393092.
(45)
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
The absolute errors in this approximation are graphed as a thick black curve in Fig. 2. The figure shows that (44) is definitely more accurate than the reference Parlange–Braddock–Sander approximation, with a maximum error that is about 1/8 of that for fxxPBS . Finally, we can place (44) in (8) to get a new approximation for f
8.2. Recurrence relation for power series In (4.2) of [1], the right-hand side should be divided by a factor of two. The correct power series recurrence is am = −
d ξ (x − ξ )
fD = A
×
0
2 (3m − 1)(3m − 2)(3m − 3)
(3j − 1) (3j − 2) aj am−j
(51)
j =1
m n × exp − (1 + pξ 3 )2/3 − (1 + qξ 3 )1/3 . 4p
2q
(46)
This yields the relative error curve marked by stars in Fig. 3, with a maximum error just above 10−5 at x = 100. 7. Improved large x asymptotics
fxx (x) ∼ Q exp(−(1/4)(x + B) ). 2
2 (47)
This is derived by substituting f ∼ B + x into 2fxxx + ffxx = 0 and solving as a first order ODE for fxx with the initial condition fxx (0) = κ . It turns out that Q = 0.2337276212850634. Integrating this twice, both with lower limits at zero, gives
− 2 exp(−(1/4)B2 ) + 2 exp(−(1/4)(x + B)2 ) .
1 120
f int ∼ −1.48925285928230415948146 (48)
We can therefore get a better asymptotic approximation by writing f asy,best ∼ (B + 1.48925285928230415948)
2
1
n (52)
36
exp(−(1 − 3β1 )y)
(53)
(1 + β1 y)3 exp(−0.452277442494834y)
(54)
(1 + 0.182574185835055y)3 exp(−(1 − 3β1 )(κ/12)x3 )
(55)
(1 + [β1 κ/12]x3 )3 exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3
.
(56)
Eq. (5.12) of [1] is in error, writing 0.4114 in the exponential instead of the correct 0.45227. However, (5.17) in the same article is correct.
− 2 exp(−(1/4)B2 ) (49)
The errors of (49) are graphed as triangles in Fig. 3. The approximation is superb for large x but poor for small x, and thus complementary to most of the other approximations derived here. 8. Corrections We are aware of several previously published typographical errors. In the spirit that repair is as important as new construction, we correct these errors in this section even though these errors are not directly related to any of our new approximations. 8.1. Large x approximation typo
9. Conclusion Padé approximants have proved so useful that they have been greatly generalized. Although we have not invoked generalized Padé theory in our derivations, we would be amiss in failing to note that many of our approximations fit in this generalized framework. A standard Padé approximant is the solution to the linear polynomial equation QN (z )fM .N − PM (z ) = 0 where the coefficient polynomials PM and QN are chosen so that the power series expansion of f[M .N ] agrees with the leading terms of f (z ). In his thesis, Padé – using ideas developed two decades earlier by his thesis advisor, Hermite – observed that there was no reason to restrict the polynomial equation to first degree. For example, one could define an approximation to f (x) by declaring that g[N ,K ,L,M ] is the solution to the cubic polynomial equation
Q(z )g[3N ,K ,L,M ] (z ) + R(z )g[2N ,K ,L,M ] + Sg[N ,K ,L,M ] − P(z ) = 0
The lowest order asymptotic approximation for large x is fxx ∼ exp(−B2 /4)Q exp(−(1/4)x(x + 2B)).
fxx = κ fxx = κ
+ (1 − .735875425275434309752)x √ π (B + x) [erf(x/2 + B/2) − erf(B/2)] +Q
x ≫ 1.
≤ cn ≤
1
8.3. PBS approximation
=κ
+ 0.735875425275434309752792 x.
n
where both the left and right hand sides are smaller by 4−n than the incorrect bounds in Weyl [9].
fxx = κ
At large x, this gives
Parlange, Braddock and Sanders give the approximation, using y = (κ/12)x3 ,
√ π (B + x) [erf(x/2 + B/2) − erf(B/2)]
+ 2 exp(−(1/4)(x + B)2 )
where the factor erroneously omitted in [1] is boxed. Weyl erroneously placed this same factor of two in the numerator instead of the denominator in [9], pg. 579. His power series coefficients cn [the coefficient of z 3n+2 ] are therefore too large by 4n . His error shrinks the radius of convergence (noting that the series proceeds in powers of z 3 ) by 41/3 (to 3.58 in our notation). His bounds on his coefficients cn should be 1
The asymptotic behavior of fxx is given by
f int ∼ Q
1
m−1
x
77
(50)
The exponential exp(−B2 /4) was incorrectly omitted from (3.2) of [2] and the error repeated in [1].
(57)
and where the coefficients in the cubic equation are polynomials in z only of degrees (N , K , L, M ) respectively. Standard Padé approximants to g (= f /(κ x2 /2)) would not be able to capture the decay of g as z 1/3 at large z. We therefore computed Padé approximants to the cube root of g. This is
78
R. Iacono, J.P. Boyd / Physica D 310 (2015) 72–78
equivalent to a cubic Hermite–Padé approximation for the special case K = L = 0 so that only the constant and the cubic term are nonzero. QN (z )g[3N ,0,0,M ] (z )] − PM (z ) = 0
↔ g[3N ,0,0,M ] (z )] =
PM QN
. (58)
Unfortunately, solutions to polynomial equations have only algebraic branch point singularities. More complex behavior can be obtained by replacing the polynomial equation by a simple differential equation. For example, one species of ‘‘Padédifferential approximant’’ can be defined as the solution of the first-order linear differential equation
A = κ exp
m
+
4p
n
≈ 0.356393092.
2q
(A.5)
f asy,best ∼ (B + 1.48925285928230415948) + (1√ − .735875425275434309752)x +Q π (B + x) [erf(x/2 + B/2) − erf(B/2)]
− 2 exp(−(1/4)B2 ) + 2 exp(−(1/4)(x + B)2 ) ,
x ≫ 1.
(A.6)
A.2. Approximations to the second derivative
x q(y) dy u = exp −
ux + q(x)u = 0
where
(59)
0
where q(x) is chosen to be a power series whose coefficients are such that the leading terms of the power series for u(x) match those of the function which is to be approximated by u(x). This is exactly what we have done with u ≡ fxx and q(x) ≡ 2f (x). Acknowledgments
fxxCG = κ exp − fxxPBS = κ
κ 12
x3 .
exp(−0.012515170232088x3 )
(1 + 0.005052091484171x3 )3 κ 3 2/3 fxxa = κ exp −5 (1 + x ) −1 . fxxb
Appendix. Summary of approximations
with A given by (A.5).
A.1. Approximations to the Blasius function
References
fA =
κ x /2 . [1 + (κ/40)x3 ]1/3
(A.1)
[fA is equivalent to (κ/2)x2 g[0/1] ] fB = κ
x
dξ (x − ξ ) exp −5 0
mx2
1+
κ 40
nx2
ξ3
2/3
−1 .
, (1 + px3 )1/3 (1 + qx3 )2/3 m = 0.198177262, n = −0.0321485939, p = 0.0077832588, q = 0.0025535952. x fD = A dξ (x − ξ ) 0 m n × exp − (1 + pξ 3 )2/3 − (1 + qξ 3 )1/3 fC =
+
4p
2q
(A.2)
(A.3)
(A.4)
.
40
We thank Tyler Nielson and his project group at Eastern Washington University for finding the error in the power series recurrence (4.2) of [1].
2
(A.7)
n m 3 2/3 3 1/3 = A exp − (1 + px ) − (1 + qx ) 4p
2q
(A.8) (A.9)
(A.10)
[1] J.P. Boyd, Themes illustrated by the blasius function: Numerical computations before computers, the value of narrow tricks, and interesting undergraduate projects and open research problems, SIAM Rev. 50 (2008) 791–804. [2] J.P. Boyd, The Blasius function in the complex plane, J. Exp. Math. 8 (1999) 381–394. [3] K. Töpfer, Bemerkung zu dem Aufsatz von H. Blasius Grenzschichten in Flüssigkeiten mit kleiner Reibung, Zeit. Math. Phys. 60 (1912) 397–398. [4] B. Punnis, Zur Differntialgleichung der Plattengrenzschicht von Blasius, Arch. d. Math. 7 (1956) 165–171. [5] L. Bairstow, J. R. Aero. Soc. 19 (1925) 3. [6] J. Parlange, R.D. Braddock, G. Sander, Analytical approximations to the solution of the Blasius equation, Acta Mech. 38 (1981) 119–125. [7] H. Weyl, Concerning the differential equations of some boundary layer problems, Proc. Natl. Acad. Sci. 28 (1942) 100–101. [8] H. Weyl, On the differential equations of the simplest boundary layer problems, Ann. of Math. 43 (1942) 381–407. [9] H. Weyl, Concerning the differential equations of some boundary layer problems, Proc. Natl. Acad. Sci. 27 (1941) 578–583. [10] G.A. Baker Jr., Padé approximants, in: Advances in Theoretical Physics, in: Advances in Theoretical Physics, vol. 1, Academic Press, New York, 1965, pp. 1–50. [11] J.P. Boyd, Padé approximant algorithm for solving nonlinear ODE boundary value problems on an unbounded domain, Comput. Phys. 11 (1997) 299–303.