Mathematical Problems in Engineering

Volume 2010 (2010), Article ID 874540, 23 pages

http://dx.doi.org/10.1155/2010/874540

## Approximate Ad Hoc Parametric Solutions for Nonlinear First-Order PDEs Governing Two-Dimensional Steady Vector Fields

Department of Engineering Sciences, University of Patras, GR 26504, Greece

Received 20 April 2010; Accepted 3 November 2010

Academic Editor: Oleg V. Gendelman

Copyright © 2010 M. P. Markakis. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Through a suitable ad hoc assumption, a nonlinear PDE governing a three-dimensional weak, irrotational, steady vector field is reduced to a system of two nonlinear ODEs: the first of which corresponds to the two-dimensional case, while the second involves also the third field component. By using several analytical tools as well as linear approximations based on the weakness of the field, the first equation is transformed to an Abel differential equation which is solved parametrically. Thus, we obtain the two components of the field as explicit functions of a parameter. The derived solution is applied to the two-dimensional small perturbation frictionless flow past solid surfaces with either sinusoidal or parabolic geometry, where the plane velocities are evaluated over the body's surface in the case of a subsonic flow.

#### 1. Introduction

First-order PDEs, which mostly appear in fluid mechanics, describe the motion of ideal as well as of real fluids [1–3] and govern even the electrostatic plasma oscillation [4]. As is well known, there is no complete general theory concerning the derivation of exact analytical solutions for such equations. However, general solutions can be obtained for the quasilinear forms by means of the subsidiary Lagrange equations ([1, Section 2.6.a], Appendix A). We also mention Charpit’s method for the general nonlinear case that yields to complete and general solutions [1, Section 2.6.b]. These solutions involve arbitrary functions of specific expressions of the dependent and independent variables. Furthermore, appropriate transformations of the dependent and (or) independent variables [1, Section 2.1], combined in several cases with the introduction of auxiliary functions (like stream functions), can occasionally linearize the original equation or more generally reduce it to a solvable form, like a quasilinear one, or even to a nonlinear ODE.

In our previous work [5], four simplified forms of the full two-dimensional nonlinear steady small perturbation equation in fluid mechanics [6] were treated analytically. As far as the three of the considered cases are concerned, closed form solutions have been derived for the two dependent variables of the equation, which represent the dimensionless velocities , of a perturbed frictionless flow past a solid body surface, while in the fourth case, a parametric solution was obtained with regard to these velocity resultants. We note that the components , are parallel to the , axes of the Cartesian plane, respectively (see Figure 1 in Section 4, where a wavy surface is represented), with being the direction of the uniform velocity of the steady flow. The extracted closed form solutions provide as a specific expression of , as well as an equation for involving an unknown arbitrary function. The analytical method was based on the introduction of a convenient *ad hoc assumption*, originally due to Pai [7], by means of which the original (simplified) equations, as well as the irrotanional relation, take a quasilinear form integrated by the Lagrange method. Thus, the above-mentioned solution (including the unknown function) for is obtained, together with an ordinary differential equation, which, after a further analytical treatment, provides the exact or approximate (depending on the case) solutions . However, it should be mentioned that only in the first, more simplified case [5, Equation (9)] of the general equation, the unknown function can be defined by the use of the boundary condition of the problem, resulting in a transcendental equation for (or ). Furthermore, no investigation has been performed in [5] with regard to the effectiveness of the obtained formulas (the expressions extracted in the application [5, Section 5] concerning the above-mentioned simplified case and the parametric solutions derived for one of the other examined cases [5, Equation (8)]) to evaluate the perturbed flow field.

In the present work (Section 2), we firstly treat a steady three-dimensional PDE concerning a general weak irrotational vector field. By taking into account the three irrotationality conditions and using the *ad hoc assumption* introduced in [5], the Lagrange method (see Appendix A) finally results in a system of two nonlinear ODEs for the two unknown functions introduced by the *ad hoc assumption.* These functions represent the field’s components , , while the first component stands for the independent variable. In Section 3, we proceed into the integration of the first ODE, which corresponds to the plane problem (the second involves also ). The herein developed methodology consists of a *functional transformation* of the dependent variable, in combination with an appropriate *split* of the resulting equation by using an arbitrary function, which eventually is eliminated. By this technique, we finally derive an Abel equation, which admits a parametric solution. Thus, we obtain the field’s components and as explicit expressions of a parameter . In several steps of the analysis developed in Section 3, the established, in Appendix C (linear), approximations based on the weakness of the field have been used. Additionally, some limitations imposed by the analysis (see Cases P-1, P-2 in Appendix D) affect the domain of the physical parameter(s) of the problem, for which the extracted solution is valid.

Then in Section 4 we apply the obtained parametric solution in the plane case of the full small perturbation equation, simplified forms of which were investigated in [5]. Here, by combining the extracted parametric formulae with the boundary condition concerning the flow tangential to the solid surface, a transcendental equation is derived, involving , , , where , represent the plane coordinates on the body’s surface. Then, for a given pair , the solution of this equation yields , and hence the “surface” perturbed flow velocity field , can be evaluated (the perturbed velocity components , refer to the , cartesian plane). Moreover, by expanding in Taylor series and taking into account the small perturbation, the perturbed velocities can be approximately obtained within a thin zone over the surface. In addition, under the mentioned limitations, we deduce that the obtained results hold true for subsonic flows as well.

Finally, by means of the extracted formulas, graphic representations of the perturbed field versus are obtained, concerning a sinusoidal as well as a parabolic boundary, and the results are compared to the solution of the linearized equation.

#### 2. The Analytical Procedure

##### 2.1. Transformation of the Governing Equations

Consider an irrotational field satisfying the following PDE: where summation convention has been adopted and with being the Cartesian space coordinates. Equation (2.1) is assumed dimensionless and properly scaled, while the coefficients (with the respective upper and subindexes) represent constants or functions of one or more parameters. In this paper, we investigate the case where as well as the case whereHowever, the proposed solution can also be applied to cases where the coefficients involved in (2.3a) and (2.3b) are sufficiently small, so that the respective terms of (2.1) can be neglected in comparison with the others. Moreover, the field is supposed to be weak in the plane, that is, In fact the approximations (see Appendix C), used in certain steps of the analytical procedure, are based on the weakness of the field under consideration.

As a first step, we make the *ad hoc *assumption that the components and are functions of the component , namely,
and thus by substituting (2.5), (2.1) (taking into account (2.3a) and (2.3b)) becomes
where has been replaced by andHere, the prime “ ' ” denotes differentiation with respect to .

On the other hand, the irrotational condition of the field is written in the form where is the well-known Levi-Civita tensor and represent the unit vectors corresponding to , , respectively. By substituting the assumption (2.5) into (2.8), we arrive at the following three equations ( is replaced by ):

With respect to the physical relevance of (2.1), as well as of the constraints imposed above, we note the following. No “mixed” nonlinear terms involving the plane components , together with are included in (2.1). Furthermore the restrictions (2.3a) and (2.3b) focus on cases where specific nonlinear terms are involved into the governing equation. More precisely, the procedure developed in this paper confronts nonlinear equations where the partial derivatives of the field components appear in products together with specific combinations of these components, of the first and the second degree. Indeed by (2.3a) and (2.3b), it is obvious that two groups of nonlinear terms are formed with respect to the variations of the plane components , , along their own axes () and the other axis (, ). This can be clearly observed in the two-dimensional steady small perturbation equation of fluid mechanics, treated in Section 4 (4.1) as an application of the present analysis.

All the above notations, as well as the *ad hoc* assumption (2.5), outline a normalized structure as regards the behavior of the field in phase space, due to a regulated physical setup. In fact the small perturbation (4.1) is representative of the imposed restrictions, since the origin of the field (the perturbed velocities due to slight “geometric perturbations” of the body’s surface) combined with the orientation of the uniform flow (with reference to the body—see Figure 1 in Section 4) can give rise to the specific nonlinear form of the governing equations (2.1), (2.3a), and (2.3b), as well as to the “weakness” and the *ad hoc* assumptions, (2.4) and (2.5), respectively.

##### 2.2. Construction of Intermediate Integrals

Now, by integrating the correspondent to (2.6), (2.9a), (2.9b), and (2.9c) subsidiary Lagrange equations (see Appendix A), we, respectively, obtain the following general solutions:where , , , and are arbitrary functions possessing continuous partial derivatives with respect to their arguments.

##### 2.3. Reduction to a System of Nonlinear ODEs

In view of (2.10), (2.11a),(2.11b), and (2.11c), we construct a first set of relations by equating identically the functions , , , and as well as their arguments. Thus, excluding the cases where in the extracted equations: we eventually obtain the following systems.

*Case 1 (). *We have

*Case 2 (). *We have

*Case 3 (). *We have

*Case 4 (). **Subcase 1 (). *We have
*Subcase 2 (). *We have

Subcases 1 and 2 are, respectively, derived by equating the arguments of and in two possible combinations. Then, in order to obtain a system of equations not containing explicitly , , and , we find that Cases 1, 3 and Subcase 2 are compatible to each other. Thus by combining their respective equations, we derive the following ODEs: Taking into account (2.7a), (2.7b), and (2.7c), we note that (2.18) contains only and , and thus it constitutes the main equation, the manipulation of which is presented in the next section.

Therefore, the ordinary differential equations (2.18) and (2.19) represent the reduced forms of the partial differential equations (2.6), (2.9a), (2.9b), and (2.9c), via assumption (2.5). Then by substituting (2.7a), (2.7b)) and replacing with and with , (2.18) becomes where denotes the derivative of with respect to and withWe note that and as well as all the other coefficients appearing in the next sessions are listed in Appendix E. Henceforth, the prime will denote differentiation with respect to the corresponding suffix.

#### 3. Integration of (2.20)

##### 3.1. Transformation of (2.20)

Introducing transformation the left hand side of (2.20) results in a nonlinear expression involving , , , , and . Thus, by taking into account this expression and setting with , as in (2.21), (2.20) takes the form Then, by substituting (3.4) becomes In addition, by substitution of (2.21) and (2.22b) into (3.3), we obtain with as in (2.22b).

##### 3.2. The Split of (3.6)

We now split (3.6) into the following system of equations: where is an unknown arbitrary function. Furthermore, after dividing (3.8b) by and setting(3.8b) will be written as where represents now the unknown arbitrary function. We see that (3.10) is an Abel equation of the second kind, and thus by following the analysis presented in [8, Chapter 1, Section 3.4] and taking into account (3.2) and (3.3), it is reduced to a simpler Abel equation, namely, with now, by differentiating , given by (3.12), with respect to and using (3.2), (3.3) (we consider the appropriate domains where is invertible and hence ) as well as the expression of provided by (3.12), we obtain , substitution of which into the left-hand side of (3.8a) results in where (3.9) has been substituted for .

Equations (3.11) and (3.13) form a new system equivalent to that of (3.8a), (3.8b), obtained by splitting (3.6). The elimination of the arbitrary function yields a nonlinear ODE, which represents the reduced form of (3.6). More precisely, after some algebra we extract Furthermore, as far as the -term is concerned, combination of (3.12) with (3.1) and (3.5) yields , . By differentiating with respect to and taking into account certain relations obtained above, as well as that , , and represent and , respectively, we conclude that is equal to , where represent expressions of the equation’s coefficients. Therefore, when the plane field’s components, as well as the variation of with respect to , are very small compared with the unit (e.g., if they denote perturbed components in a small perturbation theory), we can perfectly consider and hence we can neglect the term in the left-hand side of (3.14) in comparison with the others, as it is of . We should note here that in our previous work [5], after following a different analysis concerning two simplified forms of the full equation, an analogous to (3.15), but weaker approximation, has been applied, since the neglected term was of , yielding less accurate results compared to the obtained herein solution of (3.14) especially when takes smaller values than . Moreover by means of (3.3) and (3.12), we have Thus, by writing neglecting the term and multiplying with , then by using (3.16), (3.14) becomes The above equation is also an Abel equation of the second kind and thus we proceed as in [8, Chapter 1, Section 3.4]. More precisely, by using the formulas (D.10a), (D.10b) (see Appendix D), after some algebra, we arrive at where and being as in (2.22a). Now, by applying (C.4) (Appendix C) to both rational functions in the right-hand sides of (3.19) and (3.20), we obtain Finally, substitution of into (3.21) yields Moreover, by the followed procedure (see [8]), we have that with as in (2.22b) and given by (D.10b) (see Appendix D). Finally, the Abel equation (3.23) is solved parametrically (Appendix B, formulas (B.7)) as where In the above relations represents the parameter while is an arbitrary constant. Now, by substituting and taking into account (3.22), the above parametric solution takes the form All the coefficients appearing through the analysis are listed in Appendix E.

##### 3.3. The Parametric Solution for the Field’s Components ,

By combining (3.1), (3.5), (3.12), (3.17), and (3.24), we obtain Approximating linearly , namely, and substituting from (3.28), as well as and from (D.10a), (D.10b) (Appendix D), then (3.29) yields Furthermore, by solving the first part of (3.28) for , we have Now, approximating linearly the powers of involved into (3.31), that is and substituting (3.32) and (3.33) into (3.31), then by replacing with and with and taking also into account (3.28), we conclude thatwith as in (3.27) and given by (2.22a). Equations (3.34a), (3.34b) constitute the approximate analytical parametric solution of the problem for , . As far as the component is concerned, combination of (2.18) and (2.19) results in The above equation can be simplified a little if we neglect the last term in the left-hand side (it is of the form ) by considering . Anyhow we will not investigate (3.35) in this work.

Moreover, in order to evaluate the constant involved into the parametric solution (3.34a), (3.34b), we need a boundary condition, that is, to locate to a point where the field components , are known. Then, by solving (3.34b)) for , we extract the corresponding value of the parameter , and finally, by using (3.34a), we arrive at

In the next section we apply the derived solution in the two-dimensional case of a flow past bodies with specific boundaries.

#### 4. Parametric Solution for a 2-D Flow

As an application of the parametric solution obtained above for the plane case of (2.1), we consider the full nonlinear PDE governing the two-dimensional , steady small perturbation frictionless flow past a solid body surface [6], namely, where is the correspondent to the uniform flow Mach number, which stands for the physical parameter of the problem, and is the ratio of the specific heats usually taken equal to 1.4; hence, the respective (dimensionless) coefficients , and of (2.1) (the term vanishes) are given by Relations (2.3a), (2.3b) also hold true. As mentioned in Section 2, the above equation represents a highly appropriate case, where the physical relevance of the imposed constraints (2.3a)-(2.3b)–(2.5) can be explained by a normalized physical background like the one generated by a uniform flow passing over a slightly “perturbed” surface, according to a specific geometry (see Figure 1, and the applications at the end of this section). Here, , represent the dimensionless perturbation velocity components along the , axes (see Figure 1), normalized by the uniform velocity of the steady flow, which is parallel to the direction in the physical plane.

A wavy surface (projection in the plane) is shown in Figure 1, as a representative case able to produce small plane perturbations in the velocity field (the surface is supposed to have very small amplitude). Moreover, the irrotationality condition (2.9c) holds true and (2.10) and (2.11c) become respectively, where , denote arbitrary functions and , are as in (2.7a),(2.7b)) (we mention that . Obviously in the two-dimensional case, (2.9a),(2.9b)) and (2.11a),(2.11b)) become identities. Comparison between (4.3) and (4.4) results in (2.18).

If we refer now to the proper conditions restricted by the analysis (see Appendix D), we extract that the discriminant of ( has been replaced by ) is always positive (), and, moreover, since , by obtaining the roots of , considering the respective to the Cases P-1 and P-2 intervals for and assuming as well, then a restriction to the domain of is derived. More precisely, we find that formulae (3.34a), (3.34b) are valid for Thus, for the specific 2-D steady flow field, the obtained approximate solution can be applied only to subsonic flows. Moreover, as far as the integral , involved into the function , is concerned (see (3.26), (3.27)), the discriminant of is evaluated negative (), and therefore the integral is obtained as

Now, in order to construct an appropriate procedure to obtain , we consider the well-known boundary condition (see [6, page 208])
where is the total dimensionless velocity vector of the flow at the solid surface, while represents the equation of the *“surface line”*, that is, the section of the body’s surface with the plane. Here, , denote the plane coordinates on this line with , being the body’s length, and . Condition (4.7) states that at the surface of the body the direction of the flow must be tangential to the surface line. Developing (4.7), we arrive at
where by neglecting , we obtain
By squaring (4.9) and substituting and by their parametric expressions (3.34b),(3.34a), we derive a transcendental equation for , namely,
Thus, for a given pair on the *surface line*, the solution of (4.10) results in , substitution of which into (3.34a), (3.34b) yields the perturbed velocity vector of the flow at . In fact, in the case of the flow under consideration, only the perturbed velocity is evaluated by, use of the extracted parametric solution, since due to (4.9) simply expresses approximately the slope of the *surface line.* Furthermore, assuming that the functions , are analytic inside a domain located on any line constant with ( represents the body’s length) and , slightly different from , by developing in Taylor series around , we have
Taking into account the small perturbation theory (the derivatives involved into the series (4.11), as well as , are very small compared to unity) and also that lies close enough to , so that , all the terms after the first in the right-hand side of (4.11) can be neglected. Thus, we can approximately evaluate the perturbed plane flow field inside a thin zone over the body’s surface. Obviously, the thickness of this zone depends on the order of magnitude of . For example, if the boundary has a sinusoidal shape (one of the cases considered below), that is, , , and if we take , then within the domain (a plane zone of thickness (measured in the -direction) with parallel sinusoidal boundaries), the error in (4.11) is . Therefore, the above approximation is valid inside a zone over the solid surface of thickness less or equal to . In addition, in order to obtain and (see the end of Section 3), the axes origin is used which is located at the point where the flow arrives at the body surface and consequently , , , where is given by (4.9). Therefore, by means of (3.34b)) and (3.36) we conclude that
with provided by (3.26) and (3.27), where the integral is evaluated by (4.6).

The derived solution, constructed by relations (3.34a), (3.34b), and (4.10), is applied to the two-dimensional steady frictionless flow past a boundary of sinusoidal (wavy wall), as well as of a parabolic shape. The problem is governed by (4.1). Especially for the “sinusoidal” boundary problem, implicit solutions in the form of transcendental equations have been extracted in [5, Section 5, Equation (72) and (74), (77)], for the more simplified case of (2.1), where the only nonzero coefficients were , , and [5, Equation (9)]. Here, boundary condition (4.7) holds withwhere and describe the sinusoidal and parabolic form of the surface, respectively, while and denote the amplitude and the curvature of the *surface line* in the cases under consideration. The low magnitude of and the large magnitude of allow the small perturbation theory to be applied. Additionally, in both (4.13a), (4.13b), we have , where stands for the assumed body’s length, while in the “sinusoidal” case, the wavelength of the wavy surface is equal to .

As far as the graphs exhibited below are concerned, the “dashed” line represents the sinusoidal or the parabolic boundary, with geometries: , (Figures 2(a), 2(c)—(4.13a)) and , (Figures 3(a), 3(b), (4.13b)). Moreover, the solid blue line in Figures 2(b) and 3(a) has been obtained as the solution for of the linearized form of (4.1), where the slope of the solid surface has been substituted for the component . In both geometries, the body’s length is taken equal to (three wavelengths in the wavy case) and the correspondent to the uniform unperturbed flow Mach number is set equal to 0.7. We note that by changing the values of the geometric parameters involved in (4.13a) and (4.13b), as well as the value of the Mach number, the perturbed field presents qualitatively similar graphs to those obtained here. Finally, as mentioned above the perturbed velocity is obtained as the slope of the surface.

In Figure 2(b), we note that the linear approximation is excellent throughout the body’s length except in small intervals centered at the picks of the sinusoidal surface with radius approximately equal to . Outside these locations the maximum error of the linear approximation (with respect to the ad hoc solution) is approximately equal to , while inside these intervals the difference between the two solutions increases with moving towards the pick. Furthermore, concerning the comparison of the solutions in the case of the parabolic surface (Figure 3(a)), we find that for the considered body’s length, the maximum error of the linear approximation is approximately equal to (the error increasing with ).

#### 5. Summary and Conclusion

In this paper an *ad hoc *analytical parametric solution has been obtained, concerning a nonlinear PDE governing a two-dimensional steady irrotational vector field. However, in Section 2 of this work the three-dimensional case is treated. As a result, we obtain a system of two (nonlinear) ODEs being equivalent to that of the original PDEs (including the irrotationality conditions). The analytical tools have been used in order to integrate the first ODE (concerning the two-dimensional case), in combination with linear approximations of certain polynomial and rational expressions, succeeded in transforming the above equation to a parametrically solvable Abel form. In particular, as established in Section 3, the “splitting” technique proved excellent in manipulating and transforming strongly nonlinear ODEs to integrable equations, and hence it may be considered representative of the general pattern of the analysis. Thus, we believe that the developed methodology, possibly modified, extended and enriched with more analytical techniques, can be a powerful tool of research on nonlinear problems in mechanics and physics.

#### Appendices

#### A. Lagrange Method for Quasilinear PDEs of First Order

According to this method, a general solution of the quasilinear equation has the form where with being constants, are solutions of the subsidiary Lagrange equations and is an arbitrary function possessing continuous partial derivatives with respect to its arguments.

#### B. Analytical Parametric Solution of the Equation

It is well known that the general ODE of the first order can accept a parametric solution of the form in case where the following system can be integrated, namely,where the notation , has been adopted. The above system is obtained by the substitution of and differentiation of (B.1) with respect to . In particular, if can be eliminated from (B.2), then a closed-form solution of (B.1) is extracted.

Therefore, as far as the Abel equation is concerned, since it is solvable for , that is, then (B.3b) is considered, namely, () Integration of (B.5) in combination with (B.4) results in with being an arbitrary constant. Moreover, by substituting and taking into account [19, Integral 2.175.1], the parametric solution (B.6) takes the form

#### C. Approximations due to the Weakness of the Field

The weakness of the field under consideration, especially of the coordinate, that is, , allows us to establish the following approximations. (1)We linearly approximate all the polynomials ( represents ) of degree greater or equal than two, namely, (2)Considering the ratio of binomials we evaluate and therefore we obtain

#### D. Expressions for and

In this appendix, we extract appropriate formulas for the function , appearing in (3.2), as well as for the function , involved into the reduction procedure of the Abel equation (3.18) [8, Chapter 1, Section 3.4]. Thus, by considering the function , given from (3.2) and substituting from (2.21), by means of [9, Expression 2.175.1], we arrive at The coefficients involved in various expressions appearing in this appendix are listed in Appendix E. We mention that all these coefficients (appeared through the analytical procedure in this work) are functions of the physical parameter(s), of the problem. Therefore, for this (these) parameter(s) taking values such that the discriminant of becomes positive () if , represent the roots of ( and considering the following cases:

*Case P-1*

*Case P-2*

with , then by elementary algebra (using [9, Expression 2.172]) we can easily prove the following Lemma.

Lemma D.1. *If Case P-1 or P-2 is valid, then the integral can be written in the form:
**
with
*

The coefficients , are different as regards these two cases, while , are common (see Appendix E). Thus, substituting (D.4) into (D.1), we obtain Then by writing in the form and developing in power series (assuming that ) up to the first order, we take Furthermore, by applying (C.4) to (D.5), we arrive at Then developing and up to the second order and substituting (D.8), we conclude that Finally, substitution of (D.7) and (D.9) into (D.6) results in

#### E. List of Coefficients

#### References

- W. F. Ames,
*Nonlinear Partial Differential Equations in Engineering*, Academic Press, New York, NY, USA, 1965. - S. Golstein,
*Modern Developments in Fluid Dynamics*, vol. 1, Oxford University Press, London, UK, 1938. - S. Pai,
*Magnetogas-Dynamics and Plasma Dynamics*, Springer, Vienna, Austria, 1962. - S. Chandrasekhar,
*Plasma Physics*, University of Chicago Press, Chicago, Ill, USA, 1960. - D. E. Panayotounakos and M. Markakis, “Ad hoc closed form solutions of the two-dimensional non-linear steady small perturbation equation in fluid mechanics,”
*International Journal of Non-Linear Mechanics*, vol. 30, no. 4, pp. 597–608, 1995. View at Publisher · View at Google Scholar · View at MathSciNet - H. W. Liepmann and A. Roshko,
*Elements of Gasdynamics*, Galcit Aeronautical Series, John Wiley & Sons, New York, NY, USA, 1957. - S. I. Pai, “One-dimensional unsteady flow of magneto-gasdynamics,” in
*Proceedings of the 5th Midwestern Conference on Fluid Mechanics*, pp. 251–261, University of Michigan Press, Ann Arbor, Mich, USA, 1957. - A. D. Polyanin and V. F. Zaitsev,
*Handbook of Exact Solutions for Ordinary Differential Equations*, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2nd edition, 2003. - I. S. Gradshteyn and I. M. Ryzhik,
*Table of Integrals, Series, and Products*, Academic Press, Boston, Mass, USA, 5th edition, 1994.